search

SciTools / iris

A powerful, format-agnostic, and community-driven Python package for analysing and visualising Earth science data
https://scitools-iris.readthedocs.io/en/stable/
BSD 3-Clause "New" or "Revised" License
626 stars 283 forks source link

Error loading NSDIC-0116-nh data. CRS AXIS UNKNOWN (polar-grid) #5343

Closed rbeucher closed 11 months ago

rbeucher commented 1 year ago

🐛 Bug Report

I am trying to load the NSIDC-0116-nh raw data into a cube. Iris fails with the following error: This is a CRS error from pyproj (v3.5.0).

CRSError: Invalid projection: PROJCS["NSIDC EASE-Grid North",GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere",DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",SPHEROID["International 1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],AUTHORITY["EPSG","6053"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4053"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3408"],AXIS["X",UNKNOWN],AXIS["Y",UNKNOWN]]: (Internal Proj Error: proj_create: unhandled axis direction: UNKNOWN)

Full Traceback: traceback.txt

How To Reproduce

Steps to reproduce the behaviour:

  1. Download data from https://nsidc.org/data/nsidc-0116/versions/4
import iris 
cubes = iris.load("icemotion_daily_nh_25km_19790101_19791231_v4.1.nc")

Expected behaviour

Should load data.

Screenshots

Environment

Additional context

This is used by the NSDIC cmoriser in ESMValTool

valeriupredoi commented 1 year ago

same for sh, for this:

Basically it's like this:

import iris

c = iris.load("/home/valeriu/Desktop/nsdic/Tier3/NSIDC-0116-sh/icemotion_daily_sh_25km_19840101_19841231_v4.1.nc")

spitting this out:

Traceback (most recent call last):
  File "/home/valeriu/iris_load_nsdic.py", line 3, in <module>
    c = iris.load("/home/valeriu/Desktop/nsdic/Tier3/NSIDC-0116-sh/icemotion_daily_sh_25km_19840101_19841231_v4.1.nc")
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/__init__.py", line 319, in load
    return _load_collection(uris, constraints, callback).merged().cubes()
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/__init__.py", line 285, in _load_collection
    result = _CubeFilterCollection.from_cubes(cubes, constraints)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/cube.py", line 105, in from_cubes
    for cube in cubes:
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/__init__.py", line 270, in _generate_cubes
    for cube in iris.io.load_files(part_names, callback, constraints):
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/io/__init__.py", line 223, in load_files
    for cube in handling_format_spec.handler(
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/netcdf/loader.py", line 566, in load_cubes
    cube = _load_cube(engine, cf, cf_var, filename)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/netcdf/loader.py", line 276, in _load_cube
    engine.activate()
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/engine.py", line 97, in activate
    run_actions(self)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/actions.py", line 528, in run_actions
    action_provides_grid_mapping(engine, grid_mapping_fact)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/actions.py", line 74, in inner
    rule_name = func(engine, *args, **kwargs)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/actions.py", line 172, in action_provides_grid_mapping
    coordinate_system = builder(engine, cf_var)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/helpers.py", line 719, in build_lambert_azimuthal_equal_area_coordinate_system
    ellipsoid = _get_ellipsoid(cf_grid_var)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/iris/fileformats/_nc_load_rules/helpers.py", line 467, in _get_ellipsoid
    proj_crs = pyproj.crs.CRS.from_wkt(crs_wkt)
  File "/home/valeriu/miniconda3/envs/tool-latest2/lib/python3.10/site-packages/pyproj/crs/crs.py", line 431, in from_wkt
    raise CRSError(f"Invalid WKT string: {in_wkt_string}")
pyproj.exceptions.CRSError: Invalid WKT string:  PROJCS["NSIDC EASE-Grid South",GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere",DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",SPHEROID["International 1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],AUTHORITY["EPSG","6053"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4053"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",-90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3409"],AXIS["X",UNKNOWN],AXIS["Y",UNKNOWN]]

I can provide the file, it's only 22MB in heft. Cheers guys :beer:

remi-kazeroni commented 1 year ago

I could also reproduce this with iris 3.4.1 and 3.6.0, used in recent ESMValTool installations. The dataset could be loaded successfully with iris 3.2.1 that we were using for ESMValTool v2.7

rbeucher commented 1 year ago

[like] Romain Beucher reacted to your message:

From: Rémi Kazeroni @.> Sent: Wednesday, June 14, 2023 10:20:20 AM To: SciTools/iris @.> Cc: Romain Beucher @.>; Author @.> Subject: Re: [SciTools/iris] Error loading NSDIC-0116-nh data. CRS AXIS UNKNOWN (polar-grid) (Issue #5343)

I could also reproduce this with iris 3.4.1 and 3.6.0, used in recent ESMValTool installations. The dataset could be loaded successfully with iris 3.2.1 that we were using for ESMValTool v2.7

— Reply to this email directly, view it on GitHubhttps://github.com/SciTools/iris/issues/5343#issuecomment-1590918806, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AA4CVFUTG35OJTLOYZOLONDXLGF6JANCNFSM6AAAAAAZCTUDOE. You are receiving this because you authored the thread.Message ID: @.***>

rbeucher commented 1 year ago

Looks like a pyproj issue to me...

trexfeathers commented 1 year ago

@rbeucher just to let you know we have some scheduled time to investigate this later in the month.

rbeucher commented 1 year ago

Thanks for the update.

trexfeathers commented 1 year ago

Since it's unclear whether this problem lies in Proj/Cartopy/Iris, having the example data is important for investigation. We've so far had trouble getting hold of the data (have tried several contacts), so this is unlikely to be investigated before the release of Iris 3.7.

Do get in touch if you have a copy of a relevant file!

valeriupredoi commented 1 year ago

Hi @trexfeathers I just emailed your AVD colleague Elias a copy of the file that created troubles for me, but am also gonna attach a txt-ed extension-ed file here (fool GH to ingest it as attachement) icemotion_daily_sh_25km_19840101_19841231_v4.1.nc.txt

valeriupredoi commented 1 year ago

you guys should get your :duck: in line on the getting data business side though :grin:

valeriupredoi commented 1 year ago

OK the email bounced back from the MO email client since it complains it's too large (21MB is not a large message :laughing: ) so this is the way to go, attached above :arrow_up:

rbeucher commented 1 year ago

Let me know if I can help. This is annoying.

trexfeathers commented 1 year ago

Thanks for sending the file over. This enabled us to step through and experiment; the problem is because Proj does not accept leading spaces in WKT strings, so the crs_wkt attribute on the crs variable:

 PROJCS["NSIDC EASE-Grid South",GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere",DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",SPHEROID["International 1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],AUTHORITY["EPSG","6053"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4053"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",-90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3409"],AXIS["X",UNKNOWN],AXIS["Y",UNKNOWN]]

... is invalid. If the leading space is removed, the file loads fine.

Iris has been handling WKT for over a year (#4704) but this is the first time the problem has been reported, which suggests that either:

I'd like to draw on ESMValTool's rich experience with file formats here: would you expect this problem to crop up more often in future, or is there a specific problem with NSIDC-0116-nh files?

valeriupredoi commented 1 year ago

Glad it helped! -sh files too - we've not seen anything like this before so am gonna coin this as a very specific issue related to the way NSIDC writes their data; @remi-kazeroni have you seen this type of behaviour before? What do you propose as solution @trexfeathers ? I'd argue that this could be a nice issue at Proj upstream - if it's just a leading space that's causing it, I reckon they'd prob best take care of this? :beer:

trexfeathers commented 1 year ago

What do you propose as solution @trexfeathers ?

Well I'd always prefer to push the solution to those that write the files. The trouble with putting the fix into Proj/Iris/ESMValTool is that our users might still have trouble passing the file to collaborators using different tools.

But I know files cannot always be fixed. If that is the situation here:

trexfeathers commented 1 year ago

Actually it appears that Proj explicitly copes with leading spaces (OSGeo/PROJ@93d8fccd51a704c4accd535743a963ad0ca5a274), suggesting the problem might lie in PyProj.

I'm definitely submitting an issue to PyProj if that means it isn't aligning with Proj.

valeriupredoi commented 1 year ago

from our experience this data will probably not be regenerated for this version 4.1 (dates from April 2019, produced on me birthday, of all days) - they'll shrug it off, especially for such a minor formatting issue. Note that we have a CMORizer script exactly for this data which is probably broken now since of this issue, so things have worked in the past (prob before Iris started supporting WKTs am guessing)

valeriupredoi commented 1 year ago

Actually it appears that Proj explicitly copes with leading spaces (OSGeo/PROJ@93d8fcc), suggesting the problem might lie in PyProj.

I'm definitely submitting an issue to PyProj if that means it isn't aligning with Proj.

ah bingo! Godspeed, Martin :checkered_flag:

trexfeathers commented 1 year ago

pyproj4/pyproj#1328

valeriupredoi commented 1 year ago

that is one beautifully written/formatted issue @trexfeathers :100:

trexfeathers commented 1 year ago

https://github.com/pyproj4/pyproj/issues/1328#issuecomment-1669577334:

This issue will be addressed in the next PROJ release

rbeucher commented 12 months ago

PROJ and pyproj have new releases. I will give it another go

trexfeathers commented 11 months ago

@rbeucher any success?

bjlittle commented 11 months ago

@SciTools/peloton We're going to close this issue @rbeucher as we believe the problem has been addressed in proj.

If you manage to get around to testing and still consider this an issue, then please don't hesitate to re-open this issue.

rbeucher commented 11 months ago

Yes. Sorry for the delay. I believe the pb is fixed. Thank you for your help.