WIP: add support for UTM EPSGs

requirements.txt

@@ -1,5 +1,6 @@
 affine
 click
 geojson
+pyproj
 rasterio>=1.0.21
 shapely

tilematrix/_conf.py

@@ -1,11 +1,117 @@
 """Package configuration parameters."""
 
+from pyproj.aoi import AreaOfInterest
+from pyproj.database import query_utm_crs_info
+
+
 # round coordinates
 ROUND = 20
 
 # bounds ratio vs shape ratio uncertainty
 DELTA = 1e-6
 
+# UTM default EOX settings
+FIRST_UTM_STRIPE = 32600
+LAST_UTM_STRIPE = 60
+
+# Grid with width (x-dif) of 1310720 and height (y-dif) of 10485760
+# This leads to exactly 10[m] grid at zoom 9 and shape at z0 8 high and 1 wide
+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]
+
+# Analog logic for S2 grid but the origin is shifted to match the S2 grid
+# and the TileMatrix Definition is defined to be divisible by 60, 20 and 10
+# width 737280 to be divisible by 60, 20 and 10 and has 10[m] resolution
+# height 9584640 to be divisible by 60, 20 and 10 and has 10[m] resolution
+# S2 originally intended pyramid bounds [99960.00, 0.0000, 834000.00, 9000000]
+# For the 60[m] grid the max/optimal zoom level is 4
+# For the 10[m] and 20[m] grid the optimal zoom level is 5 and 4 respectively
+
+# 10[m] and 20[m] grid
+UTM_STRIPE_S2_SHAPE = (117, 9)
+
+# 60[m] grid
+UTM_STRIPE_S2_60M_SHAPE = (39, 3)
+
+UTM_STRIPE_S2_GRID_NORTH_BOUNDS = [145980, 0, 883260, 9584640.0]
+UTM_STRIPE_S2_GRID_SOUTH_BOUNDS = [145980, 415360, 883260, 10000000.0]
+
+# Availible UTM grids in TMX
+UTM_GRIDS = {
+    "EOX_UTM": {
+        "shape": UTM_STRIPE_SHAPE,
+        "north_bounds": UTM_STRIPE_NORTH_BOUNDS,
+        "south_bounds": UTM_STRIPE_SOUTH_BOUNDS
+    },
+    "S2_10M_UTM": {
+        "shape": UTM_STRIPE_S2_SHAPE,
+        "north_bounds": UTM_STRIPE_S2_GRID_NORTH_BOUNDS,
+        "south_bounds": UTM_STRIPE_S2_GRID_SOUTH_BOUNDS
+    },
+    "S2_60M_UTM": {
+        "shape": UTM_STRIPE_S2_60M_SHAPE,
+        "north_bounds": UTM_STRIPE_S2_GRID_NORTH_BOUNDS,
+        "south_bounds": UTM_STRIPE_S2_GRID_SOUTH_BOUNDS
+    }
+}
+
+
+def _get_utm_crs_list_from_bounds(
+        bounds=(-180., -90., 180., 90.),
+        datum_name="WGS 84"
+):
+    out_utm_epsg_list = []
+    utm_pyproj_crs_list = query_utm_crs_info(
+        datum_name=datum_name,
+        area_of_interest=AreaOfInterest(
+            west_lon_degree=bounds[0],
+            south_lat_degree=bounds[1],
+            east_lon_degree=bounds[2],
+            north_lat_degree=bounds[3],
+        )
+    )
+    for c in utm_pyproj_crs_list:
+        out_utm_epsg_list.append(c.code)
+    return out_utm_epsg_list
+
+
+def _get_utm_pyramid_config(
+        crs_epsg,
+        utm_grid="EOX_UTM"
+):
+    utm_grid_params = UTM_GRIDS[utm_grid]
+    if crs_epsg.startswith("327"):
+        grid_bounds = utm_grid_params["south_bounds"]
+    else:
+        grid_bounds = utm_grid_params["north_bounds"]
+    out_utm_config_dict = {}
+    out_utm_config_dict[f"{utm_grid}_{crs_epsg}"] = {
+        'shape': utm_grid_params['shape'],
+        'bounds': grid_bounds,
+        'srs': {"epsg": crs_epsg},
+        'is_global': False,
+    }
+    return out_utm_config_dict
+
+
+def _get_utm_pyramid_configs(grids=UTM_GRIDS, bounds=[-180, -90, 180, 90]):
+    out_dict = {}
+    utm_epsg_list = _get_utm_crs_list_from_bounds(
+        bounds=bounds
+    )
+    for utm_grid in grids.keys():
+        for utm_epsg in utm_epsg_list:
+            out_dict.update(
+                _get_utm_pyramid_config(
+                    crs_epsg=utm_epsg,
+                    utm_grid=utm_grid
+                )
+            )
+    # print(out_dict)
+    return out_dict
+
+
 # supported pyramid types
 PYRAMID_PARAMS = {
     "geodetic": {
@@ -19,11 +125,17 @@
     "mercator": {
         "shape": (1, 1),
         "bounds": (
-            -20037508.3427892, -20037508.3427892, 20037508.3427892,
-            20037508.3427892),
+            -20037508.3427892,
+            -20037508.3427892,
+            20037508.3427892,
+            20037508.3427892
+        ),
         "is_global": True,
         "srs": {
             "epsg": 3857
         }
     }
 }
