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 #460

Open
gerritholl opened this issue Oct 12, 2022 · 0 comments · May be fixed by #461
Open

Reconstruct CRS from WKT only #460

gerritholl opened this issue Oct 12, 2022 · 0 comments · May be fixed by #461
Labels
enhancement

Comments

@gerritholl
Copy link
Collaborator

Feature Request

Is your feature request related to a problem? Please describe.

For some areas, area.crs.to_wkt() returns a dictionary that contains only the crs_wkt key and nothing else. I don't know if this is CF-conform, but this is the responsibility of pyproj. For some discussion, see the comments on pytroll/satpy#2226.

For example:

In [63]: ar = create_area_def("test", 4087, resolution=1000, width=100, height=100, center=[0, 0])

In [64]: ar.crs.to_cf()
Out[64]: {'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]]'}

Currently, pyresample cannot construct an AreaDefinition using only this CRS and the corresponding coordinates. Therefore, Satpy's satpy_cf_nc reader cannot read the files written by its cf writer. See pytroll/satpy#2226 for the Satpy perspective.

The NetCDF header may look like:

netcdf GOES-16-abi-20221006120020-20221006120952 {
dimensions:
        y = 499 ;
        x = 1000 ;
variables:
        int64 wcm40km ;
                wcm40km: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]]" ;
                wcm40km:long_name = "wcm40km" ;
        double y(y) ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "m" ;
        double x(x) ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "m" ;

where trying to load the area will fail with

Traceback (most recent call last):
  File "/data/gholl/checkouts/satpy/satpy/readers/satpy_cf_nc.py", line 304, in get_area_def
    area = AreaDefinition.from_cf(self.filename)
  File "/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/pyresample/geometry.py", line 1790, in from_cf
    return load_cf_area(cf_file, variable=variable, y=y, x=x)[0]
  File "/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/pyresample/utils/cf.py", line 474, in load_cf_area
    raise ValueError("Found no AreaDefinitions in this netCDF/CF file.")
ValueError: Found no AreaDefinitions in this netCDF/CF file.

Describe the solution you'd like

I would like that, if the grid mapping variable contains only the attribute crs_wkt and no other information, pyresample should do its best to construct an area using the CRS WKT with other available information.

Describe any changes to existing user workflow

There shouldn't be any.

Additional context

My current workaround involves writing all lat/lons and using a SwathDefinition. This is inefficient, in particular since I'm then resampling back to the area I had when I was writing the NetCDF file to begin with.
@gerritholl gerritholl linked a pull request Oct 13, 2022 that will close this issue
Draft
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
enhancement
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Reconstruct CRS from WKT only gerritholl/pyresample
1 participant
@gerritholl