diff options
author | Norbert Preining <norbert@preining.info> | 2020-04-21 03:02:38 +0000 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2020-04-21 03:02:38 +0000 |
commit | 8731cc5ad70723b9a155917df509d857579c355d (patch) | |
tree | f6594e2b25c3e0aaf4838d7ada0e8885af3fd532 /macros/latex/contrib/mercatormap/mercatorpy.def | |
parent | 2a6dafa8ee740f345650d6a697665bf85a57cc77 (diff) |
CTAN sync 202004210302
Diffstat (limited to 'macros/latex/contrib/mercatormap/mercatorpy.def')
-rw-r--r-- | macros/latex/contrib/mercatormap/mercatorpy.def | 535 |
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*} + |