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

[rbeucher]

🐛 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

• Iris Version: 3.6.0

Additional context

This is used by the NSDIC cmoriser in ESMValTool

[valeriupredoi]

same for sh, for this:

• iris 3.4.1 and 3.6.0
• multiple pyproj including latest 3.6.0
• python 3,.10 and 3.11

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 🍺

[remi-kazeroni]

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]

[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 GitHub<#5343 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AA4CVFUTG35OJTLOYZOLONDXLGF6JANCNFSM6AAAAAAZCTUDOE>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[rbeucher]

Looks like a pyproj issue to me...

[trexfeathers]

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

[rbeucher]

Thanks for the update.

[trexfeathers]

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]

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]

you guys should get your 🦆 in line on the getting data business side though 😁

[valeriupredoi]

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

[rbeucher]

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

[trexfeathers]

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:

• Most crs_wkt attributes are correctly formatted for Proj to handle.
• Users don't often use Iris to load files with crs_wkt attributes.

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?