Merge pull request #29 from Nakaner/any-projection

Support any projection

nik4.py

@@ -17,17 +17,20 @@
 VERSION = '1.6'
 TILE_BUFFER = 128
 IM_MONTAGE = 'montage'
+EPSG_4326 = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
+EPSG_3857 = '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over'
 
-p3857 = mapnik.Projection('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over')
-p4326 = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
-transform = mapnik.ProjTransform(p4326, p3857)
+proj_lonlat = mapnik.Projection(EPSG_4326)
+proj_web_merc = mapnik.Projection(EPSG_3857)
+transform_lonlat_webmerc = mapnik.ProjTransform(proj_lonlat, proj_web_merc)
 
 def layer_bbox(m, names, bbox=None):
 	"""Calculate extent of given layers and bbox"""
 	for layer in (l for l in m.layers if l.name in names):
 		# it may as well be a GPX layer in WGS84
-		p = mapnik.Projection(layer.srs)
-		lbbox = layer.envelope().inverse(p).forward(p3857)
+		layer_proj = mapnik.Projection(layer.srs)
+		box_trans = mapnik.ProjTransform(layer_proj, proj_target)
+		lbbox = box_trans.forward(layer.envelope())
 		if bbox:
 			bbox.expand_to_include(lbbox)
 		else:
@@ -178,6 +181,13 @@ def add_fonts(path):
 	else:
 		raise Exception('The directory "{p}" does not exists'.format(p=path))
 
+def correct_scale(bbox_web_merc, bbox_target):
+	# correct scale if output projection is not EPSG:3857
+	x_dist_merc = bbox_web_merc.maxx - bbox_web_merc.minx
+	x_dist_target = bbox.maxx - bbox.minx
+	return scale * (x_dist_target / x_dist_merc)
+
+
 if __name__ == "__main__":
 	parser = argparse.ArgumentParser(description='Nik4 {}: Tile-aware mapnik image renderer'.format(VERSION))
 	parser.add_argument('--version', action='version', version='Nik4 {}'.format(VERSION))
@@ -199,6 +209,8 @@ def add_fonts(path):
 	parser.add_argument('--add-layers', help='Map layers to include, comma-separated')
 	parser.add_argument('--hide-layers', help='Map layers to hide, comma-separated')
 
+        parser.add_argument('-P', '--projection', help='EPSG code as 1234 (without prefix "EPSG:" or Proj4 string', default=EPSG_3857)
+
 	parser.add_argument('--url', help='URL of a map to center on')
 	parser.add_argument('--ozi', type=argparse.FileType('w'), help='Generate ozi map file')
 	parser.add_argument('--wld', type=argparse.FileType('w'), help='Generate world file')
@@ -219,6 +231,9 @@ def add_fonts(path):
 	bbox = None
 	rotate = not options.norotate
 
+        if options.ozi and options.projection.lower() != 'epsg:3857' and options.projection != EPSG_3857:
+                raise Exception('ozi map file output is only supported for Web Mercator (EPSG:3857). Please remove --projection.')
+
 	if options.url:
 		parse_url(options.url, options)
 
@@ -231,7 +246,14 @@ def add_fonts(path):
 		fmt = 'png256'
 
 	need_cairo = fmt in ['svg', 'pdf']
-
+
+	# output projection
+	if options.projection.isdigit():
+		proj_target = mapnik.Projection('+init=epsg:{}'.format(options.projection))
+	else:
+		proj_target = mapnik.Projection(options.projection)
+	transform = mapnik.ProjTransform(proj_lonlat, proj_target)
+
 	# get image size in millimeters
 	if options.paper:
 		portrait = False
@@ -278,28 +300,42 @@ def add_fonts(path):
 	if size and size[0] + size[1] <= 0:
 		raise Exception('Both dimensions are less or equal to zero')
 
+	if options.bbox:
+		bbox = options.bbox
+
 	# scale can be specified with zoom or with 1:NNN scale
 	fix_scale = False
 	if options.zoom:
 		scale = 2 * 3.14159 * 6378137 / 2 ** (options.zoom + 8) / scale_factor
 	elif options.scale:
 		scale = options.scale * 0.00028 / scale_factor
 		# Now we have to divide by cos(lat), but we might not know latitude at this point
+		#TODO division should only happen for EPSG:3857 or not at all
 		if options.center:
 			scale = scale / math.cos(math.radians(options.center[1]))
 		elif options.bbox:
 			scale = scale / math.cos(math.radians((options.bbox[3] + options.bbox[1]) / 2))
 		else:
 			fix_scale = True
 
-	if options.bbox:
-		bbox = options.bbox
 	# all calculations are in EPSG:3857 projection (it's easier)
 	if bbox:
 		bbox = transform.forward(mapnik.Box2d(*bbox))
+		bbox_web_merc = transform_lonlat_webmerc.forward(mapnik.Box2d(*(options.bbox)))
+		scale = correct_scale(bbox_web_merc, bbox)
 
 	# calculate bbox through center, zoom and target size
 	if not bbox and options.center and size and size[0] > 0 and size[1] > 0 and scale:
+		# We don't know over which latitude range the bounding box spans, so we
+		# first do everything in Web Mercator.
+		center = transform_lonlat_webmerc.forward(mapnik.Coord(*options.center))
+		w = size[0] * scale / 2
+		h = size[1] * scale / 2
+		bbox_web_merc = mapnik.Box2d(center.x-w, center.y-h, center.x+w, center.y+h)
+		bbox = transform_lonlat_webmerc.backward(bbox_web_merc)
+		bbox = transform.forward(bbox)
+		# now correct the scale
+		scale = correct_scale(bbox_web_merc, bbox)
 		center = transform.forward(mapnik.Coord(*options.center))
 		w = size[0] * scale / 2
 		h = size[1] * scale / 2
@@ -321,7 +357,7 @@ def add_fonts(path):
 	# for layer processing we need to create the Map object
 	m = mapnik.Map(100, 100) # temporary size, will be changed before output
 	mapnik.load_map_from_string(m, style_xml, False, style_path)
-	m.srs = p3857.params()
+	m.srs = proj_target.params()
 
 	# register non-standard fonts
 	if options.fonts:
@@ -334,6 +370,8 @@ def add_fonts(path):
 		# here's where we can fix scale, no new bboxes below
 		if bbox and fix_scale:
 			scale = scale / math.cos(math.radians(transform.backward(bbox.center()).y))
+		bbox_web_merc = transform_lonlat_webmerc.forward(transform.backward(bbox))
+		scale = correct_scale(bbox_web_merc, bbox)
 		# expand bbox with padding in mm
 		if bbox and options.padding and (scale or size):
 			if scale: