preventing spurious horizontal lines for ungridded pcolor(mesh) data
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?
Solution 1:
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.