Use of Custom TMS (Non-Earth)

[jlaura]

I am seeing the following error attempting to use a custom TMS (for non-earth data). I wonder if anything jumps out at anyone? I have a working example just using riotiler (many thanks @vincentsarago for seeing that issue and posting the example!). The next step, is to get a titiler example working. Then I can deploy body specific tilers to AWS (that's the plan anyway!).

Here is the error I am seeing:

rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from 'EPSG:3857' to 'PROJCRS["Equirectangular EUROPA",BASEGEOGCRS["GCS_EUROPA",DATUM["D_EUROPA",ELLIPSOID["EUROPA_localRadius",1560800,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Equidistant Cylindrical (Spherical)",ID["EPSG",1029]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",180,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

So, no transformation operations from Europa to 4326. Since both are just lat/lon, I was anticipating that this hack would work? Perhaps I am messing up the configuration of the custom WMS?

Here is a minimal reproducible app:

dependencies.py

from enum import Enum
from morecantile import tms, TileMatrixSet
from pyproj import CRS
from fastapi import Query

# 1. Create Custom TMS
europa_crs = CRS.from_authority("ESRI", 104915)
europa_tms = TileMatrixSet.custom(
    crs=europa_crs, extent=europa_crs.area_of_use.bounds, matrix_scale=[2, 1],
    name='Europa'
)

# 2. Register TMS
tms = tms.register([europa_tms,])
print(tms)
# 3. Create ENUM with available TMS
TileMatrixSetNames = Enum(  # type: ignore
    "TileMatrixSetNames", [(a, a) for a in sorted(tms.list())]
)

print('SetNames', TileMatrixSetNames)

# 4. Create Custom TMS dependency
def TMSParams(
    TileMatrixSetId: TileMatrixSetNames = Query(
        TileMatrixSetNames.WebMercatorQuad,  # type: ignore
        description="TileMatrixSet Name (default: 'WebMercatorQuad')",
    )
) -> TileMatrixSet:
    """TileMatrixSet Dependency."""
    return tms.get(TileMatrixSetId.name)

routes.py

from titiler.core.factory import TilerFactory, TMSFactory
from .dependencies import TileMatrixSetNames, TMSParams

tms = TMSFactory(supported_tms=TileMatrixSetNames, tms_dependency=TMSParams)
cog = TilerFactory(tms_dependency=TMSParams)

main.py

from titiler.core.errors import DEFAULT_STATUS_CODES, add_exception_handlers

from fastapi import FastAPI

from .routes import cog, tms

app = FastAPI(title="My simple app with custom TMS")

app.include_router(cog.router, tags=["Cloud Optimized GeoTIFF"])
app.include_router(tms.router, tags=["TileMatrixSets"])
add_exception_handlers(app, DEFAULT_STATUS_CODES)

web map (Jupyter)

bounds = (90, 7.8, 90.5, 8.3)

m = Map(
    center=(
        (bounds[1] + bounds[3]) / 2,
        (bounds[0] + bounds[2]) / 2
    ),
    zoom=4,
    crs=projections.EPSG4326,  # HACK: the europa TMS is in degree and covers -180, -90, 180, 90 like the WGS84
)

layer = TileLayer(
    url="http://127.0.0.1:8000/cog/tiles/{z}/{x}/{y}?url=https://asc-jupiter.s3.us-west-2.amazonaws.com/europa/galileo_voyager/controlled_images/equi/s0383713765/s0383713765_photr.tif",
    min_zoom=4,
    max_zoom=6,
    opacity=1,
)
m.add_layer(layer)
m

[vincentsarago]

@jlaura 👁️ your code it seems that the default TMS used is WebMercator, can you change to Europa

# 4. Create Custom TMS dependency
def TMSParams(
    TileMatrixSetId: TileMatrixSetNames = Query(
        TileMatrixSetNames. Europa,  # type: ignore
        description="TileMatrixSet Name (default: 'Europa')",
    )
) -> TileMatrixSet:
    """TileMatrixSet Dependency."""
    return tms.get(TileMatrixSetId.name)

[vincentsarago]

Oh I found why it was failing on my side. You need to make sure that rasterio/gdal use PROJ>=8.0! (FYI, homebrew compile gdal with proj7)

[jlaura]

Wonderful and thank you for trying to reproduce!

I have proj 8.2.1, libgdal 3.4.1, Rasterio 1.2.10. I am working on learning the code base. It seems to me that ti-tiler is using riotiler to do the read through the Tiler Factory. Could it be the case that the injection of the custom TMS is not being caught by the reader?

I think the error means that this is not the case, but I'm not 100% sure where the issue could be otherwise.

[vincentsarago]

Could it be the case that the injection of the custom TMS is not being caught by the reader?

I don't think, it's well tested!

[jlaura]

I have been able to get Mars and Europa support working! Whatever the issue was, the new releases have fixed it.By modifying the custom.py file I can set:

mars_crs = CRS.from_authority("ESRI", 104971)
ESRI104971 = morecantile.TileMatrixSet.custom(
    crs=mars_crs, 
    identifier="Mars",
    extent=mars_crs.area_of_use.bounds, 
    matrix_scale=[2, 2],
)

# CUSTOM TMS for Europa
europa_crs = CRS.from_authority("ESRI", 104915)
ESRI104915 = morecantile.TileMatrixSet.custom(
    crs=europa_crs, 
    identifier="Europa",
    extent=europa_crs.area_of_use.bounds, 
    matrix_scale=[2, 2],
)

tms = tms.register([ESRI104915, ESRI104971])

and then

def TMSParams(
    TileMatrixSetId: TileMatrixSetName = Query(
        TileMatrixSetName.Europa,  # type: ignore
        description="TileMatrixSet Name (default: 'Europa')",
    )
) -> TileMatrixSet:
    """TileMatrixSet Dependency."""
    return tms.get(TileMatrixSetId.name)

I then can use the /cog/tiles endpoint to get Europa tiles.

Two questions remain:

1. Is it possible to deploy with multiple different TMSParams being passed to the TileFactory? This way a single service could support multiple different bodies? Or, should I deploy multiple lambdas, one for each body?
2. The data (in S3) are 32-bit. I am having to pass the rescale=x,y parameter along with the query string to the URL. Is rio-tiler pulling any kind of no data value information or statistics off of the image before tiling? I am trying to determine if I need to do something different with my data or if I should be needing to pass the rescale bounds.

Thanks!

titiler.__version__ == 0.4.2 😉

[vincentsarago]

Is it possible to deploy with multiple different TMSParams being passed to the TileFactory? This way a single service could support multiple different bodies? Or, should I deploy multiple lambdas, one for each body?

(I'm not sure to fully understand what you want here, so I'll take a guess) you can create endpoints using different TMSParams as

TMSEuropa = Enum(  # type: ignore
    "TMSEuropa", [("Europa", "Europa")]
)

def TMSParamsEurope(
    TileMatrixSetId: TMSEuropa = Query(
        TMSEuropa.Europa,  # type: ignore
        description="TileMatrixSet Name (default: 'Europa')",
    )
) -> TileMatrixSet:
    """TileMatrixSet Dependency."""
    return tms.get(TileMatrixSetId.name)


europa = TilerFactory(tms_dependency= TMSParamsEurope, router_prefix="/europa")

TMSMars = Enum(  # type: ignore
    "TMSMars", [("Mars", "Mars")]
)

def TMSParamsMars(
    TileMatrixSetId: TMSMars= Query(
        TMSMars.Mars,  # type: ignore
        description="TileMatrixSet Name (default: 'Mars')",
    )
) -> TileMatrixSet:
    """TileMatrixSet Dependency."""
    return tms.get(TileMatrixSetId.name)


mars = TilerFactory(tms_dependency=TMSParamsMars, router_prefix="/mars")


app = FastAPI()
app.include_router(europa.router, prefix="/europa")
app.include_router(mars.router, prefix="/mars")

☝️ then you'll have 2 set of endpoints /europa/tiles/... /mars/tiles/...

The data (in S3) are 32-bit. I am having to pass the rescale=x,y parameter along with the query string to the URL. Is rio-tiler pulling any kind of no data value information or statistics off of the image before tiling? I am trying to determine if I need to do something different with my data or if I should be needing to pass the rescale bounds.

This is a tricky part with rio-tiler, we don't do any rescaling if the user don't ask to. If the data require rescaling we usually ask the user to first use the /statistics endpoint to get the min/max and then forward the rescaling parameter to the tile endpoint (or tilejson endpoint)

FYI: I'm planning to use geotiff metadata to allow automatic rescaling cogeotiff/rio-tiler#488 but I don't have timeline for this