Image and QuadMesh not projected correctly latitude-wise

[maximlt]

Image/QuadMesh elements don't seem to be correctly projected latitude-wise, at least in the example below that creates a uniform lat/lon grid whose latitude centers are (-40, -30, ..., 40). The lower bound of the most northern tiles/boxes should be 35 given the grid definition.

import numpy as np, xarray as xr, geoviews as gv, cartopy.crs as ccrs, matplotlib.pyplot as plt
gv.extension('bokeh')

x = np.linspace(0, 10, 5, endpoint=False)
y = np.linspace(-40, 50, 9, endpoint=False)
da = xr.DataArray(np.random.random((5, 9)), coords={'x': x, 'y': y})

gv.Image(da.sel(x=slice(-1, 3), y=slice(30, 40)), kdims=['x', 'y'])

I zoomed on the image below to show that the limit between the boxes is not 35.

This example shows how cartopy displays the same data, projected to Mercator (as done by geoviews under the hood for bokeh). The bound is correctly at 35.

fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.GOOGLE_MERCATOR)

p = da.plot(
    x='x',
    y='y',
    transform=ccrs.PlateCarree(),
    add_colorbar=False,
    ax=ax,
)

p.axes.set_global()

p.axes.coastlines()
p.axes.gridlines(draw_labels=True)
p.axes.set_extent([-1, 3, 32, 38])

[ahuang11]

When projection is set to ccrs.PlateCarree() it's as expected:

import numpy as np, xarray as xr, geoviews as gv, cartopy.crs as ccrs, matplotlib.pyplot as plt
gv.extension('bokeh')

x = np.linspace(0, 10, 5, endpoint=False)
y = np.linspace(-40, 50, 9, endpoint=False)
da = xr.DataArray(np.random.random((5, 9)), coords={'x': x, 'y': y})

gv.Image(da.sel(x=slice(-1, 3), y=slice(30, 40)), kdims=['x', 'y'], crs=ccrs.PlateCarree()).opts(projection=ccrs.PlateCarree())

Meaning something here after if img.crs == proj and np.isclose(src_extent, tgt_extent).all() is the issue:

    def _process(self, img, key=None):
        from cartopy.img_transform import warp_array

        if self.p.fast:
            return self._fast_process(img, key)

        proj = self.p.projection
        x0, x1 = img.range(0, dimension_range=False)
        y0, y1 = img.range(1, dimension_range=False)
        yn, xn = img.interface.shape(img, gridded=True)[:2]
        (px0, py0, px1, py1) = project_extents((x0, y0, x1, y1), img.crs, proj)

        # Some bug in cartopy is causing zero values
        eps = sys.float_info.epsilon
        src_extent = tuple(e+v if e == 0 else e for e, v in
                           zip((x0, x1, y0, y1), (eps, -eps, eps, -eps)))
        tgt_extent = (px0, px1, py0, py1)

        if img.crs == proj and np.isclose(src_extent, tgt_extent).all():
            return img

        arrays = []
        for vd in img.vdims:
            arr = img.dimension_values(vd, flat=False)
            if arr.size:
                projected, _ = warp_array(arr, proj, img.crs, (xn, yn),
                                          src_extent, tgt_extent)
            else:
                projected = arr
            arrays.append(projected)

        if xn == 0 or yn == 0:
            return img.clone([], bounds=tgt_extent, crs=proj)

        xunit = ((tgt_extent[1]-tgt_extent[0])/float(xn))/2.
        yunit = ((tgt_extent[3]-tgt_extent[2])/float(yn))/2.
        xs = np.linspace(tgt_extent[0]+xunit, tgt_extent[1]-xunit, xn)
        ys = np.linspace(tgt_extent[2]+yunit, tgt_extent[3]-yunit, yn)
        return img.clone((xs, ys)+tuple(arrays), bounds=None, kdims=img.kdims,
                         vdims=img.vdims, crs=proj, xdensity=None,
                         ydensity=None)