preventing spurious horizontal lines for ungridded pcolor(mesh) data

Question

When I have a stretch of ungridded lat/lon/data pairs that cross the antimeridian, such that longitudes swap from -180 to +180, how can I prevent cartopy with pcolor(mesh) from drawing grid cells filling the entire globe? My problem is identical to the one here, except I'm using cartopy rather than basemap. A nearly 5 year old comment to the linked question (which is about basemap) claims there is a cartopy solution but such has not been posted.

Example code:

#!/usr/bin/env python3.6

import numpy
import matplotlib.pyplot
import cartopy.crs

lons = numpy.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 = numpy.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 = numpy.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 = cartopy.crs.Mollweide() 
ax = matplotlib.pyplot.axes(projection=proj)
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats)
ax.coastlines()
ax.pcolormesh(trans[:, :, 0], trans[:, :, 1], data, transform=proj)

matplotlib.pyplot.savefig("/tmp/test.png")

Expected output would be a map with a bit of data centred somewhere in the North Pacific Ocean. In reality, I get a very elongated map spanning the width of the entire Earth:

I've limited the data to a small number of points such that I can more easily incorporate it in the question, but in reality I have a full orbit of polar satellite data that always crosses both poles, and therefore always crosses the antimeridian. The result for a real orbit may look like this:

Changing the central longitude relocates the problem. I can reduce the severity by choosing the central longitude away from where I cross the map edge. In this example, the same data as in the previous map are plotted but with a central longitude of 90°E:

This pull request from 2012 appears related, so apparently there is supposed to be a related feature, but I have no clue how to use it. The problem appears with any global map projection. I'm using cartopy 0.15.1.

How can I plot this correctly?

user_741316 · Answer

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()

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)

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()

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.

user_4124317 · Answer

Edit: Note that there are two follow-up questions on this:

why does pcolor with masked array still fill quadrangles connecting to masked points, and how do I stop this? and
Why does pcolor with masked arrays fill undesired quadrangles when projected in cartopy coordinates?

which would in total show a workaround for the given problem**

Original answer:

I'm not sure if the following is of great help, but let's try. If you are able to split your data at the position where the shift between -180 and 180 occurs, you could plot two different pcolor plots and would in this way prevent the pcolorshapes to revolve once around the earth.

Using the test data, this is rather easy; I just changed the sign of one of the data points. But one could also leave one line of data out or so. Then plotting the data seperately gives the desired result. Note that one would need a normalization to make sure the colormapping is synchronuous for both plots.

norm = plt.Normalize(data.min(), data.max())
    ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm)
    ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm)

.....

Complete code:

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

lons = numpy.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.999,  179.766,  179.266], #<- changed sign here
       [ 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 = numpy.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 = numpy.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]])
print lons.shape, lats.shape, data.shape
proj = cartopy.crs.Mollweide() 
ax = plt.axes(projection=proj)
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats)
ax.coastlines()
norm = plt.Normalize(data.min(), data.max())
ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm)
ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm)

plt.show()

preventing spurious horizontal lines for ungridded pcolor(mesh) data

2 Answers