Cartopy: coastlines not lining up with imshow projection

Question

I'm trying to produce a figure in Python that will line up the map coastlines from Cartopy with a RGB projection using satellite data (POLDER) that is fixed to Sinusoidal grid.

I have tried both Matplotlib's Basemap and Cartopy, and have had better luck with Cartopy, however even after following other people's code, the coastlines do not match up.

What I have so far:

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

proj = ccrs.Sinusoidal(central_longitude=0)
fig = plt.figure(figsize=(12, 12))

# set image extents using the min and max of lon at lat
# image_extent ~= [-73.25, -10, 59.5, 84]
min_lon = np.nanmin(subset_lon)
max_lon = np.nanmax(subset_lon)
min_lat = np.nanmin(subset_lat)
max_lat = np.nanmax(subset_lat)

extents = proj.transform_points(ccrs.Geodetic(), np.array([min_lon, max_lon]), 
                                               np.array([min_lat, max_lat]))
img_extents = [extents[0][0], extents[1][0], extents[0][1], extents[1][1]] 

ax = plt.axes(projection=proj)
# # image RGB
ax.imshow(RGB, origin='upper', extent=img_extents, transform=proj))

ax.set_xmargin(0.05)
ax.set_ymargin(0.10)
ax.coastlines(color='white')

plt.show()

Produces: this figure where the coastlines do NOT match up

Figure with central_lon as -20

Figure using only imshow

I know the projection has to be sinusoidal, so that shouldn't be the issue.

Any ideas on what else it could be or tips on how to fix it?

Here is the dataset: https://drive.google.com/file/d/1vRLKmeAXzCk5cLCJ1syOt7EJnaTmIwOa/view

And code to extract the data and make the image that I would like overlayed with the cartopy coastlines:

data = SD(path_to_file, SDC.READ)
subset_lat = data.select('subset_lat')[:]
subset_lon = data.select('subset_lon')[:]
R = data.select('mband07')[:]
G = data.select('mband02')[:]
B = data.select('mband01')[:]
RGB = np.dstack((R, G, B))
plt.imshow(RGB)

** Edited to add on two comments and code to make imshow RGB image.

Thanks!!

Answer 1

so I had a look at the data and I think the problem actually comes from the fact that your data is NOT provided in an equidistant raster! Therefore using imshow is problematic since it will introduce distortions... (note that imshow will just divide your extent into equally sized pixels based on the number of datapoints, irrespective of the actual coordinates!)

To clarify what I mean... I'm pretty sure the data you have is provided in the grid described here (page 13):

https://www.theia-land.fr/wp-content/uploads/2018/12/parasol_te_level-3_format-i2.00.pdf

I have personally never encountered this kind of grid-definition, but it seems that the number of pixels per latitude is adjusted to account for distortions... to quote from the document above:

Along a parallel, the step is chosen in order to have a resolution as constant as possible. The number of pixels from 180°W to 180°E is chosen equal to 2 x NINT[3240 cos(latitude)] where NINT stands for nearest integer.

To quickly visualize the data I'd therefore suggest to use a scatterplot + the actual coordinates instead of imshow!

---

Alternatively, I'm developing a matplotlib+ cartopy based plotting library (EOmaps) that is intended to simplify visualizations from irregular or non-uniformly gridded datasets such as yours...

After reprojection to Sinusoidal projection, one can see that the distance between the pixel-centers is at least nearly equal...

So we can use this information to visualize the data as rectangles in Sinusoidal projection with a size of 6200x6200... yielding the expected image:

... and zooming in certainly shows that those pixels have approximately the same size but they are definitely not equally distributed!

... here's the code to create the above plot with EOmaps:

from eomaps import Maps
import numpy as np
from pyhdf.SD import SD, SDC

data = SD("Copy of greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf", SDC.READ)