+
+# Update Pyramid Parameters with UTM EPSG codes and pyramid setup
+PYRAMID_PARAMS.update(_get_utm_pyramid_configs())

tilematrix/_exceptions.py

@@ -0,0 +1,4 @@
+
+
+class OutOfRangeError(ValueError):
+    """Value not in coordinate allowed values"""

tilematrix/_funcs.py

@@ -2,11 +2,14 @@
 
 from itertools import product, chain
 from rasterio.crs import CRS
-from shapely.geometry import Polygon, GeometryCollection, box
+from shapely.geometry import Point, Polygon, GeometryCollection, box
 from shapely.affinity import translate
 
-from ._conf import DELTA, ROUND
-from ._types import Bounds, Shape
+from tilematrix._conf import DELTA, ROUND
+from tilematrix._utm_coefs import *
+from tilematrix._types import Bounds, Shape
+from tilematrix._exceptions import OutOfRangeError
+from tilematrix._utm import stripe_id_from_point
 
 
 def validate_zoom(zoom):
@@ -237,3 +240,97 @@ def _tile_from_xy(tp, x, y, zoom, on_edge_use="rb"):
         raise ValueError(
                 "on_edge_use '%s' results in an invalid tile: %s" % (on_edge_use, e)
             )
+
+
+def mod_angle(value):
+    """Returns angle in radians to be between -pi and pi"""
+    return (value + math.pi) % (2 * math.pi) - math.pi
+
+
+def in_bounds(x, lower, upper, upper_strict=False):
+    return lower <= x < upper
+
+
+def check_valid_zone(zone_number, zone_letter):
+    if not 1 <= zone_number <= 60:
+        raise OutOfRangeError(
+            'zone number out of range (must be between 1 and 60)'
+        )
+
+    if zone_letter:
+        zone_letter = zone_letter.upper()
+
+        if not 'C' <= zone_letter <= 'X' or zone_letter in ['I', 'O']:
+            raise OutOfRangeError(
+                'zone letter out of range (must be between C and X)'
+            )
+
+
+def latlon_to_utm(
+        latitude,
+        longitude
+):
+    """This function converts Latitude and Longitude to UTM coordinate
+        Parameters
+        ----------
+        latitude: float or NumPy array
+            Latitude between 80 deg S and 84 deg N, e.g. (-80.0 to 84.0)
+        longitude: float or NumPy array
+            Longitude between 180 deg W and 180 deg E, e.g. (-180.0 to 180.0).
+        Returns
+        -------
+        easting: float or NumPy array
+            Easting value of UTM coordinates
+        northing: float or NumPy array
+            Northing value of UTM coordinates
+    """
+    if not in_bounds(latitude, -80, 84):
+        raise OutOfRangeError(
+            'latitude out of range (must be between 80 deg S and 84 deg N)'
+        )
+    if not in_bounds(longitude, -180, 180):
+        raise OutOfRangeError(
+            'longitude out of range (must be between 180 deg W and 180 deg E)'
+        )
+
+    lat_rad = math.radians(latitude)
+    lat_sin = math.sin(lat_rad)
+    lat_cos = math.cos(lat_rad)
+
+    lat_tan = lat_sin / lat_cos
+    lat_tan2 = lat_tan * lat_tan
+    lat_tan4 = lat_tan2 * lat_tan2
+
+    zone_number = int(stripe_id_from_point(Point(longitude, latitude))[:2])
+
+    lon_rad = math.radians(longitude)
+    central_lon = (zone_number - 1) * 6 - 180 + 3
+    central_lon_rad = math.radians(central_lon)
+
+    n = R / math.sqrt(1 - E * lat_sin**2)
+    c = E_P2 * lat_cos**2
+
+    a = lat_cos * mod_angle(lon_rad - central_lon_rad)
+    a2 = a * a
+    a3 = a2 * a
+    a4 = a3 * a
+    a5 = a4 * a
+    a6 = a5 * a
+
+    m = R * (M1 * lat_rad -
+             M2 * math.sin(2 * lat_rad) +
+             M3 * math.sin(4 * lat_rad) -
+             M4 * math.sin(6 * lat_rad))
+
+    easting = K0 * n * (a +
+                        a3 / 6 * (1 - lat_tan2 + c) +
+                        a5 / 120 * (5 - 18 * lat_tan2 + lat_tan4 + 72 * c - 58 * E_P2)) + 500000
+
+    northing = K0 * (m + n * lat_tan * (a2 / 2 +
+                                        a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c**2) +
+                                        a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)))
+
+    if latitude <= 0:
+        northing += 10000000
+
+    return easting, northing

