Example suggestion: Remap data to Robinson projection starting with cartopy area

[ninahakansson]

Here is an example of remapping data to Robinson projection. Using the cartopy definition to get the border around the globe correct. A pyresample area definition is created from the cartopy definition, for remapping. There is also a fix to remove data outside the globe borders, which is most likely caused by a bug in proj or maybe pyproj. I am hoping to get time to turn it into a notebook.

import matplotlib.pyplot as plt
from pyresample import geometry as prg
import numpy as np
from pyresample.kd_tree import resample_nearest
import copy
import matplotlib
from pyresample import load_area

# Define area from cartopy to get globe boarders
crs = ccrs.Robinson()
crs.bounds = (crs.x_limits[0], crs.x_limits[1], 
              crs.y_limits[0], crs.y_limits[1])

# Create pyresample area_def object for resampling
area_def = prg.AreaDefinition('robinson',
                              'robinson',
                              'robinson',
                              projection=crs.proj4_params,
                              width=1000, height=500, 
                              area_extent=(crs.x_limits[0],
                                           crs.y_limits[0],
                                           crs.x_limits[1],
                                           crs.y_limits[1]))

# Make test data
xi = np.linspace(-179, 179, 1000)
yi = np.linspace(-90, 90, 500)
lons, lats = np.meshgrid(xi, yi)
data = lons*6.0 + lats*3.0

# Remap to robinson projection
swath_def = prg.SwathDefinition(lons=lons, lats=lats)
result = resample_nearest(
    swath_def, data, area_def,
    radius_of_influence=500*1000*2.5, fill_value=None)

# Create colormap with nan's not visible
my_cmap = copy.copy(matplotlib.cm.BrBG)
my_cmap.set_bad('1.0', alpha=0)
my_cmap.set_over('red', alpha=1)

# Hack to remove data in corners
lon, lat = area_def.get_lonlats()
yshape, xshape = result.data.shape
result.data[:,0:xshape//2][lon[:,0:xshape//2]>0] = np.nan
result.data[:,xshape//2:][lon[:,xshape//2:]<0] = np.nan

# Plot image with repeated data in corners masked
fig = plt.figure()
ax = fig.subplots(nrows=1, subplot_kw={'projection': crs})
ax.coastlines()
ax.imshow(result, transform=crs, extent=crs.bounds, cmap=my_cmap)
plt.savefig('test_robinson.png',  bbox_inches='tight')
plt.show()

[djhoese]

Interesting. Do you have a screenshot of what the output looks like? Maybe any existing cartopy examples should be changed to use the work you've done here or have this example updated to include your example? I'd prefer to not include code that is working around a PROJ bug if possible.

[ninahakansson]

If changing to the ortho projection, the hack is not needed. However I really needed to remap and plot data at the Robinson projection. I had the same problem 4 years ago and at the time I had to settle with a bad solution. Without the hack, the main thing is that starting with the Robinson projection and turning it into a pyresample area definition, gives the borders of the globe also plotted in the image.