path: root/macros/latex/contrib/mercatormap/mercatorpy.def

\begin{filecontents*}[overwrite]{maptiles.texpy}
"""
The LaTeX package mercatormap - version 1.02 (2020/08/06)
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.02 (2020/08/06)'


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*}