issue_comments: 953852124

@jklymak Thanks for the explanation.

To get the old behaviour you simply need to do pcolormesh(x, y, Z[:-1, :-1], shading='flat') or make x and y one larger than Z in each dimension and specify the the corners of the quadrilaterals.

Method1: Subset value

This method works for the xarray tutorial data, but not for the TROPOMI polar-orbiting satellite data.

``` %matplotlib inline

import xarray as xr import cartopy.crs as ccrs import matplotlib.pyplot as plt

plt.figure(figsize=(14,6)) ax = plt.axes(projection=ccrs.PlateCarree())

ds = xr.open_dataset('./S5P_OFFL_L2__NO2____20190810T212136_20190810T230306_09456_01_010302_20190816T233944.nc', group='PRODUCT').isel(time=0)

m = ax.pcolormesh(ds['longitude'], ds['latitude'], ds['nitrogendioxide_tropospheric_column'][:-1, :-1], # ds['nitrogendioxide_tropospheric_column'], # shading='auto', transform=ccrs.PlateCarree(), vmin=0, vmax=1e-4, cmap='Spectral_r') ```

(The TROPOMI example data is uploaded to Google Drive)

Method2: bounds

This issue still exists with bounds data: ``` %matplotlib inline import numpy as np import xarray as xr import cartopy.crs as ccrs import matplotlib.pyplot as plt

def prepare_geo(bounds_data): """Prepare lat/lon bounds for pcolormesh. lat/lon bounds are ordered in the following way:: 3----2 | | 0----1 Extend longitudes and latitudes with one element to support "pcolormesh":: (X[i+1, j], Y[i+1, j]) (X[i+1, j+1], Y[i+1, j+1]) +--------+ | C[i,j] | +--------+ (X[i, j], Y[i, j]) (X[i, j+1], Y[i, j+1]) """ # Create the left array left = np.vstack([bounds_data[:, :, 0], bounds_data[-1:, :, 3]]) # Create the right array right = np.vstack([bounds_data[:, -1:, 1], bounds_data[-1:, -1:, 2]]) # Stack horizontally dest = np.hstack([left, right]) # Convert to DataArray dest = xr.DataArray(dest, dims=('y_bounds', 'x_bounds'), attrs=bounds_data.attrs ) return dest

ds = xr.open_dataset('./S5P_OFFL_L2__NO2_20190810T21213620190810T230306_09456_01_010302_20190816T233944.nc', group='PRODUCT').isel(time=0) ds_geo = xr.open_dataset('./S5P_OFFL_L2NO2____20190810T212136_20190810T230306_09456_01_010302_20190816T233944.nc', group='/PRODUCT/SUPPORT_DATA/GEOLOCATIONS').isel(time=0)

lon_bounds = prepare_geo(ds_geo['longitude_bounds']) lat_bounds = prepare_geo(ds_geo['latitude_bounds'])

plt.figure(figsize=(14,6)) ax = plt.axes(projection=ccrs.PlateCarree())

m = ax.pcolormesh(lon_bounds, lat_bounds, ds['nitrogendioxide_tropospheric_column'], # shading='auto', transform=ccrs.PlateCarree(), vmin=0, vmax=1e-4, cmap='Spectral_r') ```