Closed avalentino closed 6 months ago
The problem seems to be linked to the missing proj_id
in the AreaDefinition
loaded form the yaml file.
We don't currently have PROJ 9.3.0 available to us through conda-forge. Given that, would it be possible for you to get a full print out of the comparison? Or I guess modify self.maxDiff
as the message recommends? The proj_id
, unless I'm missing something, shouldn't be dependent on the PROJ version so I'm a bit confused.
$ python3 -m pytest -k test_def2yaml_converter
================================================================== test session starts ==================================================================
platform linux -- Python 3.11.5, pytest-7.4.0, pluggy-1.2.0
benchmark: 3.2.2 (defaults: timer=time.perf_counter disable_gc=False min_rounds=5 min_time=0.000005 max_time=1.0 calibration_precision=10 warmup=False warmup_iterations=100000)
rootdir: /home/antonio/debian/git/pyresample
plugins: benchmark-3.2.2, mock-3.11.1, remotedata-0.4.0, httpbin-1.0.0, doctestplus-1.0.0, hypothesis-6.82.6, lazy-fixture-0.6.3, recording-0.13.0, cov-4.1.0, requests-mock-1.9.3, timeout-2.1.0, xdist-3.3.1
collected 1053 items / 1052 deselected / 1 selected
pyresample/test/test_utils.py F [100%]
======================================================================= FAILURES ========================================================================
___________________________________________________________ TestMisc.test_def2yaml_converter ____________________________________________________________
self = <pyresample.test.test_utils.TestMisc testMethod=test_def2yaml_converter>
def test_def2yaml_converter(self):
self.maxDiff = None
import tempfile
from pyresample import convert_def_to_yaml, parse_area_file
def_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg')
filehandle, yaml_file = tempfile.mkstemp()
os.close(filehandle)
try:
convert_def_to_yaml(def_file, yaml_file)
areas_new = set(parse_area_file(yaml_file))
areas = parse_area_file(def_file)
areas_old = set(areas)
areas_new = {area.area_id: area for area in areas_new}
areas_old = {area.area_id: area for area in areas_old}
> self.assertEqual(areas_new, areas_old)
E AssertionError: {'ease_nh': Area ID: ease_nh
E Description: Arctic[1247 chars]714)} != {'ease_sh': Area ID: ease_sh
E Description: Antarc[1344 chars]625)}
E {'ease_nh': Area ID: ease_nh
E Description: Arctic EASE grid
E + Projection ID: ease_nh
E Projection: {'R': '6371228', 'lat_0': '90', 'lon_0': '0', 'no_defs': 'None', 'proj': 'laea', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
E Number of columns: 425
E Number of rows: 425
E Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625),
E 'ease_sh': Area ID: ease_sh
E Description: Antarctic EASE grid
E + Projection ID: ease_sh
E Projection: {'R': '6371228', 'lat_0': '-90', 'lon_0': '0', 'no_defs': 'None', 'proj': 'laea', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
E Number of columns: 425
E Number of rows: 425
E Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625),
E 'ortho': Area ID: ortho
E Description: Ortho globe
E + Projection ID: ortho_globe
E Projection: {'ellps': 'sphere', 'lat_0': '-40', 'lon_0': '40', 'no_defs': 'None', 'proj': 'ortho', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
E Number of columns: 640
E Number of rows: 480
E Area extent: (-10000000.0, -10000000.0, 10000000.0, 10000000.0),
E 'pc_world': Area ID: pc_world
E Description: Plate Carree world map
E + Projection ID: pc_world
E Projection: {'datum': 'WGS84', 'lat_0': '0', 'lat_ts': '0', 'lon_0': '0', 'no_defs': 'None', 'proj': 'eqc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
E Number of columns: 640
E Number of rows: 480
E Area extent: (-20037508.3428, -10018754.1714, 20037508.3428, 10018754.1714)}
pyresample/test/test_utils.py:232: AssertionError
=================================================================== warnings summary ====================================================================
pyresample/test/test_spherical_geometry.py:25
/home/antonio/debian/git/pyresample/pyresample/test/test_spherical_geometry.py:25: DeprecationWarning: This module will be removed in pyresample 2.0, please use the `pyresample.spherical` module functions and class instead.
from pyresample.spherical_geometry import Arc, Coordinate
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter
/home/antonio/debian/git/pyresample/pyresample/area_config.py:851: DeprecationWarning: 'create_areas_def' is deprecated. Please use `dump` instead, which also supports writing directly to a file.
yaml_file.write(area.create_areas_def())
-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html
================================================================ short test summary info ================================================================
FAILED pyresample/test/test_utils.py::TestMisc::test_def2yaml_converter - AssertionError: {'ease_nh': Area ID: ease_nh
==================================================== 1 failed, 1052 deselected, 7 warnings in 2.42s =====================================================
I've been looking at this and I am very confused. Based on the comparison done in the test (new first, old second) I would assume that the +
in the diff means that the Projection ID
exists in the old-format parsed files, but not in the new YAML generated from the old .cfg
format.
As far as I can tell proj_id
was never produced in the conversion to the YAML format. In fact, I downgraded my local copy of pyresample and ran this test with the debugger and indeed areas_old
has a .proj_id
, but areas_new
does not. So then that must mean something has changed in the comparison code...great...
Indeed the AreaDefinition comparison doesn't even look at the proj_id
:
However, since you've identified PROJ as being the cause for this error, I wonder if the CRSes are considered not equal.
Wait...how do you have PROJ 9.3.0? That isn't even completely released yet?
It is not currently listed here: https://proj.org/en/latest/download.html
Yes, it is an RC
Ok, that makes it a little more difficult for me to mix in with the rest of my dev environment. If you have any way to get a python interpreter with the PROJ RC and pyproj, then I think we could try a few things to be sure. I could maybe also look at the PROJ release notes and guess what might be breaking things. Bottom line is that pyresample is doing some pyproj CRS
-> PROJ.4 string conversion and then doing PROJ.4 string -> CRS conversion for the "new" YAML based area definitions.
So the original starting .cfg
version of the ease_nh
projection is:
proj=laea, lat_0=90, lon_0=0, a=6371228.0, units=m
But I see in the error message you posted that it says:
{'R': '6371228', 'lat_0': '90', 'lon_0': '0', 'no_defs': 'None', 'proj': 'laea', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
So we can clearly see the PROJ.4 isn't exactly the same, but that's not always a problem as PROJ figures it out anyway. I just ran it with my environment and these are equal, see:
In [1]: from pyproj import CRS
In [2]: crs1 = CRS.from_proj4("+proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m")
In [3]: crs1.to_proj4()
Out[3]: '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
In [4]: crs2 = CRS.from_user_input(crs1.to_proj4())
proj = self._crs.to_proj4(version=version)
In [5]: crs1 == crs2
Out[5]: True
>>> from pyproj import CRS
>>> crs1 = CRS.from_proj4("+proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m")
>>> crs1.to_proj4()
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
'+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
>>> crs2 = CRS.from_user_input(crs1.to_proj4())
>>> crs1 == crs2
True
>>> import pyproj
>>> pyproj.show_versions()
pyproj info:
pyproj: 3.6.0
PROJ: 9.2.1
data dir: /usr/share/proj
user_data_dir: /home/antonio/.local/share/proj
PROJ DATA (recommended version): 1.15
PROJ Database: 1.2
EPSG Database: v10.094 [2023-08-08]
ESRI Database: ArcGIS Pro 3.1 [2023-19-01]
IGNF Database: 3.1.0 [2019-05-24]
System:
python: 3.11.5 (main, Aug 29 2023, 15:31:31) [GCC 13.2.0]
executable: /usr/bin/python3
machine: Linux-6.2.0-31-generic-x86_64-with-glibc2.37
Python deps:
certifi: 2022.9.24
Cython: 0.29.36
setuptools: 68.1.2
pip: 23.2
Well wait, that says PROJ 9.2.1. I thought the problem was with 9.3.0?
Yep.
I confirm that the libproj version installed for the test is v9.3.0-rc1.
Apparently the string printed by the show_versions
refers to the libproj version used at build time.
I don't know if there is a way to query the version of the library used at runtime.
Ok, well with your True
result in my little test I'm a little confused (again) at why this is failing in your test execution.
My only other test that I can think of is to do the same test I gave you, but start with the other projection information in the test areas.cfg. So I guess:
"+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m"
"+proj=eqc"
"+proj=ortho +lon_0=40 +lat_0=-40 +ellps=sphere"
If these all say they are equal when they go through this round trip test then...I don't know.
Dear @djhoese, for all the proj strings that you suggest the result is always crs1 == crs2
.
proj-string: +proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
csr1: +proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs
crs1 == crs2: True
proj-string: +proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m
csr1: +proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs
crs1 == crs2: True
proj-string: +proj=eqc
csr1: +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs
crs1 == crs2: True
proj-string: +proj=ortho +lon_0=40 +lat_0=-40 +ellps=sphere
csr1: +proj=ortho +lat_0=-40 +lon_0=40 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs +type=crs
crs1 == crs2: True
By the way I noticed the the yaml files written in the two cases is different:
$ cat areas-libproj9.2.1.yaml
ease_sh:
description: Antarctic EASE grid
projection:
proj: laea
lat_0: -90
lon_0: 0
x_0: 0
y_0: 0
R: 6371228
no_defs: null
type: crs
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
units: m
ease_nh:
description: Arctic EASE grid
projection:
proj: laea
lat_0: 90
lon_0: 0
x_0: 0
y_0: 0
R: 6371228
no_defs: null
type: crs
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
units: m
pc_world:
description: Plate Carree world map
projection:
EPSG: 4087
shape:
height: 480
width: 640
area_extent:
lower_left_xy: [-20037508.342789244, -10018754.171394622]
upper_right_xy: [20037508.342789244, 10018754.171394622]
ortho:
description: Ortho globe
projection:
proj: ortho
lat_0: -40
lon_0: 40
x_0: 0
y_0: 0
ellps: sphere
no_defs: null
type: crs
shape:
height: 480
width: 640
area_extent:
lower_left_xy: [-10000000.0, -10000000.0]
upper_right_xy: [10000000.0, 10000000.0]
units: m
$ cat areas-libproj9.3.yaml
ease_sh:
description: Antarctic EASE grid
projection:
EPSG: 3409
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
ease_nh:
description: Arctic EASE grid
projection:
EPSG: 3408
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
pc_world:
description: Plate Carree world map
projection:
EPSG: 4087
shape:
height: 480
width: 640
area_extent:
lower_left_xy: [-20037508.342789244, -10018754.171394622]
upper_right_xy: [20037508.342789244, 10018754.171394622]
ortho:
description: Ortho globe
projection:
proj: ortho
lat_0: -40
lon_0: 40
x_0: 0
y_0: 0
ellps: sphere
no_defs: null
type: crs
shape:
height: 480
width: 640
area_extent:
lower_left_xy: [-10000000.0, -10000000.0]
upper_right_xy: [10000000.0, 10000000.0]
units: m
$ diff -u areas-libproj9.2.1.yaml areas-libproj9.3.yaml
--- areas-libproj9.2.1.yaml 2023-08-31 15:14:31.686634333 +0000
+++ areas-libproj9.3.yaml 2023-08-31 15:10:20.858926002 +0000
@@ -1,39 +1,23 @@
ease_sh:
description: Antarctic EASE grid
projection:
- proj: laea
- lat_0: -90
- lon_0: 0
- x_0: 0
- y_0: 0
- R: 6371228
- no_defs: null
- type: crs
+ EPSG: 3409
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
- units: m
ease_nh:
description: Arctic EASE grid
projection:
- proj: laea
- lat_0: 90
- lon_0: 0
- x_0: 0
- y_0: 0
- R: 6371228
- no_defs: null
- type: crs
+ EPSG: 3408
shape:
height: 425
width: 425
area_extent:
lower_left_xy: [-5326849.0625, -5326849.0625]
upper_right_xy: [5326849.0625, 5326849.0625]
- units: m
pc_world:
description: Plate Carree world map
projection:
Not sure if this is relevant or not.
Oh very very interesting. That must be it then.
In [1]: from pyproj import CRS
In [2]: crs1 = CRS.from_proj4("+proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m")
In [3]: crs2 = CRS.from_epsg(3408)
In [4]: crs1.to_proj4()
/home/davidh/miniconda3/envs/satpy_py311_2/lib/python3.11/site-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[4]: '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
In [5]: crs2.to_proj4()
/home/davidh/miniconda3/envs/satpy_py311_2/lib/python3.11/site-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[5]: '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
In [6]: crs1.to_epsg()
In [7]: crs2.to_epsg()
Out[7]: 3408
In [8]: crs1 == crs2
Out[8]: False
In [9]: crs1 == CRS.from_proj4(crs2.to_proj4())
/home/davidh/miniconda3/envs/satpy_py311_2/lib/python3.11/site-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[9]: True
So...I guess I would maybe consider this a bug in PROJ 9.3?
When you run the above code with 9.3, what do you get when just printing print(repr(crs1))
? In 9.2.x I get a lot more information from the crs2
EPSG-based one (note the below output is from the -90 latitude version of the CRS):
In [9]: crs1
Out[9]:
<Projected CRS: +proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units ...>
Name: unknown
Axis Info [cartesian]:
- E[north]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Lambert Azimuthal Equal Area (Spherical)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich
In [10]: crs2
Out[10]:
<Projected CRS: EPSG:3409>
Name: NSIDC EASE-Grid South
Axis Info [cartesian]:
- X[north]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: Southern hemisphere.
- bounds: (-180.0, -90.0, 180.0, 0.0)
Coordinate Operation:
- name: US NSIDC Equal Area south projection
- method: Lambert Azimuthal Equal Area (Spherical)
Datum: Not specified (based on International 1924 Authalic Sphere)
- Ellipsoid: International 1924 Authalic Sphere
- Prime Meridian: Greenwich
The two crs are considered equal, which seems to be the correct behaviour, but the representation is still different (like in libproj v9.2.1).
In [1]: from pyproj import CRS
In [2]: crs1 = CRS.from_proj4("+proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m")
In [3]: crs2 = CRS.from_epsg(3408)
In [4]: crs1.to_proj4()
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[4]: '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
In [5]: crs2.to_proj4()
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[5]: '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs'
In [6]: crs1.to_epsg()
Out[6]: 3408
In [7]: crs2.to_epsg()
Out[7]: 3408
In [8]: crs1 == crs2
Out[8]: False
In [9]: crs1 == CRS.from_proj4(crs2.to_proj4())
/usr/lib/python3/dist-packages/pyproj/crs/crs.py:1293: 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
proj = self._crs.to_proj4(version=version)
Out[9]: True
In [10]: print(repr(crs1))
<Projected CRS: +proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units= ...>
Name: unknown
Axis Info [cartesian]:
- E[south]: Easting (metre)
- N[south]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Lambert Azimuthal Equal Area (Spherical)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich
In [11]: print(repr(crs2))
<Projected CRS: EPSG:3408>
Name: NSIDC EASE-Grid North
Axis Info [cartesian]:
- X[south]: Easting (metre)
- Y[south]: Northing (metre)
Area of Use:
- name: Northern hemisphere.
- bounds: (-180.0, 0.0, 180.0, 90.0)
Coordinate Operation:
- name: US NSIDC Equal Area north projection
- method: Lambert Azimuthal Equal Area (Spherical)
Datum: NSIDC International 1924 Authalic Sphere
- Ellipsoid: International 1924 Authalic Sphere
- Prime Meridian: Greenwich
Ok, but now in PROJ 9.3 (as you showed with your diff) the to_epsg()
produces an actual EPSG code. Pyresample assumes that the EPSG code is more accurate so it uses that when converting the AreaDefinition (the one loaded from PROJ.4 string in the .cfg
file) to the YAML file. The CRS
that is then loaded from YAML is equivalent to the crs2
in our last test so crs1
no longer equals crs2
. I think I'd consider that a bug in PROJ 9.3. I'll see if I can come up with a C example and file a bug with them, but I might need help from Alan Snow (pyproj maintainer).
@avalentino I think a newer version of PROJ fixes the issue(s) identified here. Can you confirm that builds are passing for you or do you think there is more work to do?
Nevermind, PROJ 9.3.1 isn't released yet.
FYI it looks like PROJ 9.3.1 was released in December. The bug fix to make pyresample tests pass should be included in that.
I confirm that all works now. Please feel free to close the issue.
Code Sample, a minimal, complete, and verifiable piece of code
Problem description
Unittests fail with PROJ 9.3.0. Error detected on Debian sid. More details available in https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=1050774 and https://bugs.debian.org/cgi-bin/bugreport.cgi?att=1;bug=1050774;filename=pyresample_1.27.1-1.1_amd64.build;msg=5
Expected Output
All tests pass successfully.
Actual Result, Traceback if applicable
Versions of Python, package at hand and relevant dependencies