# ---- get the actual coordinates of your datapoints
lon, lat = data.select("subset_lon")[:], data.select("subset_lat")[:]
# mask nan-values
mask = np.logical_and(np.isfinite(lon), np.isfinite(lat))
lon, lat = lon[mask], lat[mask]

# ---- get the colors you want to use
R = data.select("mband07")[:][mask]
G = data.select("mband02")[:][mask]
B = data.select("mband01")[:][mask]
RGB = list(zip(R, G, B)) # a list of rgb-tuples

# ---- plot the data
m = Maps(Maps.CRS.Sinusoidal())
m.add_feature.preset.coastline()
m.set_shape.rectangles(radius=3100, radius_crs=Maps.CRS.Sinusoidal())
# use "dummy" values since we provide explicit colors
m.set_data(data=np.empty(lon.shape), xcoord=lon, ycoord=lat, crs=4326)
m.plot_map(set_extent=True, color=RGB)

... and with this approach also the coastlines are located where they are supposed to be!

Answer 2

The issue consists of two parts:

• You need to compute the extents in the projected space and not in the latitude/longitude space, because there is not a direct correspondence between limit coordinates.
• You need to know the reference longitude of your sinusoidal projection, but unfortunately it is not available in your dataset metadata. A manual inspection showed me that it is somewhere around -14 or -15 degrees.

To solve this by means of basemap, you need a bit of hacking because at the moment Basemap does not support extents for sinusoidal projection (i.e. it is always set to the whole globe). But you can update the extent by hand both in the Basemap object and the underlying Axes object and then it should work:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from pyhdf.SD import SD
from pyhdf.SD import SDC

# Open dataset file.
path = "Greenland_PM01_L2.2007061414150564.058140.1.46-public-release.hdf"
data = SD(path, SDC.READ)

# Read latitude/longitude coordinates.
subset_lat = data.select("subset_lat")[:]
subset_lon = data.select("subset_lon")[:]

# Read RGB image.
R = data.select("mband07")[:]
G = data.select("mband02")[:]
B = data.select("mband01")[:]
RGB = np.dstack((R, G, B))

# Create `Basemap` object to plot the RGB image.
fig1 = plt.figure()
ax1 = fig1.add_axes([0, 0, 1, 1])
bmap1 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax1)

# Set image extents in the projected space. The hack here is that we
# update the corners by hand in both the `Axes` and `Basemap` objects.
x, y = bmap1(subset_lon, subset_lat)
llcrnrx, urcrnrx = np.nanmin(x), np.nanmax(x)
llcrnry, urcrnry = np.nanmin(y), np.nanmax(y)
bmap1.llcrnrx, bmap1.llcrnry = llcrnrx, llcrnry
bmap1.urcrnrx, bmap1.urcrnry = urcrnrx, urcrnry
ax1.set_xlim(llcrnrx, urcrnrx)
ax1.set_ylim(llcrnry, urcrnry)

# Now we can call `imshow`. The keyword argument `origin` does not work
# here, so we can simply reverse the image rows.
bmap1.drawcoastlines(color="white", linewidth=0.5)
bmap1.imshow(RGB[::-1])

# We can also play with the land cover image and check how well the
# coastlines are fitting.
landcover = data.select("subset_land")[:]

fig2 = plt.figure()
ax2 = fig2.add_axes([0, 0, 1, 1])
bmap2 = Basemap(projection="sinu", resolution="i", lon_0=-14, ax=ax2)
bmap2.llcrnrx, bmap2.llcrnry = llcrnrx, llcrnry
bmap2.urcrnrx, bmap2.urcrnry = urcrnrx, urcrnry
ax2.set_xlim(llcrnrx, urcrnrx)
ax2.set_ylim(llcrnry, urcrnry)

bmap2.drawcoastlines(color="red", linewidth=0.5)
bmap2.imshow(landcover[::-1])

which shows this RGB:

and this land cover:

The coastlines match better in the South than in the North, but I am not sure if this is just a limitation of the instrument itself (plus the uncertainty in the correct reference longitude).