[FIXED] Matplotlib Contour/Contourf in Cartopy singularity at North Pole/ dateline

Issue

I am trying to produce heatmaps showing atmospheric attenuation values for a RF link to a satellite above the North Pole, but I have issues with the interpolation done by the Matplotlib contour/contourf functions.

The linear interpolation done by the contourf function does not work well around the N.Pole, as I suspect it does not know to interpolate between values which go from (-180 deg to +180 deg) - i.e. cross the dateline, or cross the pole.

Any suggestions on a different approach to generate the heatmap, to avoid this horrible hole at the centre?!

Code below to generate plot.

import cartopy.crs as ccrs
import cartopy.feature


plt.figure(figsize=(10,10))

# Initialise Cartopy Axes.
proj=ccrs.LambertAzimuthalEqualArea(central_longitude=0, central_latitude=90)

ax = plt.axes(projection = proj)

ax.set_extent([-180,180,45,90], ccrs.PlateCarree())

ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.gridlines(ls=":",color="grey",lw=0.5)


x0,x1 = attenuation_df.lon.min(), attenuation_df.lon.max()
y0,y1 = attenuation_df.lat.min(), attenuation_df.lat.max()

x,y = np.linspace(x0,x1,1000), np.linspace(y0,y1,1000)
X,Y = np.meshgrid(x,y)

Z = scipy.interpolate.griddata(
    attenuation_df[["lon","lat"]],
    attenuation_df["attenuation"],
    (X,Y),
    method="linear",
)

plt.contourf(X,Y,Z,transform=ccrs.PlateCarree(),alpha=0.5)
plt.colorbar(shrink=0.5)

plt.title("Attenuation")

plt.show()

Attenuation_df is a Pandas Dataframe which contains an attenuation value at approximately 3500 sample points, which are equally spaced around the globe. Here is the location of the sample points:

Here is the header of attenuation_df:

	lon	lat	attenuation
0	-30.8538	48.8813	0.860307
1	-29.0448	49.5026	0.783662
2	-27.2358	50.1317	0.720165
3	-32.6628	48.2676	0.947662
4	37.4226	46.0322	0.27495

The link to the csv of attenuation_df is here: https://pastebin.com/NYA1jFgt

Solution

A solution is to reproject your data to a different coordinate system, my suggestion is to use a Polar Stereographic system. However, the large "hole" centered at the North Pole is not coming from the coordinate system in use but to the presence of some nans in your dataset, so you first have to remove those values.

Here a working solution:

from pyproj import Proj

# Define a pyproj function to reproject data
def coordinate_conv(x, y, inverse = True):
    p = Proj('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs')
    return p(x, y, inverse = inverse)

# Drop null values
attenuation_df.dropna(how = 'any', inplace = True)

# Reproject data
rpjx, rpjy = coordinate_conv(attenuation_df.lon, attenuation_df.lat, False)
rpj_cord = pd.DataFrame({'x': rpjx, 'y': rpjy})

# Interpoolate data
x,y = np.linspace(rpjx.min(),rpjx.max(),1000), np.linspace(rpjy.min(),rpjy.max(),1000)
X,Y = np.meshgrid(x,y)

Z = interpolate.griddata(
    rpj_cord,
    attenuation_df["attenuation"],
    (X,Y),
    method="linear",
)

# Figure
plt.figure(figsize=(10,10))

# Initialise Cartopy Axes.
proj=ccrs.LambertAzimuthalEqualArea(central_longitude=0, central_latitude=90)

ax = plt.axes(projection = proj)

ax.set_extent([-180,180,45,90], ccrs.PlateCarree())

ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.gridlines(ls=":",color="grey",lw=0.5)

kw = dict(central_latitude=90, central_longitude=-45, true_scale_latitude=70)
plt.contourf(X,Y,Z, transform=ccrs.Stereographic(**kw),alpha=0.5)
plt.colorbar(shrink=0.5)

plt.title("Attenuation")

And this is the output figure:

Answered By - baccandr