Plotting data on faces with pcolormesh (and contourf) fails when data wrap around edge of plot

elenamarie commented 5 years ago

Description

I have model output from the sea ice-ocean model MITgcm which is on a lat-lon-cap grid (similar to a cubed-sphere grid). The model domain can be subdivided into 13 faces with 90 x 90 grid points. When trying to plot the individual faces with cartopy I encounter problems when the data wrap around the edges of the plot. Then the data are smeared across the whole plot from one site to the other. One 'solution' is to change the central longitude of the projection, but then other faces fail when they cross the edge of the plot.

Here a very simplified example:

Plotting with contourf works - except when faces cross the dateline and are on a non-rectangular grid (see the example in the jupyter notebook link below). This is probably related to this issue: #895

Code to reproduce

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm,colors,colorbar
# This is needed because of some error with cartopy and matplotlib axes 
# See also [Axes issue](https://github.com/SciTools/cartopy/issues/1120)
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh
import cartopy.crs as ccrs

X1,Y1 = np.meshgrid(np.linspace(140,229,90),np.linspace(10,70,90))
X2,Y2 = np.meshgrid(np.linspace(-40,49,90),np.linspace(10,70,90))
data1 = np.exp(np.sin(np.deg2rad(X1)) + np.cos(np.deg2rad(Y1)))

fig = plt.figure(figsize=(15, 8))

ax1 = fig.add_subplot(2,2,1, projection=ccrs.PlateCarree(central_longitude=0.))
ax1.stock_img()
ax1.set_global()
ax1.pcolormesh(X1, Y1, data1, transform=ccrs.PlateCarree())
ax1.set_title('Face 1, Pcolormesh, Center 0$\degree$E')

ax2 = fig.add_subplot(2,2,2, projection=ccrs.PlateCarree(central_longitude=180.))
ax2.stock_img()
ax2.set_global()
ax2.pcolormesh(X1, Y1, data1, transform=ccrs.PlateCarree())
ax2.set_title('Face 1, Pcolormesh, Center 180$\degree$E')

ax3 = fig.add_subplot(2,2,3, projection=ccrs.PlateCarree(central_longitude=0.))
ax3.stock_img()
ax3.set_global()
ax3.pcolormesh(X2, Y2, data1, transform=ccrs.PlateCarree())
ax3.set_title('Face 2, Pcolormesh, Center 0$\degree$E')

ax4 = fig.add_subplot(2,2,4, projection=ccrs.PlateCarree(central_longitude=180.))
ax4.stock_img()
ax4.set_global()
ax4.pcolormesh(X2, Y2, data1, transform=ccrs.PlateCarree())
ax4.set_title('Face 2, Pcolormesh, Center 180$\degree$E')

A more detailed example with real data can be found here: jupyter notebook with example

Operating system

macOS Mojave

Cartopy version

0.16.0

conda list (extract)

cartopy                   0.16.0           py36h9263bd1_0  
conda                     4.5.11                py36_1000    conda-forge
ipykernel                 4.10.0                   py36_0  
ipython                   6.5.0                    py36_0  
ipython_genutils          0.2.0            py36h241746c_0
jupyter                   1.0.0                    py36_7  
jupyter_client            5.2.3                    py36_0  
jupyter_console           5.2.0                    py36_1  
jupyter_core              4.4.0                    py36_0
numpy                     1.15.1           py36h6a91979_0  
numpy-base                1.15.1           py36h8a80b8c_0  
python                    3.6.6                hc167b69_0  
scipy                     1.1.0            py36h28f7352_1
xarray                    0.10.9                   py36_0  
xmitgcm                   0.2.2                      py_0    conda-forge

QuLogic commented 5 years ago

# This is needed because of some error with cartopy and matplotlib axes
# See also [Axes issue](https://github.com/SciTools/cartopy/issues/1120)
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh

Your patch here is what's breaking it. You should remove it and downgrade Matplotlib or wait for Cartopy 0.17.0 (planned for the next two weeks.)

elenamarie commented 5 years ago

Thanks for your reply, I should have checked this before.

But downgrading to matplotlib 2.2.2 does not completely solve my problem with data wrapping around the edges of the plot. Due to the geometry of my grid the coordinate arrays of some faces are transposed (an example of the MITgcm lat-lon-cap grid geometry can be found on the ECCO webpage). And with transposed coordinate arrays the data are again smeared across the whole plot when they cross the edge of the plot. Additionally when the x-coordinate array goes from -180 to 180 ( and not 0 to 360) the data are smeared across the plot as well. (Some of my faces cross the 180 meridian and the x-coordinate array go from positive to negative).

Contourf can handle transposed coordinates but fails as well when coordinates go from -180 to 180. Would it then be necessary to shift the coordinates to [0,360]?

Code to reproduce (simplified example, no real data)

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm,colors,colorbar
import cartopy.crs as ccrs

XX,YY = np.meshgrid(np.linspace(140.5,229.5,90),np.linspace(10,70,90))

