Feature request: Support non-standard CRS in mapping functions

[s-m-t-c]

I would like to plot EPSG:3031 data using something like data.odc.add_to(map).

I have been experimenting with solutions using something like below which currently does not work.

Others have encountered similar issues jupyter-widgets/ipyleaflet#550 though without the same projection requirements.

For an example input we can use a sea ice concentration image

from ipyleaflet import Map, basemaps, ImageOverlay
import folium
import xarray as xr
import rioxarray as rio
from numpy.random import uniform
from odc.geo.data import country_geom
from odc.geo.xr import rasterize
from osgeo import gdal,ogr,osr

def GetExtent(ds):
    """ Return list of corner coordinates from a gdal Dataset """
    xmin, xpixel, _, ymax, _, ypixel = ds.GetGeoTransform()
    width, height = ds.RasterXSize, ds.RasterYSize
    xmax = xmin + width * xpixel
    ymin = ymax + height * ypixel

    return (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)

path = 'S_201905_concentration_v3.0.tif'

src = gdal.Open(path)

bounds = GetExtent(src)

spsLayout=Layout(width='800px', height='800px')

m = Map(basemap=basemaps.NASAGIBS.BlueMarble3031, center=(-90, 0), zoom=1, crs=projections.EPSG3031, layout=spsLayout)

image = ImageOverlay(
    url=path,
    bounds=(bounds[1], bounds[3])
)

m.add_layer(image)

m

[Kirill888]

1. Detect Map projection inside add_to from passed in map, and use that instead of assuming epsg:3857
2. Add optional crs argument to .map_bounds(crs) and figure out what leaflet.js expects to see there. Looking through javascript code I think it expects location of two corner points of the image projected to latlon and NOT actual bounding box of the image in lat/lon. The two corners I believe are (0,0) and (W,H).

So, if you have gbox of the image you can do this:

(x1, y1), _, (x2, y2) = gbox.extent.exterior.to_crs("epsg:4326").points[:3]
bounds = [[y1, x1], [y2, x2]]

[Kirill888]

@s-m-t-c this below works:

import numpy as np
from ipyleaflet import Map, basemaps, ImageOverlay, projections
import rioxarray as rio
import odc.geo.xr

# wget https://masie_web.apps.nsidc.org/pub//DATASETS/NOAA/G02135/south/monthly/geotiff/08_Aug/S_201908_concentration_v3.0.tif

path = 'S_201908_concentration_v3.0.tif'
yy = rio.open_rasterio(path).isel(band=0)

# 1. bring to CRS of the map
xx = yy.odc.reproject("epsg:3031")
display(xx.odc.geobox)

# 2. Convert to RGB and compress into data url
vmin, vmax = np.nanpercentile(xx.data, [2, 98])
rgba = xx.odc.colorize(vmin=vmin, vmax=vmax)
rgba.plot.imshow(aspect=1, size=4);
data_url = rgba.odc.compress(as_data_url=True)

# 3. TL and BR corner coordinates projected to LatLon makes bounds
(x1, y1), _, (x2, y2) = rgba.odc.geobox.extent.exterior.to_crs("epsg:4326").points[:3]
bounds = [[y1, x1], [y2, x2]]

# 4. Make map and display
m = Map(basemap=basemaps.NASAGIBS.BlueMarble3031, center=(-90, 0), zoom=1, crs=projections.EPSG3031)
m.add_layer(ImageOverlay(url=data_url, bounds=bounds, opacity=0.7))
display(m)

so coordinates of top-left and bottom-right corners in lat/lon order is what bounds= actually expects, it just happens to match bounding box in axis aligned images in epsg:3857 and epsg:4326.

[Kirill888]

This was solved by #54, should work in odc-geo==0.2.1