tilematrix/_utm.py

@@ -0,0 +1,55 @@
+import math
+
+STRIPE_WIDTH_DEG = 6.
+UTM_TOP = 84.
+UTM_BOTTOM = -80
+UTM_STRIPE_COUNT = 60
+
+
+def utm_stripe_epsg(stripe_id):
+    """
+    Return EPSG code from UTM stripe ID.
+
+    stripe_id: e.g. 33N --> EPSG:32633
+    32: UTM EPSG prefix
+    6 or 7: N or S
+    33: UTM stripe
+    """
+    return "EPSG:32%s%s" % (
+        "6" if _hemisphere(stripe_id) == "N" else "7",
+        _stripe(stripe_id)
+    )
+
+
+def stripe_id_from_point(point):
+    """
+    Return UTM stripe ID from point.
+
+    Parameters
+    ----------
+    point : Point geometry in EPSG:4326
+    """
+    stripe = int(math.ceil((180. + point.x) / STRIPE_WIDTH_DEG))
+    if 0. <= point.y <= UTM_TOP:
+        hemisphere = "N"
+    elif UTM_BOTTOM <= point.y < 0.:
+        hemisphere = "S"
+    else:
+        raise ValueError("point outside UTM bounds: %s", repr(point))
+    return "%s%s" % (str(stripe), hemisphere)
+
+
+def _stripe(stripe_id):
+    s = stripe_id[:2]
+    if int(s) > UTM_STRIPE_COUNT:
+        raise ValueError("invalid UTM stripe ID: %s" % stripe_id)
+    return s
+
+
+def _hemisphere(stripe_id):
+    if not isinstance(stripe_id, str) or len(stripe_id) != 3:
+        raise TypeError("invalid UTM stripe ID: %s" % stripe_id)
+    hemisphere = stripe_id[2]
+    if hemisphere not in ["N", "S"]:
+        raise ValueError("invalid UTM stripe ID: %s" % stripe_id)
+    return hemisphere

tilematrix/_utm_coefs.py

@@ -0,0 +1,27 @@
+import math
+
+K0 = 0.9996
+
+E = 0.00669438
+E2 = E * E
+E3 = E2 * E
+E_P2 = E / (1 - E)
+
+SQRT_E = math.sqrt(1 - E)
+_E = (1 - SQRT_E) / (1 + SQRT_E)
+_E2 = _E * _E
+_E3 = _E2 * _E
+_E4 = _E3 * _E
+_E5 = _E4 * _E
+
+M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
+M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
+M3 = (15 * E2 / 256 + 45 * E3 / 1024)
+M4 = (35 * E3 / 3072)
+
+P2 = (3 / 2 * _E - 27 / 32 * _E3 + 269 / 512 * _E5)
+P3 = (21 / 16 * _E2 - 55 / 32 * _E4)
+P4 = (151 / 96 * _E3 - 417 / 128 * _E5)
+P5 = (1097 / 512 * _E4)
+
+R = 6378137