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.