Enable tiling for TMS `WebMercatorQuad` (and other TMS's accepted by the `morecantile` library)

[julietcohen]

When setting the config to tile for the TMS WebMercatorQuad, which differs from the default TMS WGS1984Quad, the x and y identifiers for each tile cannot be identified (only the z component can be identified) when we save the tiles.

Before we save the tiles, when we extract the EPSG code pulled from the default TMS WGS1984Quad here, it looks like: urn:ogc:def:crs:EPSG::4326. When we extract the EPSG code from TMS WebMercatorQuad, it is in the same format: urn:ogc:def:crs:EPSG::3857. This means identifying the CRS of the TMS is not likely the issue.

When we save tiles, we use the method tile_from_str() which pulls the tile x, y, and z identifiers.

The following output identifies 2 options for the source of the error:

1. the staged file is not correctly re-projected into the TMS
   • I added logging to check the CRS after re-projecting, and it is reported as urn:ogc:def:crs:EPSG::3857 which is what we want
   • the point identified as (-18034999.95487016, 10070194.973863265) looks like the units are meters, while the TMS bounds [-180.0, -85.0511287798066, 180.0, 85.0511287798066] are in degrees
2. the tile is re-projected correctly but we cannot identify the x and y

/home/jcohen/anaconda3/envs/arcade_layer/lib/python3.9/site-packages/morecantile/models.py:468: PointOutsideTMSBounds: Point (-18034999.95487016, 10070194.973863265) is outside TMS bounds [-180.0, -85.0511287798066, 180.0, 85.0511287798066].
  warnings.warn(
Traceback (most recent call last):
  File "/home/jcohen/IWP_validation_app/workflow.py", line 177, in <module>
    stager.stage(iwp_paths[0])
  File "/home/jcohen/viz-staging/pdgstaging/TileStager.py", line 139, in stage
    self.save_tiles(gdf)
  File "/home/jcohen/viz-staging/pdgstaging/TileStager.py", line 566, in save_tiles
    self.save_new_tile(data = data,
  File "/home/jcohen/viz-staging/pdgstaging/TileStager.py", line 628, in save_new_tile
    logger.info(f"input tiles_str into tile_from_str is {tiles_str}.")
  File "/home/jcohen/viz-staging/pdgstaging/TileStager.py", line 628, in <listcomp>
    logger.info(f"input tiles_str into tile_from_str is {tiles_str}.")
  File "/home/jcohen/viz-staging/pdgstaging/TilePathManager.py", line 279, in tile_from_str
    x, y, z = [int(i) for i in regex.findall(tile_str)]
ValueError: not enough values to unpack (expected 3, got 1)

config

``` { "deduplicate_clip_to_footprint": True, "dir_input": "/home/pdg/data/ice-wedge-polygon-data/version_2023-01-31/iwp_detections/iwp_files/high/alaska/", "ext_input": ".shp", "ext_footprints": ".shp", "dir_footprints": "/home/pdg/data/ice-wedge-polygon-data/version_2023-01-31/footprints/footprint_files_with_date_20230119/high/alaska/", "dir_staged": "staged/", "dir_geotiff": "geotiff/", "dir_web_tiles": "web_tiles/", "filename_staging_summary": "staging_summary.csv", "filename_rasterization_events": "raster_events.csv", "filename_rasters_summary": "raster_summary.csv", "filename_config": "config", "simplify_tolerance": 0.1, "tms_id": "WebMercatorQuad", "input_crs": None, "z_range": [ 0, 15 ], "geometricError": 57, "z_coord": 0, "statistics": [ { "name": "iwp_coverage", "weight_by": "area", "property": "area_per_pixel_area", "aggregation_method": "sum", "resampling_method": "average", "val_range": [ 0, 1 ], "palette": [ "#f8ff1f1A", # 10% alpha yellow "#f8ff1f" # solid yellow ], "nodata_val": 0, "nodata_color": "#ffffff00" }, ], "deduplicate_at": [ "raster" ], "deduplicate_keep_rules": [ [ "Date", "larger" ] ], "deduplicate_method": "footprints" } ```

Other TMS's that we should be able to process can be found here.

[julietcohen]

When we save any tile, regardless if it will be a new tile or appended to an existing tile, we extract the tile_path here using path_from_tile(tile, base_dir='staged').

Then, if the tile is a new tile, we use that tile_path as an input to save_new_tile() which uses it like so: data.to_file(tile_path, mode=mode), and mode is set to 'w'. The input argument tile_path looks like: staged/WebMercatorQuad/15/nan/nan.gpkg so the issue must occur before we execute save_new_tile().

Because of those nan values, when we move on to define tiles_str with tiles_str = data[self.props['tile']].copy(), tile_str looks like Tile(x=nan, y=nan, z=15).

[julietcohen]

Notes:

• seems like the re-projecting to the CRS of the TMS is working, because the logged CRS after re-projecting is urn:ogc:def:crs:EPSG::3857 (as noted in the first comment) and Grid.py includes a check that the CRS if the TMS is the same as the geodataframe's CRS: __check_crs_match(self, gdf), so the issue is more likely in the extraction of the x and y from the TMS
• Within Grid.py, we execute indices_from_xy(self, x, y) which returns "The row indices and column indices, respectively, that each of the given point coordinates falls within."

When I run the workflow with the default TMS, WGS1984, logging for this step looks like:

INFO:logger:before removing points outside grid:
row_ind is [ 1  1  1 ... 15 15 15]
col_ind is [23 23 23 ... 56 56 56]
INFO:logger:after removing points outside grid:
row_ind is [4241, 4241, 4241, 4241, 4223, 4223, 4223,...
col_ind is [3241, 3241, 3241, 3241, 3218, 3219, 3218,...

When I run the workflow with the desired TMS, WorldMercatorQuad, logging for this step looks like:

INFO:logger:before removing points outside grid:
row_ind is [1 1 1 ... 1 1 1]
col_ind is [-1 -1 -1 ... -1 -1 -1]
INFO:logger:after removing points outside grid:
row_ind is [None, None, None, None, None, None, None,...
col_ind is [None, None, None, None, None, None, None,...

Seems the None values result from this operation. I am not sure if the repeated 1 and -1 values are valid. In contrast, the default TMS shows a large diversity in values for these lists.