preventing spurious horizontal lines for ungridded pcolor(mesh) data

Firstly, thanks for providing some data and a piece of code to reproduce – it meant that I could quickly focus on the issue itself, and not on reproducing the problem.

The major difference between cartopy and basemap is that cartopy can handle vector/raster transformations for you. It is entirely possible to get cartopy to operate in basemap’s fashion where it is beholden on the user to transform their data themselves. The example you have provided is doing precisely this by transforming the lats/lons manually into the target projection. Without a great deal of care, you are going to quickly find antimeridian issues such as the one you have encountered. Thankfully, cartopy has taken a great deal of care with regards to data transformation, and I encourage you to make use of it.

In pseudo code, your code does:

create a mollweide map
convert your lats/lons to mollweide coordinate system
plot newly converted mollweide data on mollweide map

In practice, we want to change the paradigm with cartopy and do:

create a mollweide map
plot lat/lon data on mollweide map

By doing so, we are giving cartopy the necessary context to transform your data correctly.

The major change to your code is to plot the original data (in lats/lons), not the coordinates you transformed by hand:

ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree())

In this instance I’ve used the PlateCarree projection not the Geodetic coordinate system as we don’t currently implement geodetic pcolormesh boxes (i.e. with great circles) and are essentially producing boxes of constant lat/lon.

Using this, we end up producing a plot very similar to the first image in your question, which isn’t exactly what you wanted. The reason for this is that some of the boxes you are defining have a width of ~360 degrees in the PlateCarree projected space (which is a flat piece of paper, and knows nothing of wrap-arounds/the antimeridian).

Let’s take a look at a contrived example. If you think in terms of geodesics, you might expect the following code to produce two small boxes on either side of the map:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=170, maxx=-170, miny=40, maxy=60)

proj = ccrs.Mollweide()

ax = plt.axes(projection=proj)
ax.coastlines()
ax.add_geometries([box], ccrs.PlateCarree(), facecolor="coral", 
                  edgecolor="black", alpha=0.5)

plt.show()

Big box, when we wanted two little ones...

Alas, that is not what we get. This makes sense if we remember that the Plate Carree projection is a 2d cartesian projection where the only valid line between two points is a straight line – it knows nothing about wrapping across the antimeridian.

(it is worth noting: if we were to change the geometry projection to geodetic then we draw great circles between the given points and get the desired boxes)

So to produce the desired boxes we need the box’s coordinates to have a small x range, not one nearing 360 degrees. Thankfully cartopy does allow us to define PlateCarree coordinate values beyond 180 degrees – this is the key to being able to define a PlateCarree box with a small x range.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=170, maxx=190, miny=40, maxy=60)

proj = ccrs.Mollweide()

ax = plt.axes(projection=proj)
ax.coastlines()
ax.add_geometries([box], ccrs.PlateCarree(), facecolor="coral", 
                  edgecolor="black", alpha=0.5)

two boxes - hurrah!

So going back to your example – we have a bunch of lat/lons, which are really defining geodetic patches. Cartopy can’t pcolormesh geodetic coordinates yet – the workaround is to pcolormesh PlateCarree coordinates. Despite geodetic coordinates and PlateCarree coordinate points being interchangeable, they have fundamentally different topology.

In the example you have given it is possible to convert your data to valid PlateCarree topology by adding 360 to the values below 0. Unfortunately this wouldn’t work for geometries that cross the central meridian – that would be a little more involved, and would be a useful extension to cartopy IMO.

The final code now looks like:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

lons = np.array([[-174.719, -175.297, -175.883],
       [-175.164, -175.734, -176.312],
       [-175.594, -176.164, -176.734],
       [-176.016, -176.578, -177.148],
       [-176.43 , -176.984, -177.547],
       [-176.836, -177.383, -177.938],
       [-177.227, -177.773, -178.312],
       [-177.609, -178.148, -178.688],
       [-177.984, -178.516, -179.047],
       [-178.352, -178.875, -179.398],
       [-179.727,  179.766,  179.266],
       [ 179.945,  179.445,  178.945],
       [ 179.625,  179.133,  178.641],
       [ 179.312,  178.828,  178.336],
       [ 179.008,  178.523,  178.039],
       [ 178.711,  178.234,  177.75 ],
       [ 178.414,  177.945,  177.469],
       [ 178.133,  177.656,  177.188],
       [ 177.844,  177.383,  176.914],
       [ 177.57 ,  177.109,  176.648]])

lats = np.array([[ 67.391,  67.492,  67.586],
       [ 67.055,  67.148,  67.25 ],
       [ 66.711,  66.812,  66.906],
       [ 66.375,  66.469,  66.562],
       [ 66.031,  66.125,  66.219],
       [ 65.688,  65.781,  65.875],
       [ 65.344,  65.438,  65.523],
       [ 65.   ,  65.094,  65.18 ],
       [ 64.656,  64.742,  64.836],
       [ 64.312,  64.398,  64.484],
       [ 62.922,  63.   ,  63.086],
       [ 62.57 ,  62.648,  62.734],
       [ 62.219,  62.297,  62.383],
       [ 61.867,  61.945,  62.023],
       [ 61.516,  61.594,  61.672],
       [ 61.164,  61.242,  61.32 ],
       [ 60.812,  60.891,  60.961],
       [ 60.812,  60.891,  60.961],
       [ 60.461,  60.531,  60.609],
       [ 60.102,  60.18 ,  60.25 ]])

data = np.array([[ 231.73,  231.56,  231.22],
       [ 231.72,  231.72,  231.72],
       [ 232.24,  232.73,  233.37],
       [ 233.22,  233.69,  234.01],
       [ 234.33,  234.94,  235.39],
       [ 234.5 ,  235.11,  235.71],
       [ 235.41,  235.71,  236.  ],
       [ 235.27,  235.72,  236.31],
       [ 234.67,  235.43,  235.73],
       [ 235.43,  236.17,  235.88],
       [ 236.18,  236.18,  236.18],
       [ 236.07,  236.36,  236.79],
       [ 235.8 ,  236.1 ,  235.8 ],
       [ 236.84,  236.84,  236.55],
       [ 238.27,  238.27,  238.54],
       [ 237.72,  237.44,  237.72], 
       [ 238.42,  238.28,  238.28],
       [ 238.57,  238.57,  238.43],
       [ 240.17,  240.04,  239.65],
       [ 241.21,  241.21,  241.09]])

proj = ccrs.PlateCarree(central_longitude=180) 
ax = plt.axes(projection=proj)
ax.coastlines('50m')
ax.margins(0.3)

lons[lons < 0] += 360
ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree())

plt.show()

the resulting geometry

If interested, I encourage you to open a cartopy feature request to add a function that generally converts geodetic pcolormesh bounds into plate carree ones. The cartopy tracker can be found at https://github.com/SciTools/cartopy/issues/new.

Leave a Comment