[FIXED] Cartopy: coastlines not lining up with imshow projection

Issue

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!!

Solution

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!

Answered By - raphael