search

pytroll / pyresample

Geospatial image resampling in Python
http://pyresample.readthedocs.org
GNU Lesser General Public License v3.0
350 stars 94 forks source link

Reconstruct CRS from WKT only #460

Open gerritholl opened 2 years ago

gerritholl commented 2 years ago

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 https://github.com/pytroll/satpy/issues/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 https://github.com/pytroll/satpy/issues/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.