Allow Non-Earth Coordinates System

[jlaura]

While testing TiTiler with an S3 hosted COG and I am running into some (not unexpected) issues with the fact that the data are of another body (Europa in this case. Also tested Mars).

I am seeing the following error when attempting to hit the /cog/bounds endpoint.

File "rasterio/_err.pyx", line 215, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from '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]]]]' to 'EPSG:4326'

Speaking with the devs over on that project it sounds like this could be an issues with how the bounds are being computed. The morecantile project pulls a transformation object from pyproj. I wonder if a similar solution might be viable here so that the coordinate transformation can be generated?

Also, here is the direct link to the data if anyone is interested! https://asc-jupiter.s3-us-west-2.amazonaws.com/europa/galileo_voyager/controlled_mosaics/11ESCOLORS01-02_GalileoSSI_Equi-cog.tif

[vincentsarago]

👋 @jlaura

The problem is that by default rio-tiler tries to get the bounding box of the data in WGS84 (earth)

rio-tiler/rio_tiler/io/cogeo.py

Lines 103 to 105 in 809ee70

	self.bounds = transform_bounds(
	self.dataset.crs, WGS84_CRS, *self.dataset.bounds, densify_pts=21
	)

and also get the min/max zoom in selected TMS (default is WebMercator)

rio-tiler/rio_tiler/io/cogeo.py

Lines 131 to 151 in 809ee70

	if self.dataset.crs != self.tms.crs:
	dst_affine, w, h = calculate_default_transform(
	self.dataset.crs,
	self.tms.crs,
	self.dataset.width,
	self.dataset.height,
	*self.dataset.bounds,
	)
	else:
	dst_affine = list(self.dataset.transform)
	w = self.dataset.width
	h = self.dataset.height
	resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
	maxzoom = self.tms.zoom_for_res(resolution)
	overview_level = get_maximum_overview_level(w, h, minsize=tilesize)
	ovr_resolution = resolution * (2 ** overview_level)
	minzoom = self.tms.zoom_for_res(ovr_resolution)
	return minzoom, maxzoom

both will fail for non-earth data (there is no possible coordinates transformation between earth and non-earth projection).

I'm currently working on this but it will require a lot of change, which will then be added to #426 for our next major release.

As mentioned over TiTiler, in order to fully utilise titiler/rio-tiler you'll need a TileMatrixSet for Europa (you won't be able to use WebMercator TMS to display your data). I was wondering if you had already created one. I don't know the specific of the coordinate system uses for your data or if there are standard TMS for non-earth body 🤷‍♂️ (having a TMS could help me during testing).

[jlaura]

@vincentsarago We must have a TMS defined for at least Mars because another provider serves up a WMTS (https://api.nasa.gov/mars-wmts/catalog/). The data are all in the projection that I am used to working in IAU2000:49900. For the moon, I would anticipate using IAU2000:30100.

Here are all of the IAU projection codes in wkt format.

I believe that the IAU projections should be getting into the proj library at some point, but I do not see them in the 8.1.1 release. I'll check with our liaison to that project. I also started to dig into the OGC TMS specification to see what all is required and will ping @thareUSGS to see if he has any additional information.

[jlaura]

Minor addendum - the codes are in proj (looking in 8.1.1), but under the ESRI namespace (not, IAU). For example, ESRI:104905 "GCS_Mars_2000" is, I believe, IAU2000:49900. Likewise, Europa is ESRI:104915 "GCS_Europa_2000".

The 2015 IAU codes look like they will be in Proj 8.2 🚀 (OSGeo/PROJ#2876)

Does having these help at all?

[vincentsarago]

@jlaura that's excellent info

we can then construct a custom TMS:

from pyproj import CRS

europa_crs = CRS.from_authority("ESRI", 104915)
print(europa_crs)

>>> <Geographic 2D CRS: ESRI:104915>
Name: GCS_Europa_2000
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: D_Europa_2000
- Ellipsoid: Europa_2000_IAU_IAG
- Prime Meridian: Reference_Meridian

from morecantile import TileMatrixSet

tms = TileMatrixSet.custom(
    crs=europa_crs,
    extent=europa_crs.area_of_use.bounds,
    matrix_scale=[2, 1],
)

# Checking if rasterio can do transformation from TMS crs and Image crs (the one you gave)
import rasterio
from rasterio.warp import transform_bounds

with rasterio.open("cog_nonearth.tif") as src:
    transform_bounds(src.crs, tms.rasterio_crs, *src.bounds)

>>> (129.36834223297478, 13.985559117409744, 138.90253908503576, 23.13673177454536)

Short story, I think it's possible to make rio-tiler/morecantile/titiler work with the data :-)

[jlaura]

@vincentsarago Awesome! Looks like I need to use the docs here - using a custom tms with the example that you used above and I should be good to go.

I'll give that a look this evening (UTC-7) and report back on the discussion over on TiTiler.

Much appreciate the information and rapid assist!

[vincentsarago]

FYI, for now rio-tiler will still error on non-earth images. I hope to finish my work over #429 later today.

[vincentsarago]

FYI, creating a tiler using #429 and the custom TMS shared earlier should work fine

I've added a notebook in #429 https://github.com/cogeotiff/rio-tiler/blob/00d7f283c47d1ca9d5df5b3ca68e5e350c9a6576/docs/examples/Using-nonEarth-dataset.ipynb