area not equal after storing and loading

[gerritholl]

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

from pyproj import CRS
from pyresample import create_area_def
from pyresample.area_config import load_area_from_string
ar = create_area_def("test", CRS.from_authority("IAU_2015", 39916), area_extent=[-10_000_000, -10_000_000, 10_000_000, 10_000_000], resolution=10000)
ar2 = load_area_from_string(ar.dump())
print(ar == ar2)

Problem description

After storing and loading, the areas should be equal. Even considering floating point precision.

Expected Output

True

Actual Result, Traceback if applicable

/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/pyproj/crs/crs.py:1256: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  return self._crs.to_proj4(version=version)
False

Versions of Python, package at hand and relevant dependencies

pyresample v1.23.0.

[djhoese]

Can you show what the .dump() output is?

[gerritholl]

ar.dump() and ar2.dump() are actually equal:

test:
  description: test
  projection:
    proj: eqc
    lat_ts: 0
    lat_0: 0
    lon_0: 180
    x_0: 0
    y_0: 0
    a: 6378136.6
    rf: 298.257006177324
    no_defs: null
    type: crs
  shape:
    height: 2000
    width: 2000
  area_extent:
    lower_left_xy: [-10000000.0, -10000000.0]
    upper_right_xy: [10000000.0, 10000000.0]
    units: m

test:
  description: test
  projection:
    proj: eqc
    lat_ts: 0
    lat_0: 0
    lon_0: 180
    x_0: 0
    y_0: 0
    a: 6378136.6
    rf: 298.257006177324
    no_defs: null
    type: crs
  shape:
    height: 2000
    width: 2000
  area_extent:
    lower_left_xy: [-10000000.0, -10000000.0]
    upper_right_xy: [10000000.0, 10000000.0]
    units: m

[gerritholl]

It seems information about the CRS has gone missing:

In [255]: ar.crs
Out[255]:
<Derived Projected CRS: IAU_2015:39916>
Name: Earth (2015) / Ographic / Equirectangular, clon = 180
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: Equirectangular, clon = 180
- method: Equidistant Cylindrical
Datum: Earth (2015)
- Ellipsoid: Earth (2015)
- Prime Meridian: Reference Meridian

In [256]: ar2.crs
Out[256]:
<Derived Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unk ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Equidistant Cylindrical
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich

This also impacts files saved, at least with the ninjogeotiff writer:

 diff -u <(gdalinfo 202104100050-GOES-17-abi-test-colorized_ir_clouds.tif) <(gdalinfo 202104100050-GOES-17-abi-test-colorized_ir_clouds-2.tif)
--- /dev/fd/63  2022-04-08 16:51:19.588313642 +0200
+++ /dev/fd/62  2022-04-08 16:51:19.588313642 +0200
@@ -1,14 +1,14 @@
 Driver: GTiff/GeoTIFF
-Files: 202104100050-GOES-17-abi-test-colorized_ir_clouds.tif
+Files: 202104100050-GOES-17-abi-test-colorized_ir_clouds-2.tif
 Size is 2000, 2000
 Coordinate System is:
-PROJCRS["Earth (2015) / Ographic / Equirectangular, clon = 180",
-    BASEGEOGCRS["Earth (2015) / Ographic",
-        DATUM["Earth (2015)",
-            ELLIPSOID["Earth (2015)",6378136.6,298.257006177324,
+PROJCRS["unknown",
+    BASEGEOGCRS["unknown",
+        DATUM["unknown",
+            ELLIPSOID["unknown",6378136.6,298.257006177324,
                 LENGTHUNIT["metre",1,
                     ID["EPSG",9001]]]],
-        PRIMEM["Reference Meridian",0,
+        PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433,
                 ID["EPSG",9122]]]],
     CONVERSION["Equidistant Cylindrical",

[djhoese]

So I think this is one of those cases where PROJ.4 really does lose information from the original authority CRS description. This is one of the reasons (I'm guessing) that PROJ no longer consideres PROJ.4 strings/dicts the best way to describe a projection. I think you have a real case here where we need to change dump to use WKT instead of a PROJ.4 dict. If not by default then at least as an option.

[gerritholl]

Is there an issue yet for changing dump to use WKT instead of a proj4 dict?

[djhoese]

I think this is the best/closest one: https://github.com/pytroll/pyresample/issues/335