data1 = np.exp(np.sin(np.deg2rad(XX)) + np.cos(np.deg2rad(YY)))

# X coordinates from 140.5 to 179.5 and -179.5 to -130.5
XX_180 = np.where(XX <= 180, XX, XX-360)

# Transpose the coordinates
XX_T = XX.T
YY_T = YY.T

fig1 = plt.figure(figsize=(15, 12))

ax1 = fig1.add_subplot(3,2,1, projection=ccrs.PlateCarree(central_longitude=0.))
ax1.stock_img()
ax1.set_global()
ax1.pcolormesh(XX, YY, data1, transform=ccrs.PlateCarree())
ax1.set_title('Pcolormesh, Center 0$\degree$E')

ax2 = fig1.add_subplot(3,2,2, projection=ccrs.PlateCarree(central_longitude=180.))
ax2.stock_img()
ax2.set_global()
ax2.pcolormesh(XX, YY, data1, transform=ccrs.PlateCarree())
ax2.set_title('Pcolormesh, Center 180$\degree$E')

ax3 = fig1.add_subplot(3,2,3, projection=ccrs.PlateCarree(central_longitude=0.))
ax3.stock_img()
ax3.set_global()
ax3.pcolormesh(XX_T, YY_T, data1, transform=ccrs.PlateCarree())
ax3.set_title('Transposed coordinates, Pcolormesh, Center 0$\degree$E')

ax4 = fig1.add_subplot(3,2,4, projection=ccrs.PlateCarree(central_longitude=180.))
ax4.stock_img()
ax4.set_global()
ax4.pcolormesh(XX_T, YY_T, data1, transform=ccrs.PlateCarree())
ax4.set_title('Transposed coordinates, Pcolormesh, Center 180$\degree$E')

ax5 = fig1.add_subplot(3,2,5, projection=ccrs.PlateCarree(central_longitude=0.))
ax5.stock_img()
ax5.set_global()
ax5.pcolormesh(XX_180, YY, data1, transform=ccrs.PlateCarree())
ax5.set_title('Coordinates -180 to 180, Pcolormesh, Center 0$\degree$E')

ax6 = fig1.add_subplot(3,2,6, projection=ccrs.PlateCarree(central_longitude=180.))
ax6.stock_img()
ax6.set_global()
ax6.pcolormesh(XX_180, YY, data1, transform=ccrs.PlateCarree())
ax6.set_title('Coordinates -180 to 180, Pcolormesh, Center 180$\degree$E')

fig2 = plt.figure(figsize=(15, 12))

ax1 = fig2.add_subplot(3,2,1, projection=ccrs.PlateCarree(central_longitude=0.))
ax1.stock_img()
ax1.set_global()
ax1.contourf(XX, YY, data1, transform=ccrs.PlateCarree())
ax1.set_title('contourf, Center 0$\degree$E')

ax2 = fig2.add_subplot(3,2,2, projection=ccrs.PlateCarree(central_longitude=180.))
ax2.stock_img()
ax2.set_global()
ax2.contourf(XX, YY, data1, transform=ccrs.PlateCarree())
ax2.set_title('contourf, Center 180$\degree$E')

ax3 = fig2.add_subplot(3,2,3, projection=ccrs.PlateCarree(central_longitude=0.))
ax3.stock_img()
ax3.set_global()
ax3.contourf(XX_T, YY_T, data1, transform=ccrs.PlateCarree())
ax3.set_title('Transposed coordinates, contourf, Center 0$\degree$E')

ax4 = fig2.add_subplot(3,2,4, projection=ccrs.PlateCarree(central_longitude=180.))
ax4.stock_img()
ax4.set_global()
ax4.contourf(XX_T, YY_T, data1, transform=ccrs.PlateCarree())
ax4.set_title('Transposed coordinates, contourf, Center 180$\degree$E')

ax5 = fig2.add_subplot(3,2,5, projection=ccrs.PlateCarree(central_longitude=0.))
ax5.stock_img()
ax5.set_global()
ax5.contourf(XX_180, YY, data1, transform=ccrs.PlateCarree())
ax5.set_title('Coordinates -180 to 180, contourf, Center 0$\degree$E')

ax6 = fig2.add_subplot(3,2,6, projection=ccrs.PlateCarree(central_longitude=180.))
ax6.stock_img()
ax6.set_global()
ax6.contourf(XX_180, YY, data1, transform=ccrs.PlateCarree())
ax6.set_title('Coordinates -180 to 180, contourf, Center 180$\degree$E')
plt.show()

QuLogic commented 5 years ago

I think the situation with transposed coordinates is related to this masking: https://github.com/SciTools/cartopy/blob/3c379ea50459a6b538529523c2bafe14c1c8dfc7/lib/cartopy/mpl/geoaxes.py#L1584-L1588 which appears on first glance to assume that the first dimension is the one that varies vertically.

QuLogic commented 5 years ago

But also, things have been reshaped a bit by that point, so the question is whether this worked in 0.15.0? Because that part was refactored a bit for 0.16.0.

SciTools / cartopy