Search code examples
cartopy

Cartopy fails to correctly contour data on rotated grid

when contouring a dataset that is defined on a rotated pole grid, I get the following result:

enter image description here

This is only a problem when using contour and contourf, not pcolormesh.

Anyone have any idea what could be going on here? Should I file a bug report, and if yes, should I do so with cartopy or with matplotlib?

Reproducing

Solution

  • Data for contouring should be split into 2 parts to avoid errors you found. I choose longitude=0 as the dividing line. Numpy's masked array technique is used to achieve the data manipulation. Here is the working code that produces a useful plot.

    import numpy as np
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy.ma as ma
from netCDF4 import Dataset

nc = Dataset('./data/snow_rlat_rlon.nc')

# Prep values for contouring
snow_2d_array = nc.variables[u'snowfall'][:]   # need *(86400*365); mm/s-> mm/yr
lat_2d_array = nc.variables[u'lat2d'][:]
lon_2d_array = nc.variables[u'lon2d'][:]

# do masked-array on the lon_2d
lon2d_greater = ma.masked_greater(lon_2d_array, -0.01)
lon2d_lesser = ma.masked_less(lon_2d_array, 0)

# apply masks to other associate arrays: lat_2d
lat2d_greater = ma.MaskedArray(lat_2d_array, mask=lon2d_greater.mask)
lat2d_lesser = ma.MaskedArray(lat_2d_array, mask=lon2d_lesser.mask)
# apply masks to other associate arrays: snow_2d
snow_2d_greater = ma.MaskedArray(snow_2d_array, mask=lon2d_greater.mask)
snow_2d_lesser = ma.MaskedArray(snow_2d_array, mask=lon2d_lesser.mask)

# set levels for contouring of snow_2d
levels = (0, 25, 50, 75, 100, 200, 400, 600, 800, 1000, 2000, 4000)

# get snow_2d value-limits for use with colormap
vmax, vmin = snow_2d_array.max()*86400*365, snow_2d_array.min()*86400*365
cmap1 = "viridis"
norm1 = colors.BoundaryNorm(boundaries=levels, ncolors=16)
norm2 = colors.Normalize(vmin=vmin, vmax=vmax/4)

# setup fig+axes, specifying projection to use
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.SouthPolarStereo()})
fig.set_size_inches([10, 10])

ax.coastlines(color="red", linewidth=2)  # draw coastlines in red

# plot contour using each part of the 2 masked data sets
ct1 = ax.contour(lon2d_greater, lat2d_greater, snow_2d_greater*86400*365, \
                 norm=norm2, levels=levels, \
                 transform=ccrs.PlateCarree())

ct2 = ax.contour(lon2d_lesser, lat2d_lesser, snow_2d_lesser*86400*365, \
                 norm=norm2, levels=levels, \
                 transform=ccrs.PlateCarree()) 

#plt.colorbar(ct1, shrink=0.85)
plt.show()

    The output plot:

    spole_plot

    For filled-contour, replace ax.contour() with ax.contourf() and add:

    ax.set_xlim((-4052327.4304452268, 4024164.250636036))

    in front of plt.show().

    contourf

    Hope it is useful to your project.