Deployed Planetcantile 49900 - Unable to get Tiles

[jlaura]

See the following example:

import math

def get_zxy(lat, lon, zoom):
    x = math.floor((lon + 180) / 360 * (1 << zoom))
    rlat = math.radians(lat)
    y = math.floor((1 - math.log(math.tan(rlat)) + 1 / math.cos(rlat)) / math.pi / 2 * (1 << zoom))
    return z, x, y

lat = 20.7242
lon = 78.8926
z, x, y = get_zxy(lat, lon, z)

api = 'https://astrogeology.usgs.gov/titiler'
cog = 'https://astrogeo-ard.s3-us-west-2.amazonaws.com/mars/mro/hirise/uncontrolled/ESP_077499_2010_COLOR/ESP_077499_2010_COLOR.tif'

query_url=f"{api}/cog/tiles/IAU_2015_49900/{z}/{x}/{y}?url={cog}"
print(query_url)

import requests
r = requests.get(query_url, stream=True)
print(r.status_code)

The response is a 404. The COG is projected to "Mars (2015) - Sphere / Ocentric / Equirectangular, clon = 0" according to gdalinfo (don't take my word for it!).

Then take this example of a THEMIS Coprates quad. Note that is this passing 15 to the get_zxy func because we can't take the log of math.radians(-15).

import math

def get_zxy(lat, lon, zoom):
    x = math.floor((lon + 180) / 360 * (1 << zoom))
    rlat = math.radians(lat)
    y = math.floor((1 - math.log(math.tan(rlat)) + 1 / math.cos(rlat)) / math.pi / 2 * (1 << zoom))
    return z, x, y

lat = 15
lon = 22.5
z, x, y = get_zxy(lat, lon, z)

api = 'https://astrogeology.usgs.gov/titiler'
cog = 'https://astrogeo-ard.s3.us-west-2.amazonaws.com/mars/mo/themis/controlled_mosaics/mc18_coprates_nir/mc18_coprates_nir.tif'

query_url=f"{api}/cog/tiles/IAU_2015_49900/{z}/{x}/{y}?url={cog}"
print(query_url)

import requests
r = requests.get(query_url, stream=True)
print(r.status_code)

The result of the second code block is a 200 (nice looking png). In terms of projection, the COGs are identical. This one has me stumped.

This is a version of titiler using main off of this repo, so completely up to date.

[AndrewAnnex]

try using https://github.com/developmentseed/morecantile/blob/main/morecantile/models.py#L687 to get the correct tile index for a given lat/lon/zoom. Likely titiler is responding with a error message about the tile being out of bounds in some way

[jlaura]

Zoom can be anything that the TMS supports. 4 or 9 are good starting points.

[AndrewAnnex]

@jlaura as I suggested your get_zxy function is wrong. Using the tms.tile method for that lon/lat I get a valid tile, if you had looked at your response content you would have seen b'{"detail":"Tile 9/368/247 is outside https://astrogeo-ard.s3-us-west-2.amazonaws.com/mars/mro/hirise/uncontrolled/ESP_077499_2010_COLOR/ESP_077499_2010_COLOR.tif bounds"}'.

Using morecantile's tile function I instead get the following url https://astrogeology.usgs.gov/titiler/cog/tiles/IAU_2015_49900/9/736/197?url=https://astrogeo-ard.s3-us-west-2.amazonaws.com/mars/mro/hirise/uncontrolled/ESP_077499_2010_COLOR/ESP_077499_2010_COLOR.tif which returns a valid (unstyled) tile

[jlaura]

I don't know if the issue is here or if the issue is in leaflet, but these TMS are definitely not working as expected. The x,y,z tile coordinates that you used only appear in the API logs via leaflet if one specifies a 128x128 tile size. TiTiler supports queries on 256+ tile sizes. This mismatch does not make sense.

Below is the example I have been testing in a notebook. If I change the tile size to be 256 and see that TiTiler is responding with 404 (and the x,y,z tile coordinates are wrong). I am also not convinced that setting tilesize to 128 is actually correct. Passing rescale=358,705 to the API (which are the 2%,98% clips on the unscaled data) still returns black tiles.

api = 'https://astrogeology.usgs.gov/titiler'
url = 'https://stac.astrogeology.usgs.gov/api/collections/mro_hirise_uncontrolled_observations/items/ESP_077499_2010_COLOR'
resp = requests.get(f'{api}/stac/bounds?url={url}')
bounds = resp.json()['bounds']
center = (bounds[1] + (bounds[3] - bounds[1])/2, bounds[0] + (bounds[2] - bounds[0]) / 2)
print(center)

from ipyleaflet import Map, TileLayer
m = Map(
    zoom=7,
    center=center,
)

cog = 'https://astrogeo-ard.s3-us-west-2.amazonaws.com/mars/mro/hirise/uncontrolled/ESP_077499_2010_COLOR/ESP_077499_2010_COLOR.tif'
query_url=f"{api}/cog/tiles/IAU_2015_49915/{{z}}/{{x}}/{{y}}?url={cog}&return_mask=false&rescale=358,705"

layer = TileLayer(
    url=query_url,
    opacity=1,
    tile_size=128,
    show_Loading=True
)
m.add_layer(layer)
m

[AndrewAnnex]

@jlaura it's a leaflet issue. Try setting the map to use EPSG4326 CRS explicitly like:

m = Map(
    center=center,
    zoom=7,
    crs='EPSG4326',
)

leaflet by default queries mercator tile ZXY coordinates, and the tiles for mercator TMSs do not match global equirectangular tile indexes used for 49900. Once we get merctor codes and TMSs titiler will be able to work with leaflet defaults more easily.

Also past the response contents next time for the 404 errors, titiler will tell you if the tile bounds are invalid and you can debug that by trying to get the xy bounds for the tile queried by leaflet