From b4f71af55bc8eb58a543b7d253eb774d315621ab Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 13 Oct 2022 10:05:19 +0200 Subject: [PATCH 1/2] Add test for loading CRS WKT Add a test for loading an area from a CF NetCDF file encoded only with a CRS WKT. The test currently fails as I did not yet implement anything. --- pyresample/test/test_utils.py | 76 +++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py index d28bfdcb9..e8c0c41e1 100644 --- a/pyresample/test/test_utils.py +++ b/pyresample/test/test_utils.py @@ -683,6 +683,82 @@ def validate_llnocrs(adef, cfinfo, lat='lat', lon='lon'): adef, cf_info = load_cf_area(cf_file) validate_llnocrs(adef, cf_info) + def test_load_from_crs_wkt(self): + """Test that we can load with CRS WKT only. + + Test that we can reconstruct the area when the grid mapping variable + only contains the CRS WKT attribute. This is the case for + equirectangular projections, for example, ``pyproj.CRS(4087).to_cf()``. + + See also https://github.com/pytroll/pyresample/issues/460. + """ + import xarray as xr + + from pyresample import create_area_def + from pyresample.utils import load_cf_area + + width = height = 100 + crs_wkt = """ +PROJCRS["WGS 84 / World Equidistant Cylindrical", + BASEGEOGCRS["WGS 84", + ENSEMBLE["World Geodetic System 1984 ensemble", + MEMBER["World Geodetic System 1984 (Transit)"], + MEMBER["World Geodetic System 1984 (G730)"], + MEMBER["World Geodetic System 1984 (G873)"], + MEMBER["World Geodetic System 1984 (G1150)"], + MEMBER["World Geodetic System 1984 (G1674)"], + MEMBER["World Geodetic System 1984 (G1762)"], + MEMBER["World Geodetic System 1984 (G2139)"], + ELLIPSOID["WGS 84",6378137,298.257223563, + LENGTHUNIT["metre",1]], + ENSEMBLEACCURACY[2.0]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433]], + ID["EPSG",4326]], + CONVERSION["World Equidistant Cylindrical", + METHOD["Equidistant Cylindrical", + ID["EPSG",1028]], + PARAMETER["Latitude of 1st standard parallel",0, + ANGLEUNIT["degree",0.0174532925199433], + ID["EPSG",8823]], + PARAMETER["Longitude of natural origin",0, + ANGLEUNIT["degree",0.0174532925199433], + ID["EPSG",8802]], + PARAMETER["False easting",0, + LENGTHUNIT["metre",1], + ID["EPSG",8806]], + PARAMETER["False northing",0, + LENGTHUNIT["metre",1], + ID["EPSG",8807]]], + CS[Cartesian,2], + AXIS["easting (X)",east, + ORDER[1], + LENGTHUNIT["metre",1]], + AXIS["northing (Y)",north, + ORDER[2], + LENGTHUNIT["metre",1]], + USAGE[ + SCOPE["Graticule coordinates expressed in simple Cartesian form."], + AREA["World."], + BBOX[-90,-180,90,180]], + ID["EPSG",4087]] +""" + ds = xr.Dataset( + {"data": ( + ('y', 'x'), + np.zeros((width, height)), + {"grid_mapping": "area"}), + "area": ( + (), + 0, + {"crs_wkt": crs_wkt, "long_name": "test"})}, + coords={'x': np.linspace(-7049500, -6950500, width), + 'y': np.linspace(-49500, 49500, height)}) + (adef, cf_info) = load_cf_area(ds) + assert adef == create_area_def( + "test", 4087, resolution=1000, width=100, + height=100, center=(-7000000, 0)) + class TestLoadCFArea_Private(unittest.TestCase): """Test private routines involved in loading an AreaDefinition from netCDF/CF files.""" From 3cf9b5ce75e30de7886b1010e383ba9eb6601431 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 13 Oct 2022 10:48:13 +0200 Subject: [PATCH 2/2] refactor in cf.py Refactor out _get_type_of_grid_mapping in utils/cf.py in preparation for getting the grid mapping from the WKT / CRS directly. Tolerate absence of grid mapping name if there is a crs_wkt attribute. --- pyresample/utils/cf.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/pyresample/utils/cf.py b/pyresample/utils/cf.py index 23871d4f2..f8778daba 100644 --- a/pyresample/utils/cf.py +++ b/pyresample/utils/cf.py @@ -92,7 +92,9 @@ def _load_crs_from_cf_gridmapping(nc_handle, grid_mapping_varname): raise ValueError("Not a valid CF grid_mapping variable ({})".format(grid_mapping_varname)) except AttributeError: # no :grid_mapping_name thus it cannot be a valid grid_mapping variable - raise ValueError("Not a valid CF grid_mapping variable ({})".format(grid_mapping_varname)) + # but let's try to get everything from the CRS if we can + if "crs_wkt" not in v.attrs: + raise ValueError("Not a valid CF grid_mapping variable ({})".format(grid_mapping_varname)) # use pyproj to load the CRS return pyproj.CRS.from_cf(v.attrs) @@ -334,15 +336,8 @@ def _load_cf_area_one_variable(nc_handle, variable, y=None, x=None): crs, grid_mapping_variable, variable_is_itself_gridmapping = _load_cf_area_one_variable_crs(nc_handle, variable) # the type of grid_mapping (its grid_mapping_name) impacts several aspects of the CF reader - if grid_mapping_variable == 'latlon_default': - type_of_grid_mapping = 'latitude_longitude' - else: - try: - type_of_grid_mapping = nc_handle[grid_mapping_variable].grid_mapping_name - except AttributeError: - raise ValueError( - ("Not a valid CF grid_mapping variable ({}):" - "it lacks a :grid_mapping_name attribute").format(grid_mapping_variable)) + type_of_grid_mapping = _get_type_of_grid_mapping( + nc_handle, grid_mapping_variable) cf_info['grid_mapping_variable'] = grid_mapping_variable cf_info['type_of_grid_mapping'] = type_of_grid_mapping @@ -380,6 +375,18 @@ def _load_cf_area_one_variable(nc_handle, variable, y=None, x=None): return area_def, cf_info +def _get_type_of_grid_mapping(nc_handle, grid_mapping_variable): + """Get the type of grid mapping.""" + if grid_mapping_variable == 'latlon_default': + return 'latitude_longitude' + try: + return nc_handle[grid_mapping_variable].grid_mapping_name + except AttributeError: + raise ValueError( + ("Not a valid CF grid_mapping variable ({}):" + "it lacks a :grid_mapping_name attribute").format(grid_mapping_variable)) + + def _load_cf_area_several_variables(nc_handle): """Load the AreaDefinition corresponding to several netCDF variables/fields.""" def _indices_unique_AreaDefs(adefs):