Issue with resample_geometry + regridded Image/RGB with non-lat/lon bounds

[TheoMathurin]

ALL software version info

GeoViews 1.9.5, HoloViews 1.15.1 on Linux

The code provided below is executed in a notebook (jupyter_client 7.0.6, jupyter_core 4.8.1, jupyter_server 1.11.1, notebook 6.4.4) but I'm getting the same bug in a server (Panel 0.14.0).

This is not a regression as I've experienced this issue for more than a year with various versions.

Description of expected behavior and the observed behavior

When working in geostationary projection (or orthographic and presumably others) with Image or RGB elements, it is preferable to provide bounds calculated in the projection coordinates instead of lat/lon. This is because the corners of the element may not fit within the globe and therefore may not correspond to a valid (lat, lon) pair:

When the image is large, regridding using hd.regrid allows to drastically reduce rendering time. Likewise, using resample_geometry on Feature objects (coastline typically) also helps.

However you can't use both at the same time. If you do, the plot does render, but the moment you pan or zoom, one of the element disappears (geometry or image, depending on their order in the overlay). The same thing happens if you use hv.Image/hv.RGB instead of their gv counterparts, which I do (it's a bit faster).

So in this case you have to choose between image regridding and feature resampling. Apparently the bug does not occur if you provide bounds in lat/lon.

Complete, minimal, self-contained example code that reproduces the issue

import numpy as np
import geoviews as gv
import holoviews.operation.datashader as hd
import cartopy.crs as ccrs

gv.extension('bokeh')

coastline = gv.feature.coastline.geoms('10m').opts(color='black')
resampled_coastline = gv.operation.resample_geometry(coastline)

projection = ccrs.Orthographic()
opts = dict(projection=projection, width=400, height=400)

ls = np.linspace(0, 10, 2000)
xx, yy = np.meshgrid(ls, ls)

# In separate cells:
latlon_bounds = (-20, -10, 20, 10)
latlon_img = gv.Image(np.sin(2*xx)*np.cos(yy), bounds=latlon_bounds)
regrid_latlon_img = hd.regrid(latlon_img)
(regrid_latlon_img*resampled_coastline).opts(**opts) # Works fine

# and
geo_bounds = (-2000000, -1000000, 2000000, 1000000)
geo_img = gv.Image(np.sin(2*xx)*np.cos(yy), bounds=geo_bounds, crs=projection)
# or geo_img = hv.Image(np.sin(2*xx)*np.cos(yy), bounds=geo_bounds)
regrid_geo_img = hd.regrid(geo_img)
(regrid_geo_img*resampled_coastline).opts(**opts) # Coastline will disappear as soon as you pan or zoom

Stack traceback and/or browser JavaScript console output

Nothing is logged in standard output or in the browser console

Screenshots or screencasts of the bug in action

coastline_bug.mp4

[TheoMathurin]

Looking at the _process method in resample_geometry, it turns out the view port is expressed in geostationary coordinates (from self.p.x_range and self.p.y_range) when regridding is enabled, instead of lat/lon when it is not.

As a result the assumed view port is wrong with respect to the vertices (lat/lon coordinates) and the new_geoms list is not properly populated.

When the view port is artificially set to something that makes sense in lat/lon, the geometries persist. Therefore it seems that the key would be to make sure self.p.x_range/self.p.y_range return lat/lon coordinates.

[TheoMathurin]

Conversely, when the order is reversed in the overlay (coastline*img and the image element disappears), then in the regrid function the x_range/y_range are wrong (the bounds are the same i.e. the ranges are zero).

[TheoMathurin]

I have tried various combinations of the apply_extents parameter, to no avail. I thought it was relevant, from its description in the docs. Specifically I thought that by disabling extents override on the geometry element, the geostationary bounds of the regridded image would not be applied. But apparently I was misunderstood.

[TheoMathurin]

I've just found a workaround, using a DynamicMap driven by a RangeXY stream which returns gv.operation.resample_geometry with dynamic=False. In the function, you can convert the projected coordinates back to lat/lon and feed that to resample_geometry.

xy_range = RangeXY()    

def resample_coastline(x_range, y_range):
    geometry = gv.feature.coastline.geoms('50m')
    # Convert the projected coordinates provided by xy_range back to latitude/longitude
    x_range, y_range = convert_coords(x_range, y_range)
    return gv.operation.resample_geometry(geometry, x_range=x_range, y_range=y_range, dynamic=False)

coastline = hv.DynamicMap(resample_coastline, streams=[xy_range])

[MoritzImendoerffer]

Hi @TheoMathurin, out of curiosity, how does your function convert_coords look like?

[TheoMathurin]

Hi @MoritzImendoerffer

I use the Transformer object in pyproj, something like:

from pyproj import CRS, Transformer

def convert_coords(x, y, geo_satellite):
    proj4_string = ... # build string based on geo_satellite information
    crs_geo = CRS.from_proj4(proj4_string)
    transformer = Transformer.from_crs(4326, crs_geo)
    return transformer.transform(x, y)

When it returns invalid results because latitude/longitude are undefined (points in space), I manually set them to default values.

[TheoMathurin]

I hoped that #722 and/or #738 would fix this issue but unfortunately it's still there on 1.13.0:

import numpy as np
import geoviews as gv
import holoviews.operation.datashader as hd
import cartopy.crs as ccrs

gv.extension('bokeh')

coastline = gv.feature.coastline.geoms('10m').opts(color='black')
resampled_coastline = gv.operation.resample_geometry(coastline)

projection = ccrs.Orthographic()
opts = dict(projection=projection, width=400, height=400)

ls = np.linspace(0, 10, 2000)
xx, yy = np.meshgrid(ls, ls)

# and
geo_bounds = (-2000000, -1000000, 2000000, 1000000)
geo_img = gv.Image(np.sin(2*xx)*np.cos(yy), bounds=geo_bounds, crs=projection)
# or geo_img = hv.Image(np.sin(2*xx)*np.cos(yy), bounds=geo_bounds)
regrid_geo_img = hd.regrid(geo_img)
(regrid_geo_img*resampled_coastline).opts(**opts) # Coastline will disappear as soon as you pan or zoom