summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/mercatormap/mercatorpy.def
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2020-04-21 03:02:38 +0000
committerNorbert Preining <norbert@preining.info>2020-04-21 03:02:38 +0000
commit8731cc5ad70723b9a155917df509d857579c355d (patch)
treef6594e2b25c3e0aaf4838d7ada0e8885af3fd532 /macros/latex/contrib/mercatormap/mercatorpy.def
parent2a6dafa8ee740f345650d6a697665bf85a57cc77 (diff)
CTAN sync 202004210302
Diffstat (limited to 'macros/latex/contrib/mercatormap/mercatorpy.def')
-rw-r--r--macros/latex/contrib/mercatormap/mercatorpy.def535
1 files changed, 535 insertions, 0 deletions
diff --git a/macros/latex/contrib/mercatormap/mercatorpy.def b/macros/latex/contrib/mercatormap/mercatorpy.def
new file mode 100644
index 0000000000..fe608f47b4
--- /dev/null
+++ b/macros/latex/contrib/mercatormap/mercatorpy.def
@@ -0,0 +1,535 @@
+\begin{filecontents*}[overwrite]{maptiles.texpy}
+"""
+The LaTeX package mercatormap - version 1.00 (2020/04/20)
+maptiles.texpy: Python script for map tile download
+
+-------------------------------------------------------------------------------------------
+Copyright (c) 2020 by Prof. Dr. Dr. Thomas F. Sturm <thomas dot sturm at unibw dot de>
+-------------------------------------------------------------------------------------------
+
+This work may be distributed and/or modified under the
+conditions of the LaTeX Project Public License, either version 1.3
+of this license or (at your option) any later version.
+The latest version of this license is in
+ http://www.latex-project.org/lppl.txt
+and version 1.3 or later is part of all distributions of LaTeX
+version 2005/12/01 or later.
+
+This work has the LPPL maintenance status `author-maintained'.
+
+This work consists of all files listed in README
+"""
+
+import argparse
+import math
+import os
+import requests
+from pathlib import Path
+from PIL import Image
+
+
+packageversion = '1.00 (2020/04/20)'
+
+
+def gd(x):
+ return math.atan(math.sinh(x))
+
+
+def arggd(x):
+ return math.asinh(math.tan(x))
+
+
+def rad2planar(lon_rad, lat_rad):
+ x = 0.5*lon_rad/math.pi + 0.5
+ y = 0.5*arggd(lat_rad)/math.pi + 0.5
+ return (x, y)
+
+
+def deg2planar(lon_deg, lat_deg):
+ return rad2planar(math.radians(lon_deg), math.radians(lat_deg))
+
+
+def planar2rad(x, y):
+ lon_rad = math.pi*(2*x-1)
+ lat_rad = gd(math.pi*(2*y-1))
+ return (lon_rad, lat_rad)
+
+
+def planar2deg(x, y):
+ (lon_rad, lat_rad) = planar2rad(x, y)
+ return (math.degrees(lon_rad), math.degrees(lat_rad))
+
+
+def download(session, url, filename, overwrite=False):
+ local_file = Path(filename)
+ if local_file.is_file() and not overwrite:
+ return
+ print(url+' -> '+filename)
+ r = session.get(url)
+ r.raise_for_status()
+ with open(filename, 'wb') as handler:
+ handler.write(r.content)
+ return
+
+
+def makeparentdirs(afile):
+ directory = Path(os.path.abspath(afile)).parent
+ if not directory.is_dir():
+ os.makedirs(directory)
+
+
+class MapDefinition(object):
+ def __init__(self, args, west, east, north, south, zoom, x_1, y_1, x_2, y_2,
+ width, height, west_offset, north_offset, south_offset):
+ self.west = west
+ self.east = east
+ self.north = north
+ self.south = south
+ self.zoom = zoom
+ self.x_1 = x_1
+ self.y_1 = y_1
+ self.x_2 = x_2
+ self.y_2 = y_2
+ self.width = width
+ self.height= height
+ self.west_offset = west_offset
+ self.north_offset = north_offset
+ self.south_offset = south_offset
+ self.n = 2 ** self.zoom
+ self.local_file_prefix = 'tiles/tile_'
+ self.attribution = args.attribution
+ self.attributionprint = args.attributionprint
+ self.pixel = args.pixel
+ self.pixel_width = round(self.width*self.pixel)
+ self.pixel_height = round(self.height*self.pixel)
+ self.tilesize = args.tilesize
+
+
+ @classmethod
+ def construct_from_boundaries_rad(cls, args, zoom, west_rad, east_rad, north_rad, south_rad):
+ if zoom<0:
+ zoom = 0
+ if west_rad>=east_rad:
+ east_rad += 2*math.pi*math.ceil(0.5*(west_rad-east_rad)/math.pi+1e-10)
+ max_lat = gd(math.pi)-1e-10
+ if north_rad<south_rad:
+ (north_rad, south_rad) = (south_rad, north_rad)
+ if north_rad > max_lat:
+ north_rad = max_lat
+ if south_rad < -max_lat:
+ south_rad = -max_lat
+ P = 2**(zoom-1)/math.pi
+ a1 = (west_rad+math.pi)*P
+ b1 = (math.pi-arggd(north_rad))*P
+ a2 = (east_rad+math.pi)*P
+ b2 = (math.pi-arggd(south_rad))*P
+ x1 = math.floor(a1)
+ y1 = math.floor(b1)
+ x2 = math.ceil(a2)-1
+ y2 = math.ceil(b2)-1
+ width = a2-a1
+ height = b2-b1
+ west_offset = a1-x1
+ north_offset = b1-y1
+ south_offset = y2-b2+1
+ return cls(args, west_rad, east_rad, north_rad, south_rad, zoom, x1, y1, x2, y2,
+ width, height, west_offset, north_offset, south_offset )
+
+
+ @classmethod
+ def construct_from_boundaries_deg(cls, args, zoom, west_deg, east_deg, north_deg, south_deg):
+ return cls.construct_from_boundaries_rad(args,
+ zoom,math.radians(west_deg),math.radians(east_deg),
+ math.radians(north_deg),math.radians(south_deg))
+
+
+ @classmethod
+ def construct_from_boundaries_args(cls, args):
+ return cls.construct_from_boundaries_rad(args,
+ args.zoom,math.radians(args.west),math.radians(args.east),math.radians(args.north),math.radians(args.south))
+
+
+ @classmethod
+ def construct_from_reference_rad(cls, args, zoom, width, height,
+ ref_lon_rad, ref_lat_rad, xalign, yalign):
+ if zoom<0:
+ zoom = 0
+ P = 2**(zoom-1)/math.pi
+ ac = (ref_lon_rad+math.pi)*P
+ bc = (math.pi-arggd(ref_lat_rad))*P
+ a1 = ac - (1+xalign)/2*width
+ a2 = a1 + width
+ b1 = bc - (1+yalign)/2*height
+ b2 = b1 + height
+ west_rad = a1/P - math.pi
+ north_rad = gd(math.pi-b1/P)
+ east_rad = a2/P - math.pi
+ south_rad = gd(math.pi-b2/P)
+ x1 = math.floor(a1)
+ y1 = math.floor(b1)
+ x2 = math.ceil(a2)-1
+ y2 = math.ceil(b2)-1
+ west_offset = a1-x1
+ north_offset = b1-y1
+ south_offset = y2-b2+1
+ return cls(args, west_rad, east_rad, north_rad, south_rad, zoom, x1, y1, x2, y2,
+ width, height, west_offset, north_offset, south_offset )
+
+
+ @classmethod
+ def construct_from_reference_deg(cls, args, zoom, width, height,
+ ref_lon_deg, ref_lat_deg, xalign, yalign):
+ return cls.construct_from_center_rad(args, zoom, width, height,
+ math.radians(ref_lon_deg), math.radians(ref_lat_deg), xalign, yalign)
+
+
+ @classmethod
+ def construct_from_reference_args(cls, args):
+ return cls.construct_from_reference_rad(args, args.zoom, args.width, args.height,
+ math.radians(args.longitude), math.radians(args.latitude), args.xalign, args.yalign)
+
+
+ @classmethod
+ def construct_to_fit_region_rad(cls, args, width, height, reg_west_rad, reg_east_rad,
+ reg_north_rad, reg_south_rad, xalign, yalign):
+ if reg_west_rad>=reg_east_rad:
+ reg_east_rad += 2*math.pi*math.ceil(0.5*(reg_west_rad-reg_east_rad)/math.pi+1e-10)
+ max_lat = gd(math.pi)-1e-10
+ if reg_north_rad<reg_south_rad:
+ (reg_north_rad, reg_south_rad) = (reg_south_rad, reg_north_rad)
+ if reg_north_rad > max_lat:
+ reg_north_rad = max_lat
+ if reg_south_rad < -max_lat:
+ reg_south_rad = -max_lat
+ zoom = max( 0, 1 + math.floor( math.log2( math.pi *
+ min( width/(reg_east_rad-reg_west_rad), height/(arggd(reg_north_rad)-arggd(reg_south_rad)))
+ )))
+ P = 2**(zoom-1)/math.pi
+
+ if xalign<0:
+ a1 = (reg_west_rad+math.pi)*P
+ a2 = a1 + width
+ elif xalign>0:
+ a2 = (reg_east_rad+math.pi)*P
+ a1 = a2 - width
+ else:
+ a1 = ((reg_west_rad+reg_east_rad)/2+math.pi)*P - width/2
+ a2 = a1 + width
+
+ if yalign<0:
+ b1 = (math.pi-arggd(reg_north_rad))*P
+ b2 = b1 + height
+ elif yalign>0:
+ b2 = (math.pi-arggd(reg_south_rad))*P
+ b1 = b2 - height
+ else:
+ b1 = (math.pi-(arggd(reg_north_rad)+arggd(reg_south_rad))/2)*P - height/2
+ b2 = b1 + height
+ west_rad = a1/P - math.pi
+ north_rad = gd(math.pi-b1/P)
+ east_rad = a2/P - math.pi
+ south_rad = gd(math.pi-b2/P)
+ x1 = math.floor(a1)
+ y1 = math.floor(b1)
+ x2 = math.ceil(a2)-1
+ y2 = math.ceil(b2)-1
+ west_offset = a1-x1
+ north_offset = b1-y1
+ south_offset = y2-b2+1
+ return cls(args, west_rad, east_rad, north_rad, south_rad, zoom, x1, y1, x2, y2,
+ width, height, west_offset, north_offset, south_offset )
+
+
+ @classmethod
+ def construct_to_fit_region_deg(cls, args, width, height,
+ reg_west_deg, reg_east_deg, reg_north_deg, reg_south_deg,
+ xalign, yalign):
+ return cls.construct_to_fit_region_rad(args, width, height,
+ math.radians(reg_west_deg),math.radians(reg_east_deg),
+ math.radians(reg_north_deg),math.radians(reg_south_deg),
+ xalign, yalign)
+
+
+ @classmethod
+ def construct_to_fit_region_args(cls, args):
+ return cls.construct_to_fit_region_rad(args, args.width, args.height,
+ math.radians(args.west),math.radians(args.east),
+ math.radians(args.north),math.radians(args.south),
+ args.xalign, args.yalign)
+
+
+ def get_tile_quantity(self):
+ return (self.x_2-self.x_1+1)*(self.y_2-self.y_1+1)
+
+
+ def get_pixel_quantity(self):
+ return self.pixel_width*self.pixel_height
+
+
+ def download_tiles(self, url_format, local_file_prefix):
+ self.local_file_prefix = local_file_prefix
+ makeparentdirs(local_file_prefix)
+ with requests.Session() as s:
+ s.headers.update(dict(referer = 'mercatormap/'+packageversion))
+ try:
+ for xx in range(self.x_1, self.x_2+1):
+ x = xx % self.n
+ for y in range(self.y_1, self.y_2+1):
+ if (y>=0) and (y<self.n):
+ url = url_format.format(x=x, y=y, z=self.zoom)
+ filename = local_file_prefix + '_{z}_{x}_{y}.png'.format(x=x, y=y, z=self.zoom)
+ download(s,url,filename)
+ except requests.exceptions.RequestException as err:
+ print('Download error:', err)
+ return False
+ return True
+
+
+ def merge_map(self, map_definition, local_file_prefix):
+ makeparentdirs(map_definition)
+ new_im = Image.new('RGB', (self.pixel_width, self.pixel_height), 'white')
+ xoff = -round(self.west_offset*self.pixel)
+ for xx in range(self.x_1, self.x_2+1):
+ x = xx % self.n
+ yoff = -round(self.north_offset*self.pixel)
+ for y in range(self.y_1, self.y_2+1):
+ if (y>=0) and (y<self.n):
+ filename = local_file_prefix + '_{z}_{x}_{y}.png'.format(x=x, y=y, z=self.zoom)
+ new_im.paste(Image.open(filename), (xoff,yoff))
+ yoff+=self.pixel
+ xoff+=self.pixel
+ new_im.save(map_definition+'.png', optimize=True)
+ return True
+
+
+ def download_wms_map(self, url_format, map_definition):
+ makeparentdirs(map_definition)
+ # bbox in tiles for given zoom
+ xmin = self.x_1+self.west_offset
+ xmax = xmin+self.width
+ ymin = self.y_1+self.north_offset
+ ymax = ymin+self.height
+
+ # global scale bbox
+ f = math.pi*6378137
+ xmin = (2*xmin/self.n-1)*f
+ xmax = (2*xmax/self.n-1)*f
+ (ymin, ymax) = (ymax, ymin) # switch ymin/ymax
+ ymin = (1-2*ymin/self.n)*f
+ ymax = (1-2*ymax/self.n)*f
+
+ url = url_format.format(xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax,
+ width=self.pixel_width, height=self.pixel_height)
+ filename = map_definition+'.png'
+ with requests.Session() as s:
+ s.headers.update(dict(referer = 'mercatormap/'+packageversion))
+ try:
+ download(s,url,filename,True)
+ except requests.exceptions.RequestException as err:
+ print('Download error:', err)
+ return False
+ return True
+
+
+ def save_map(self, map_definition, tiles=False, merge=False, wms=False):
+ makeparentdirs(map_definition)
+ with open(map_definition+'.def', "w", encoding="utf-8") as desc:
+ desc.write(f'% generated by mercatormap maptiles {packageversion}\n')
+ desc.write('\\mrcdefinemap{\n')
+ desc.write(f' west={math.degrees(self.west)},\n')
+ desc.write(f' east={math.degrees(self.east)},\n')
+ desc.write(f' north={math.degrees(self.north)},\n')
+ desc.write(f' south={math.degrees(self.south)},\n')
+ desc.write(f' zoom={self.zoom},\n')
+ desc.write(f' xmin={self.x_1},\n')
+ desc.write(f' ymin={self.y_1},\n')
+ desc.write(f' xmax={self.x_2},\n')
+ desc.write(f' ymax={self.y_2},\n')
+ desc.write(f' pixelwidth={self.pixel_width},\n')
+ desc.write(f' pixelheight={self.pixel_height},\n')
+ desc.write(f' westoffset={self.west_offset},\n')
+ desc.write(f' northoffset={self.north_offset},\n')
+ desc.write(f' southoffset={self.south_offset},\n')
+ desc.write(f' basename={self.local_file_prefix},\n')
+ if self.attribution!=None:
+ desc.write(' attribution={'+self.attribution+'},\n')
+ if self.attributionprint!=None:
+ desc.write(' attribution print={'+self.attributionprint+'},\n')
+ if self.tilesize!=None:
+ desc.write(f' tile size={self.tilesize},\n')
+ if wms:
+ desc.write(f' resource=wmsmap,\n')
+ elif merge:
+ desc.write(f' resource=mergedmap,\n')
+ elif tiles:
+ desc.write(f' resource=tiles,\n')
+ else:
+ desc.write(f' resource=none,\n')
+ desc.write('}\n')
+ return
+
+
+ def print_description(self):
+ print('------ map description (class) ------')
+ print(f'west = {math.degrees(self.west)}')
+ print(f'east = {math.degrees(self.east)}')
+ print(f'north = {math.degrees(self.north)}')
+ print(f'south = {math.degrees(self.south)}')
+ print(f'zoom = {self.zoom}')
+ print(f'X_1 = {self.x_1}')
+ print(f'Y_1 = {self.y_1}')
+ print(f'X_2 = {self.x_2}')
+ print(f'Y_2 = {self.y_2}')
+ print(f'width = {self.width}')
+ print(f'height = {self.height}')
+ print(f'd_west = {self.west_offset}')
+ print(f'd_north = {self.north_offset}')
+ print(f'd_south = {self.south_offset}')
+ return
+
+
+def map_tile_info(Z, X, Y):
+ N = 2.0 ** Z
+ lon_min = math.pi*(2*X-N)/N
+ lon_max = math.pi*(2*X+2-N)/N
+ lat_min = gd(math.pi*(N-2*Y-2)/N)
+ lat_max = gd(math.pi*(N-2*Y)/N)
+ r = 6371000
+ height = r*(lat_max-lat_min)
+ rdel_lon = r*(lon_max-lon_min)
+ width_n = rdel_lon*math.cos(lat_max)
+ width_s = rdel_lon*math.cos(lat_min)
+ area = r*rdel_lon*(math.sin(lat_max)-math.sin(lat_min))
+ print(f'(Z,X,Y) = ({Z},{X},{Y})')
+ print(f'west = {math.degrees(lon_min)}')
+ print(f'east = {math.degrees(lon_max)}')
+ print(f'south = {math.degrees(lat_min)}')
+ print(f'north = {math.degrees(lat_max)}')
+ print(f'height = {height/1000}')
+ print(f'width(N)= {width_n/1000}')
+ print(f'width(S)= {width_s/1000}')
+ print(f'area = {area/1000000}')
+ return
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Compute a full map description and download map tiles')
+
+ parser.add_argument('type', type=str, help='type of map computation',
+ choices=['boundaries', 'reference', 'areafit'])
+
+ parser.add_argument('-z', '--zoom', type=int, default=9,
+ help='zoom level (boundaries/reference)')
+ parser.add_argument('-w', '--west', type=float, default=11,
+ help='longitude west (boundaries/areafit)')
+ parser.add_argument('-e', '--east', type=float, default=13,
+ help='longitude east (boundaries/areafit)')
+ parser.add_argument('-n', '--north', type=float, default=48,
+ help='latitude north (boundaries/areafit)')
+ parser.add_argument('-s', '--south', type=float, default=50,
+ help='latitude south (boundaries/areafit)')
+ parser.add_argument('-lon', '--longitude', type=float, default=12,
+ help='longitude (reference)')
+ parser.add_argument('-lat', '--latitude', type=float, default=49,
+ help='latitude (reference)')
+ parser.add_argument('-mw', '--width', type=float, default=4,
+ help='map width as multiplicity of map tiles (reference/areafit)')
+ parser.add_argument('-mh', '--height', type=float, default=4,
+ help='map height as multiplicity of map tiles (reference/areafit)')
+ parser.add_argument('-a', '--align', type=str, default='center',
+ help='alignment of reference or area (reference/areafit)',
+ choices=['northwest', 'north', 'northeast',
+ 'west', 'center', 'east',
+ 'southwest', 'south', 'southeast'] )
+ parser.add_argument('-d', '--definition', type=str, default='mapdefinition',
+ help='name of the produced map definition for LaTeX')
+ parser.add_argument('-t', '--target', type=str, default='tiles',
+ help='target to create tiles or a merged map',
+ choices=['none', 'tiles', 'mergedmap', 'wmsmap'] )
+ parser.add_argument('-u', '--url', type=str, default=None,
+ help='url format with placeholder {z}{x}{y} for downloading tiles')
+ parser.add_argument('-at', '--attribution', type=str, default=None,
+ help='attribution text of the tile producer')
+ parser.add_argument('-atp', '--attributionprint', type=str, default=None,
+ help='attribution text of the tile producer')
+ parser.add_argument('-f', '--fileprefix', type=str, default='tiles/unknown',
+ help='prefix for local tile files, e.g. "tiles/map" for "tiles/map_10_10_10.png"')
+ parser.add_argument('-c', '--confirm', action="store_true",
+ help='confirm bulk download for map tiles')
+ parser.add_argument('-p', '--pixel', type=int, default=256,
+ help='pixel width/height of a tile')
+ parser.add_argument('-ts', '--tilesize', type=str, default=None,
+ help='TeX tile size')
+
+ args = parser.parse_args()
+
+ print('')
+ print('** This is maptiles.py version', packageversion,
+ 'for generating', args.definition)
+
+ switcher = {
+ 'northwest': (-1,-1),
+ 'north': ( 0,-1),
+ 'northeast': ( 1,-1),
+ 'west': (-1, 0),
+ 'center': ( 0, 0),
+ 'east': ( 1, 0),
+ 'southwest': (-1, 1),
+ 'south': ( 0, 1),
+ 'southeast': ( 1, 1)
+ }
+ (args.xalign, args.yalign) = switcher.get(args.align)
+
+ switcher = {
+ 'boundaries': MapDefinition.construct_from_boundaries_args,
+ 'reference': MapDefinition.construct_from_reference_args,
+ 'areafit': MapDefinition.construct_to_fit_region_args
+ }
+ mapdef = switcher.get(args.type)(args)
+
+ switcher = {
+ 'none': 0,
+ 'tiles': 1,
+ 'mergedmap': 2,
+ 'wmsmap': 3
+ }
+ target = switcher.get(args.target)
+
+ tiles_success = False
+ merge_success = False
+ wms_success = False
+
+ if target>0:
+ if args.url!=None:
+ if target<3:
+ if mapdef.get_tile_quantity()>400:
+ if args.confirm:
+ tiles_success = mapdef.download_tiles(args.url, args.fileprefix)
+ else:
+ print(f'Too many tiles ({mapdef.get_tile_quantity()}). Add confirm option to proceed anyway.');
+ else:
+ tiles_success = mapdef.download_tiles(args.url, args.fileprefix)
+ if tiles_success:
+ print('** Download of tiles successful for', args.definition)
+ if target>1:
+ merge_success = mapdef.merge_map(args.definition, args.fileprefix)
+ if merge_success:
+ print('** Merging tiles successful for', args.definition)
+ else:
+ print('** WARNING: Merging tiles failed for', args.definition)
+ else:
+ print('** WARNING: Map tiles download failed for', args.definition)
+ else:
+ if mapdef.get_pixel_quantity()>100000000:
+ if args.confirm:
+ wms_success = mapdef.download_wms_map(args.url, args.definition)
+ else:
+ print(f'Too many pixels ({mapdef.get_pixel_quantity()}). Add confirm option to proceed anyway.');
+ else:
+ wms_success = mapdef.download_wms_map(args.url, args.definition)
+ else:
+ print('No url for download given')
+
+ mapdef.save_map(args.definition, tiles_success, merge_success, wms_success)
+\end{filecontents*}
+