Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Reconstruct CRS from WKT only #461

Draft
wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Reconstruct CRS from WKT only #461

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b4f71af
Add test for loading CRS WKT
gerritholl Oct 13, 2022
3cf9b5c
refactor in cf.py
gerritholl Oct 13, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 76 additions & 0 deletions pyresample/test/test_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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."""
Expand Down
27 changes: 17 additions & 10 deletions pyresample/utils/cf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down