Duplicate data when resampling to Eckert IV projection and area spans 180th meridian

[pnuu]

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

from satpy import Scene
import glob

fnames = glob.glob('IMG_DK*')
glbl = Scene(reader='ahi_hrit', filenames=fnames)
glbl.load(['dust'])
lcl = glbl.resample("eckert_iv_global_2km", radius_of_influence=30e3)
lcl.save_datasets(base_dir='/tmp')

And the area definition:

eckert_iv_global_2km:
  description: Eckert IV global area with 2 km resolution
  projection:
    ESRI:54012
  shape:
    height: 8313
    width: 16921
  area_extent:
    lower_left_xy: [-16921201.240, -8245471.428]
    upper_right_xy: [16916798.760, 8250528.572]

Problem description

Some of the data are duplicated outside the valid area.

Expected Output

This image is created with gdalwarp from a full disk image in geos projection with Satpy.

Actual Result, Traceback if applicable

This is the (scaled) output from the above code and area definition.

Note in particular the (duplicate) cloud formation in top-left and top-right regions.

Versions of Python, package at hand and relevant dependencies

Satpy: 0.16.2.dev149+g4313bb4f Pyresample: 1.15.0+2.g3d6d38d Pyproj: 2.4.1 Proj: 6.2.1

[pnuu]

This might explain the duplication: The plot shows longitudes of the created area definition. The corners should not have any coordinates assigned to them.

[mraspaud]

I think this is a problem with several projections in proj. Mollweide had the same problem until I made a pr for proj. We discovered in Copenhagen that the Robin projection also has this problem...

[pnuu]

Yeah, I suspected this is something more general than just the Eckert IV. I'll see if I can hack something together.

[pnuu]

NOTE: this used incorrect syntax (should be inverse=True)

Interestingly, pyproj.Proj seems to be working as expected:

from pyproj import Proj
prj = Proj("ESRI:54012")
print(prj(-16921202.922943164, -8313010.558165222, inv=True))
(1e+30, 1e+30)

The inputs are the bottom left corner extents.

[djhoese]

I would expect Proj and GDAL to behave the same since (I think) they are both using PROJ underneath.

[pnuu]

Indeed, both https://github.com/OSGeo/PROJ/pull/304 and https://github.com/pytroll/pyresample/issues/232 seem to be the same case as this. I think I have a fix (the same as in the PROJ patch @mraspaud made), but don't have a way to test test with Satpy/Pyresample to verify.

[pnuu]

And the pyproj example above has an incorrect kwarg. It should be inverse=True, and then we get an incorrect result. So the bug is in PROJ.