From 4e097e1dd7800017621d53788418deccac34c8da Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 10 Jan 2018 16:33:42 +0100 Subject: [PATCH 1/9] support any projection This is work in progress. --- nik4.py | 48 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 41 insertions(+), 7 deletions(-) diff --git a/nik4.py b/nik4.py index 79295c5..0963edd 100755 --- a/nik4.py +++ b/nik4.py @@ -18,16 +18,18 @@ TILE_BUFFER = 128 IM_MONTAGE = 'montage' -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_target = 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') +proj_lonlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') +proj_web_merc = 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') +transform = mapnik.ProjTransform(proj_lonlat, proj_target) +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) + lbbox = layer.envelope().inverse(p).forward(proj_target) if bbox: bbox.expand_to_include(lbbox) else: @@ -178,6 +180,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 +208,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 EPSG:1234 or Proj4 string', default='+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') + 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') @@ -231,6 +242,13 @@ def add_fonts(path): fmt = 'png256' need_cairo = fmt in ['svg', 'pdf'] + + # output projection + if options.projection.lower().startswith('epsg:'): + proj_target = mapnik.Projection('+init={}'.format(options.projection.lower())) + else: + proj_target = mapnik.Projection(options.projection.lower()) + transform = mapnik.ProjTransform(proj_lonlat, proj_target) # get image size in millimeters if options.paper: @@ -278,6 +296,9 @@ 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: @@ -285,6 +306,7 @@ def add_fonts(path): 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: @@ -292,14 +314,24 @@ def add_fonts(path): 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.inverse(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 +353,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 +366,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.inverse(layer_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: From 54937d4f326a8a5fbd4ce91939680596ccfec1b5 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Thu, 11 Jan 2018 12:15:27 +0100 Subject: [PATCH 2/9] use correct methods for transformation --- nik4.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nik4.py b/nik4.py index 0963edd..b111b06 100755 --- a/nik4.py +++ b/nik4.py @@ -328,7 +328,7 @@ def correct_scale(bbox_web_merc, bbox_target): 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.inverse(bbox_web_merc) + bbox = transform_lonlat_webmerc.backward(bbox_web_merc) bbox = transform.forward(bbox) # now correct the scale scale = correct_scale(bbox_web_merc, bbox) @@ -366,7 +366,7 @@ def correct_scale(bbox_web_merc, bbox_target): # 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.inverse(layer_bbox)) + bbox_web_merc = transform_lonlat_webmerc.forward(transform.backward(layer_bbox)) scale = correct_scale(bbox_web_merc, bbox) # expand bbox with padding in mm if bbox and options.padding and (scale or size): From bb793cd883a8f27befe2fdbd6a7ece7ab41af313 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Mon, 15 Jan 2018 13:36:27 +0100 Subject: [PATCH 3/9] use ProjTransform, not Projection for transformations Projection does only project geographical coordinates but ProjTransform takes also care of datum shifts --- nik4.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nik4.py b/nik4.py index b111b06..5b8c844 100755 --- a/nik4.py +++ b/nik4.py @@ -28,8 +28,9 @@ 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(proj_target) + layer_proj = mapnik.Projection(layer.srs) + box_trans = mapnik.ProjTransform(layer_proj, proj_target) + lbbox = box_trans.forward(mapnik.Box2d.from_string('995585 6291591 998380 6293955')) if bbox: bbox.expand_to_include(lbbox) else: @@ -366,8 +367,8 @@ def correct_scale(bbox_web_merc, bbox_target): # 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(layer_bbox)) - scale = correct_scale(bbox_web_merc, bbox) + 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: From ebff0c4966d23ccd2f7f2ab4e7eee02d0b7c57f2 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Mon, 15 Jan 2018 17:32:45 +0100 Subject: [PATCH 4/9] disable ozi map output if --projection is used --- nik4.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nik4.py b/nik4.py index 5b8c844..ebece32 100755 --- a/nik4.py +++ b/nik4.py @@ -231,6 +231,9 @@ def correct_scale(bbox_web_merc, bbox_target): bbox = None rotate = not options.norotate + if options.ozi and options.projection: + 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) From 29f7e64736b70390bb309c208bb84474752d8bde Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 14 Feb 2018 10:38:32 +0100 Subject: [PATCH 5/9] remove forgotten debugging stuff --- nik4.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nik4.py b/nik4.py index ebece32..77cdda4 100755 --- a/nik4.py +++ b/nik4.py @@ -30,7 +30,7 @@ def layer_bbox(m, names, bbox=None): # it may as well be a GPX layer in WGS84 layer_proj = mapnik.Projection(layer.srs) box_trans = mapnik.ProjTransform(layer_proj, proj_target) - lbbox = box_trans.forward(mapnik.Box2d.from_string('995585 6291591 998380 6293955')) + lbbox = box_trans.forward(layer.envelope()) if bbox: bbox.expand_to_include(lbbox) else: From 00f4b0b5f8694d714866ee3600d349c07bb8e6a8 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 14 Feb 2018 10:40:29 +0100 Subject: [PATCH 6/9] remove repeated Proj4 strings by variables --- nik4.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/nik4.py b/nik4.py index 77cdda4..e630733 100755 --- a/nik4.py +++ b/nik4.py @@ -17,10 +17,12 @@ 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' -proj_target = 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') -proj_lonlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') -proj_web_merc = 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') +proj_target = mapnik.Projection(EPSG_3857) +proj_lonlat = mapnik.Projection(EPSG_4326) +proj_web_merc = mapnik.Projection(EPSG_3857) transform = mapnik.ProjTransform(proj_lonlat, proj_target) transform_lonlat_webmerc = mapnik.ProjTransform(proj_lonlat, proj_web_merc) @@ -209,7 +211,7 @@ def correct_scale(bbox_web_merc, bbox_target): 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 EPSG:1234 or Proj4 string', default='+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') + parser.add_argument('-P', '--projection', help='EPSG code as EPSG:1234 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') From baad3f711fa90d26a9585add6ebe3267427f3cc9 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 14 Feb 2018 10:43:05 +0100 Subject: [PATCH 7/9] --projection has a default value and will therefore not be None --- nik4.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nik4.py b/nik4.py index e630733..f40f94b 100755 --- a/nik4.py +++ b/nik4.py @@ -233,7 +233,7 @@ def correct_scale(bbox_web_merc, bbox_target): bbox = None rotate = not options.norotate - if options.ozi and options.projection: + 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: From 30e50f81f2d3a6fdbc730880df5de22d6db6a544 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 14 Feb 2018 15:20:52 +0100 Subject: [PATCH 8/9] remove unnecessary early definition of variables --- nik4.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nik4.py b/nik4.py index f40f94b..31b08a9 100755 --- a/nik4.py +++ b/nik4.py @@ -20,10 +20,8 @@ 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' -proj_target = mapnik.Projection(EPSG_3857) proj_lonlat = mapnik.Projection(EPSG_4326) proj_web_merc = mapnik.Projection(EPSG_3857) -transform = mapnik.ProjTransform(proj_lonlat, proj_target) transform_lonlat_webmerc = mapnik.ProjTransform(proj_lonlat, proj_web_merc) def layer_bbox(m, names, bbox=None): From da68a76f6fba7d4535ef8d2a3adb94379f339916 Mon Sep 17 00:00:00 2001 From: Michael Reichert Date: Wed, 14 Feb 2018 16:12:12 +0100 Subject: [PATCH 9/9] don't expect a prefix EPSG: before a EPSG code --- nik4.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/nik4.py b/nik4.py index 31b08a9..aaaf2ec 100755 --- a/nik4.py +++ b/nik4.py @@ -209,7 +209,7 @@ def correct_scale(bbox_web_merc, bbox_target): 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 EPSG:1234 or Proj4 string', default=EPSG_3857) + 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') @@ -248,12 +248,12 @@ def correct_scale(bbox_web_merc, bbox_target): need_cairo = fmt in ['svg', 'pdf'] # output projection - if options.projection.lower().startswith('epsg:'): - proj_target = mapnik.Projection('+init={}'.format(options.projection.lower())) + if options.projection.isdigit(): + proj_target = mapnik.Projection('+init=epsg:{}'.format(options.projection)) else: - proj_target = mapnik.Projection(options.projection.lower()) + proj_target = mapnik.Projection(options.projection) transform = mapnik.ProjTransform(proj_lonlat, proj_target) - + # get image size in millimeters if options.paper: portrait = False