WIP: add support for UTM EPSGs

[Scartography]

• this enables all 120 UTM stripes over the globe in tmx

[Scartography]

Based and integrated on this simple logic:

import os
import geojson

from tilematrix._grid import GridDefinition
from tilematrix._tilepyramid import TilePyramid

start_utm_stripe = 32600

UTM_DEFAULT_DICT = {
    '32600': {
        'grid': None,
        'shape': [8, 1],
        'bounds': [166021.4431, 0.0000, 1476741.4431, 10485760],
        'crs': {"epsg": 32600},
        'is_global': False
    }
}

# CRS_GEOJSON_DICT = {
#     "crs": {
#         "type": "name",
#         "properties": {
#         "name": "EPSG:32601"
#         }
#     }
# }

OUT_UTM_DICT = UTM_DEFAULT_DICT
OUT_UTM_GRIDS = {}
OUT_UTM_PYRAMIDS = {}


def get_utm_tmx_grids(
    utm_start_col=1,
    utm_end_col=60,
    tile_size=256,
    metatiling=16,
):
    for i in range(utm_start_col, utm_end_col+1):
        utm_stripe = start_utm_stripe + i
        OUT_UTM_DICT[utm_stripe] = {
                    'shape': [8, 1],
                    'bounds': [166020.0, -50000.0, 1476740, 10435760.0],
                    'crs': {"epsg": utm_stripe},
                    'is_global': False
                }
        OUT_UTM_GRIDS[utm_stripe] = GridDefinition(
            grid=None,
            shape=OUT_UTM_DICT[utm_stripe]['shape'],
            bounds=OUT_UTM_DICT[utm_stripe]['bounds'],
            srs=OUT_UTM_DICT[utm_stripe]['crs'],
            is_global=OUT_UTM_DICT[utm_stripe]['is_global']
        )
        OUT_UTM_PYRAMIDS[utm_stripe] = TilePyramid(
            grid=OUT_UTM_GRIDS[utm_stripe],
            tile_size=tile_size,
            metatiling=metatiling,
        )

    return OUT_UTM_DICT, OUT_UTM_GRIDS, OUT_UTM_PYRAMIDS


def get_utm_tmx_pyramid_geojson(
        crs_epsg,
        pyramid,
        zoom=0,
        bounds=(166020.0, -50000.0, 1476740, 10435760.0),
        pixel_buffer=0,
        clip_to_utm_bouds=True,
):
    if clip_to_utm_bouds:
        bounds = (166021.4431, 0.0000, 833978.5569, 9329005.18)
    out_dict = {}
    tiles = []
    for i, tile in enumerate(list(
            pyramid.tiles_from_bounds(bounds=tuple(bounds), zoom=zoom))
    ):
        out_dict.update(
            id=i,
            geometry=tile.bbox(),
            properties={}
        )
        tiles.append(geojson.Feature(
                    geometry=tile.bbox(pixelbuffer=pixel_buffer),
                    properties=dict(
                        zoom=tile.zoom,
                        row=tile.row,
                        col=tile.col
                    )
                )
        )
    # Update geojson with its CRS
    out_tiles = geojson.FeatureCollection(tiles)

    out_tiles["crs"] = {
        "type": "name",
        "properties": {
            "name": f"EPSG:{crs_epsg}"
        }
    }
    return out_tiles


# This will take all UTM North Zones and write the tilepyramid
# output to geojsons

# By default the output is clipped to the UTM bounds as per (left+top bounds):
# https://spatialreference.org/ref/epsg/32601/
for nr, pyramid in enumerate(get_utm_tmx_grids()[2], start=1):
    tiles = get_utm_tmx_pyramid_geojson(
        crs_epsg=pyramid,
        pyramid=get_utm_tmx_grids()[2][pyramid],
        zoom=2
    )
    if os.path.isfile('out_stuff/'+str(pyramid)+'.geojson'):
        os.remove('out_stuff/'+str(pyramid)+'.geojson')
    with open('out_stuff/'+str(pyramid)+'.geojson', 'w') as outfile:
        geojson.dump(tiles, outfile)
    print('DONE: {} of {}'.format(nr, len(get_utm_tmx_grids()[2].keys())))

[Scartography]

Use OGC naming for UTM grid name

[Scartography]

Figure out how to work with southern hemisphere

[Scartography]

This TMX UTM Grid Size is defined to match 10[m] pixel_size/resolution at zoom 9!!! Which is the OGC definition: https://docs.opengeospatial.org/is/17-083r2/17-083r2.html#65

The 10[m] is to aim at sensible the Sentinel-2 regrinding and targeting exactly the projected 10 meters.

The top left origin of the southern UTM zones and the bottom left origin of the northern UTM is the one of the EPSG bounds for those to have an exact grid at equator, this is not OGC compatibel:

• https://epsg.io/32633
• https://epsg.io/32733
  It is subject to be done adaptively if we decide to keep pyproj as dependency of the tmx!

# Grid with width (x-dif) of 1310720 and height (y-dif) of 10485760
# This leads to exactly 10[m] grid at zoom 9
UTM_STRIPE_SHAPE = (8, 1)
UTM_STRIPE_NORTH_BOUNDS = [166021.4431, 0.0000, 1476741.4431, 10485760]
UTM_STRIPE_SOUTH_BOUNDS = [441867.78, -485760.00, 1752587.78, 10000000.00]