\begin{filecontents}[force,noheader,nowarn]{maptiles.texpy} """ The LaTeX package mercatormap - version 1.2.0 (2024/08/05) maptiles.texpy: Python script for map tile download ------------------------------------------------------------------------------------------- Copyright (c) 2020-2024 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 https://www.latex-project.org/lppl.txt and version 1.3c or later is part of all distributions of LaTeX version 2008-05-04 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.2.0 (2024/08/05)' 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}