Spurious lines of the pcolormesh example

[zxdawn]

What happened:
The example of plotting data using pcolormesh has some spurious lines:

What you expected to happen:
No spurious lines.

Minimal Complete Verifiable Example:

%matplotlib inline
import numpy as np
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
from matplotlib import pyplot as plt

ds = xr.tutorial.open_dataset('rasm').load()

plt.figure(figsize=(14,6))

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()

ds.Tair[0].plot.pcolormesh(ax=ax,
                           transform=ccrs.PlateCarree(),
                           x='xc',
                           y='yc',
                           shading='auto',
                           add_colorbar=False)

ax.coastlines()
ax.set_ylim([0,90]);

Anything else we need to know?:
pcolormesh doc

Environment:

Output of xr.show_versions()

INSTALLED VERSIONS

commit: None
python: 3.8.10 | packaged by conda-forge | (default, Sep 13 2021, 21:46:58)
[GCC 9.4.0]
python-bits: 64
OS: Linux
OS-release: 5.11.0-34-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.12.1
libnetcdf: 4.8.1

xarray: 0.19.0
pandas: 1.3.3
numpy: 1.21.2
scipy: 1.7.1
netCDF4: 1.5.7
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.5.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: 2021.10.0
distributed: 2021.10.0
matplotlib: 3.4.3
cartopy: 0.19.0.post1
seaborn: None
numbagg: None
pint: None
setuptools: 58.0.4
pip: 21.2.4
conda: None
pytest: None
IPython: 7.28.0
sphinx: None

[TomNicholas]

Thanks for flagging this @zxdawn!

😕 We should probably work out whether this is happening because of something we did or something cartopy/matplotlib did...

When did this start? I wonder if there is an easy way to see all the previous editions of xarray's documentation - might help us quickly narrow down which version first had this issue?

[zxdawn]

BTW, the question on StackOverflow, which was raised by @gerritholl a long time ago, looks similar. I'm not sure whether this is the cartopy issue, CC @QuLogic, and @greglucas.

[zxdawn]

@TomNicholas I checked the doc and this issue begins from v0.16.1. Note that there're also small spurious lines after v0.10.9. Before v0.10.9, the figure looks fine. It's may be also related to matplotlib ... CC @jklymak and @timhoffm.

[QuLogic]

Have you tried with the latest Cartopy?

[zxdawn]

@QuLogic Ha, it looks well with the latest cartopy (0.20.1). Thanks a lot.

@TomNicholas So, is it better to keep this open until the doc is updated?

[mathause]

Thanks for the report. This is related to #4364 SciTools/cartopy#1638 and matplotlib/matplotlib#18317.

[jklymak]

Just to explain the issue a bit: pcolormesh(x, y, Z, shading='flat') was the default. If shape(Z) used to be (len(y), len(x)) matplotlib would drop, without warning, the last row and column of Z so that y and x were one larger than Z.

The new default is 'auto' where in the case above it uses shading='nearest' to create new x and y co-ordinates with the dumbest algorithm possible of placing them in the midpoints, and adding two new points on the outside. So rather than drop data from Z, we add to the x and y co-ordinates.

That is usually fine, except the midpoint between 0 and 360 is 180, and we get a weird wrap. Cartopy (xarray) accounted for this with the old shading by adding a NaN. I've not followed what they are doing now.

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.

Sorry this is a pain, but the old behaviour of silently dropping data, while having a long lineage going back to Matlab, was deemed unacceptable. However, if this continues to be a nuisance to folks, I'm sure matplotlib would consider a new shading argument, something like 'flat-drop', or somesuch that would do the data dropping for you.

[zxdawn]

@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____20190810T212136_20190810T230306_09456_01_010302_20190816T233944.nc', group='PRODUCT').isel(time=0)
ds_geo = xr.open_dataset('./S5P_OFFL_L2__NO2____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')

[mathause]

I am closing this as fixed upstream.