Possibly an issue with get_lonlats() for robinson projection

[ninahakansson]

Code Sample, a minimal, complete, and verifiable piece of code

import matplotlib.pyplot as plt
from pyresample import load_area
import copy
import matplotlib

# Create colormap with high values red
my_cmap = copy.copy(matplotlib.cm.BrBG)
my_cmap.set_over('red', alpha=1)

# Create axis 
fig = plt.figure()
ax, ax2 = fig.subplots(nrows=2)

# Plot ortho longitudes
area_def = load_area('areas.yaml', 'ortho')
lon, lat = area_def.get_lonlats()
ax.imshow(lon,  vmin=-180, vmax=180, cmap=my_cmap)

#Plot robinson longitudes
area_def = load_area('areas.yaml', 'robinson')
lon, lat = area_def.get_lonlats()
ax2.imshow(lon,  vmin=-180, vmax=180, cmap=my_cmap)

# Save or show
plt.savefig('test_robinson.png')
plt.show()

ortho:
  description: Ortho globe
  projection:
    proj: ortho
    lon_0: 40.
    lat_0: -40.
    a: 6370997.0
  shape:
    height: 480
    width: 640
  area_extent:
    lower_left_xy: [-10000000, -10000000]
    upper_right_xy: [10000000, 10000000]

robinson:
  description: robinson
  projection:
    proj: robin
    lat_0: 0.0
    lon_0: 0.0
    ellps: WGS84
  shape:
    height: 480
    width: 640
  area_extent:
    lower_left_xy: [-17900000.0, -8800000.0]
    upper_right_xy: [17900000.0,  8800000.0]

Problem description

I am trying to resample and plot data on a Robinson projection.
However when plotting I get data also in the 4 corners. I have found away around it
and was thinking of adding an example to the docs.
But when using the ortho projection there is no problem with data in the corners.
Longitudes retrieved by get_lonlats() are for the ortho projection very large 10^30
in the corners outside the globe. However for the robinson projection they are not.

Expected Output

I would have expected longitude outside the globe to be very large or nan or nodata
for both the robinson and ortho projection. In the attached picture I expected red color outside the globe.

Actual Result, Traceback if applicable

Versions of Python, package at hand and relevant dependencies

Python 3.6.7
Pyresample 1.13.2

[mraspaud]

Could this be related to OSGeo/PROJ#304 ?

[djhoese]

So the question is is what version of proj and pyproj are you using @ninahakansson and @mraspaud was your fix only for the mollweide projection?

[mraspaud]

Yes, my fix was only for mollweide. Looking at the robinson code, it doesn't seem to be returning any HUGE_VAL when it's out of bounds, as was the case with mollweide before my patch. In other words, it's probably a bug in proj.

[mraspaud]

Ok, here is an example of a problematic behaviour:

> invproj +proj=robin +a=6371000
16986804.16 100000
179d59'21.026"W	0d56'9.769"N
> proj +proj=robin +a=6371000   
179d59'21.026"W	0d56'9.769"N
-16984760.94	99999.99
> invproj +proj=robin +a=6371000
-16984760.94	99999.99
179d59'21.026"W	0d56'9.769"N

[mraspaud]

So one lon/lat has multiple possible projection coordinates...

[mraspaud]

PR made to PROJ: OSGeo/PROJ#2151

[mraspaud]

This is now fixed and merged in PROJ, closing this issue.