pythongistiffholoviewsdatashader

Geotiff overlay position is slightly off on Holoviews/Bokeh tilemap


I have a Geotiff that I display on a tile map, but it's slightly off to the south. For example, on this screenshot the edge of the image should be where the country border is, but it's a bit to the south:

screenshot of the problem

Here's the relevant part of the code:

tiff_rio_500 = rioxarray.open_rasterio('/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif')
dataarray_500 = tiff_rio_500[0]
dataarray_500_meters = dataarray_500.copy()
dataarray_500_meters['x'], dataarray_500_meters['y'] = ds.utils.lnglat_to_meters(dataarray_500.x, dataarray_500.y)
hv_dataset_500_meters = hv.Dataset(dataarray_500_meters, name='nightlights', vdims='cumulative_cost')
hv_tiles_osm_bokeh = hv.element.tiles.OSM().opts(width=1000, height=800)
hv_image_500_meters_bokeh = hv.Image(hv_dataset_500_meters, kdims=['x', 'y'], vdims=['cumulative_cost'], rtol=1).opts(cmap='inferno_r')
hv_combined_osm_500_meters_bokeh = hv_tiles_osm_bokeh * hv_image_500_meters_bokeh
hv_combined_osm_500_meters_bokeh

You can see the live notebook on google colab.

Now this is not the usual "everything is way off" problem that occurs when one doesn't convert the map to Web Mercator. It is almost perfect, it just isn't.

The Geotiff is an Earth Engine export. This is how it looked originally in Earth Engine (see live code):
screenshot of everything working fine in Earth Engine
As you can see, the image follows the borders everywhere.

At first, I suspected that maybe the export went wrong, or the google map tileset is somewhat different, but no, if I open the same exported Tiff in the QGis application on my windows laptop and view it on the same OSM tilemap as I do in the colab notebook, it looks fine:
screenshot of everything being fine in QGis
Okay, the image does not follow the borders perfectly, but I know why and that's unrelated (I oversimplified the country border geometry). The point is, that it is projected to the correct location. So based on that, the tiff contains the correct information, it can be displayed at the same location as the borders are in the OSM tilemap, but still in my Holoviews-Datashader-Bokeh project it is slightly off.

Any idea why this happens?


Solution

  • I've got the answer on the Holoviz Discourse from one of the developers. Seeing how the recommended function is practically undocumented, I copy it here in case somebody looks for an easy way to load a geotiff and add to a tilemap in Holoviews/Geoviews:

    https://discourse.holoviz.org/t/geotiff-overlay-position-is-slightly-off-on-holoviews-bokeh-tilemap/2071

    philippjfr
    I wouldn’t expect manually transforming the coordinates to work particularly well. While it’s a much heavier weight dependency for accurate coordinate transforms I’d recommend using GeoViews.

    img = gv.util.load_tiff( '/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif' ) gv.tile_sources.OSM() * img.opts(cmap='inferno_r')

    Edit: Now it is possible one doesn't want to use Geoviews as it has a pretty heavy dependency chain that requires a lot of patience and luck to set it up right. Fortunately rioxarray (through rasterio) has a tool to reproject, just append ".rio.reproject('EPSG:3857')" to the first line and then you don't have to use the lnglat_to_meters which is not intended for this purpose.

    So the corrected code becomes:

    tiff_rio_500 = rioxarray.open_rasterio('/content/mw/mw_dist_to_light_at_all_from_light_mask_mw_cut_s3_500.tif').rio.reproject('EPSG:3857')
    hv_dataset_500_meters = hv.Dataset(tiff_rio_500[0], name='nightlights', vdims='cumulative_cost')
    hv_tiles_osm_bokeh = hv.element.tiles.OSM().opts(width=1000, height=800)
    hv_image_500_meters_bokeh = hv.Image(hv_dataset_500_meters, kdims=['x', 'y'], vdims=['cumulative_cost'], rtol=1).opts(cmap='inferno_r')
    hv_combined_osm_500_meters_bokeh = hv_tiles_osm_bokeh * hv_image_500_meters_bokeh
    hv_combined_osm_500_meters_bokeh
    

    Now compared to the Geoviews solution (that supposedly handles everything automatically), this solution has a downside that if you use a Hover Tooltip to display the values and coordinates under the mouse cursor, the coordinates are showing up in the newly projected web mercator system in millions of meters instead of the expected degrees. The solution for that is outside the scope of this answer, but I'm just finishing a detailed step by step guide that contains a solution for that too, and I will link that here as soon as it is published. Of course if you don't use Hover Tooltip, the code above will be perfect for you without any more tinkering.