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.

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

satpy_cf_nc reader fails to read area written by cf writer for projections not supported by CF conventions #2226

Open
gerritholl opened this issue Oct 7, 2022 · 19 comments

Comments

@gerritholl
Copy link
Collaborator

Describe the bug

I'm writing a data variable to a NetCDF file using the cf writer, then reading it with the satpy_cf_nc reader. Upon reading, no area information is available.

To Reproduce

from satpy import Scene
from glob import glob
from satpy.utils import debug_on; debug_on()
filenames = glob("/media/nas/x21308/scratch/wcm/in/*C13_G16_*")
sc = Scene(filenames=filenames, reader=["abi_l1b"])
sc.load(["C13"])
ls = sc.resample("wcm40km")
ls.save_datasets(
        writer="cf",
        filename="/tmp/{platform_name}-{sensor}-{start_time:%Y%m%d%H%M%S}-{end_time:%Y%m%d%H%M%S}.nc",
        encoding={
            "C13":
            {"zlib": True,
             "complevel": 9,
             "scale_factor": 0.1,
             "dtype": "int16",
             "_FillValue": 0}},
        include_lonlats=False)
sc2 = Scene(filenames=["/tmp/GOES-16-abi-20221006120020-20221006120952.nc"], reader=["satpy_cf_nc"])
sc2.load(["C13"])
print(sc2["C13"].attrs["area"])

Expected behavior

I expect all area/geolocation information to be retained.

Actual results

The writtes has written a NetCDF file that contains a data variable wcm40km with an attribute crs_kwt that contains area metadata, but the satpy_cf_nc reader does not appear to find this. Full output from the MCVE:

/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 15:42:00 : satpy.readers.yaml_reader] Reading ('/data/gholl/checkouts/satpy/satpy/etc/readers/abi_l1b.yaml',)
[DEBUG: 2022-10-07 15:42:00 : satpy.readers.yaml_reader] Assigning to abi_l1b: ['/media/nas/x21308/scratch/wcm/in/OR_ABI-L1b-RadF-M6C13_G16_s20222791200206_e20222791209526_c20222791209590.nc']
[DEBUG: 2022-10-07 15:42:00 : satpy.composites.config_loader] Looking for composites config file abi.yaml
/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 15:42:01 : satpy.composites.config_loader] Looking for composites config file visir.yaml
[DEBUG: 2022-10-07 15:42:01 : satpy.readers.abi_l1b] Reading in get_dataset C13.
[DEBUG: 2022-10-07 15:42:01 : satpy.scene] Resampling DataID(name='C13', wavelength=WavelengthRange(min=10.1, central=10.35, max=10.6, unit='µm'), resolution=2000, calibration=<calibration.brightness_temperature>, modifiers=())
[INFO: 2022-10-07 15:42:01 : satpy.resample] Using default KDTree resampler
[DEBUG: 2022-10-07 15:42:01 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2022-10-07 15:42:01 : satpy.resample] Resampling Rad
[DEBUG: 2022-10-07 15:42:01 : satpy.writers] Reading ['/data/gholl/checkouts/satpy/satpy/etc/writers/cf.yaml']
[INFO: 2022-10-07 15:42:01 : satpy.writers.cf_writer] Saving datasets to NetCDF4/CF.
/data/gholl/checkouts/satpy/satpy/writers/cf_writer.py:571: FutureWarning: The default behaviour of the CF writer will soon change to not compress data by default.
  warnings.warn("The default behaviour of the CF writer will soon change to not compress data by default.",
/data/gholl/checkouts/satpy/satpy/writers/cf_writer.py:197: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  if LooseVersion(pyproj.__version__) < LooseVersion('2.4.1'):
[WARNING: 2022-10-07 15:42:01 : satpy.writers.cf_writer] No time dimension in datasets, skipping time bounds creation.
/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in cos
  return func(*(_execute_task(a, cache) for a in args))
/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in sin
  return func(*(_execute_task(a, cache) for a in args))
/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 15:42:08 : satpy.readers.yaml_reader] Reading ('/data/gholl/checkouts/satpy/satpy/etc/readers/satpy_cf_nc.yaml',)
[DEBUG: 2022-10-07 15:42:08 : satpy.readers.yaml_reader] Assigning to satpy_cf_nc: ['/tmp/GOES-16-abi-20221006120020-20221006120952.nc']
/data/gholl/checkouts/satpy/satpy/readers/satpy_cf_nc.py:241: DeprecationWarning: The truth value of an empty array is ambiguous. Returning False, but in future this will result in an error. Use `array.size > 0` to check that an array is not empty.
  if 'modifiers' in ds_info and not ds_info['modifiers']:
[DEBUG: 2022-10-07 15:42:08 : satpy.readers.satpy_cf_nc] Getting data for: C13
[DEBUG: 2022-10-07 15:42:08 : satpy.readers.satpy_cf_nc] No AreaDefinition to load from nc file. Falling back to SwathDefinition.
[DEBUG: 2022-10-07 15:42:08 : satpy.readers.yaml_reader] No coordinates found for DataID(name='C13', wavelength=WavelengthRange(min=10.1, central=10.35, max=10.6, unit='µm'), resolution=2000, calibration=<calibration.brightness_temperature>, modifiers=())
Traceback (most recent call last):
  File "/data/gholl/checkouts/protocode/mwe/test-nc-area.py", line 21, in <module>
    print(sc2["C13"].attrs["area"])
KeyError: 'area'

Headers of the NetCDF file:

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" ;
        short C13(y, x) ;
                C13:_FillValue = 0s ;
                C13:calibration = "brightness_temperature" ;
                C13:cell_methods = "t: point area: point" ;
                C13:end_time = "2022-10-06 12:09:52.600000" ;
                C13:grid_mapping = "wcm40km" ;
                C13:instrument_ID = "FM1" ;
                C13:long_name = "Brightness Temperature" ;
                C13:modifiers = "" ;
                C13:observation_type = "Rad" ;
                C13:orbital_parameters = "{\"projection_longitude\": -75.0, \"projection_latitude\": 0.0, \"projection_altitude\": 35786023.0, \"satellite_nominal_latitude\": 0.0, \"satellite_nominal_longitude\": -75.19999694824219, \"satellite_nominal_altitude\": 35786023.4375, \"yaw_flip\": false}" ;
                C13:orbital_slot = "GOES-East" ;
                C13:platform_name = "GOES-16" ;
                C13:platform_shortname = "G16" ;
                C13:production_site = "RBU" ;
                C13:reader = "abi_l1b" ;
                C13:resolution = 2000LL ;
                C13:scan_mode = "M6" ;
                C13:scene_abbr = "F" ;
                C13:scene_id = "Full Disk" ;
                C13:sensor = "abi" ;
                C13:sensor_band_bit_depth = 12b ;
                C13:standard_name = "toa_brightness_temperature" ;
                C13:start_time = "2022-10-06 12:00:20.600000" ;
                C13:units = "K" ;
                string C13:wavelength = "10.35 µm (10.1-10.6 µm)" ;
                C13:scale_factor = 0.1 ;

// global attributes:
                :history = "Created by pytroll/satpy on 2022-10-07 13:42:01.922409" ;
                :Conventions = "CF-1.7" ;
}

Environment Info:

  • OS: openSUSE Leap 15.3
  • Satpy Version: v0.37.1-116-g7a30eb87
@gerritholl
Copy link
Collaborator Author

If I try eurol instead of wcm40km, the bug does not occur and the area is available.

@gerritholl
Copy link
Collaborator Author

If I try the built-in area worldeqc30km instead of our local wcm40km, the bug still occurs and the area is unavailable.

@gerritholl
Copy link
Collaborator Author

Temporarily replacing logger.debug with logger.exception in SatpyCFFileHandler.get_area_def reveals a ValueError is raised in pyresample.utils.cf.load_cf_area (which is then caught in get_area_def):

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.

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

That grid mapping variable looks like it is missing a lot of parameters it is supposed to have.

@gerritholl
Copy link
Collaborator Author

Right, when I try to use eurol the grid mapping variable has much more metadata:

        int64 eurol ;
                eurol:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"Unknown based on WGS84 ellipsoid\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7030]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Polar Stereographic (variant B)\",ID[\"EPSG\",9829]],PARAMETER[\"Latitude of standard parallel\",60,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8832]],PARAMETER[\"Longitude of origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8833]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",south,MERIDIAN[90,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",south,MERIDIAN[180,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ;
                eurol:false_easting = 0. ;
                eurol:false_northing = 0. ;
                eurol:geographic_crs_name = "unknown" ;
                eurol:grid_mapping_name = "polar_stereographic" ;
                eurol:horizontal_datum_name = "Unknown based on WGS84 ellipsoid" ;
                eurol:inverse_flattening = 298.257223563 ;
                eurol:long_name = "eurol" ;
                eurol:longitude_of_prime_meridian = 0. ;
                eurol:prime_meridian_name = "Greenwich" ;
                eurol:projected_crs_name = "unknown" ;
                eurol:reference_ellipsoid_name = "WGS 84" ;
                eurol:semi_major_axis = 6378137. ;
                eurol:semi_minor_axis = 6356752.31424518 ;
                eurol:standard_parallel = 60. ;
                eurol:straight_vertical_longitude_from_pole = 0. ;

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

Before saving to netcdf, what do ls['C13'].attrs['area'].crs.to_cf() give you?

@gerritholl
Copy link
Collaborator Author

gerritholl commented Oct 7, 2022

Before saving to netcdf, what do ls['C13'].attrs['area'].crs.to_cf() give you?

For worldeqc30km, it gives a dictionary with only the crs_wkt key:

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'}

For eurol, it gives a dictionary with many keys:

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]', 'semi_major_axis': 6378137.0, 'semi_minor_axis': 6356752.314245179, 'inverse_flattening': 298.257223563, 'reference_ellipsoid_name': 'WGS 84', 'longitude_of_prime_meridian': 0.0, 'prime_meridian_name': 'Greenwich', 'geographic_crs_name': 'unknown', 'horizontal_datum_name': 'Unknown based on WGS84 ellipsoid', 'projected_crs_name': 'unknown', 'grid_mapping_name': 'polar_stereographic', 'standard_parallel': 60.0, 'straight_vertical_longitude_from_pole': 0.0, 'false_easting': 0.0, 'false_northing': 0.0}

@gerritholl
Copy link
Collaborator Author

Is all this extra information actually necessary, though? Isn't that redundant compared to a WKT that should fully describe the area?

In [4]: pyproj.CRS(4087).to_cf()
Out[4]: {'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]]'}

also contains only one key in the grid mapping dict.

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

Interesting. That definitely seems like a bug in pyproj or maybe something explicitly not supported, but the proj dict for that worldeqc30km is nothing special.

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

@snowman2 I can't seem to get pyproj CRS to give me a full grid mapping for the eqc proj...Oh! is it not supported by CF?

In [1]: from pyproj import CRS

In [2]: crs = CRS.from_dict({"proj": "eqc", "ellps": "WGS84"})

In [3]: crs.to_cf()
Out[3]: {'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'}

In [4]: CRS.from_dict({"proj": "eqc"}).to_cf()
Out[4]: {'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'}

In [5]: CRS.from_dict({"proj": "eqc", "lon_0": -45.0, "datum": "WGS84"}).to_cf()
Out[5]: {'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",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",-45,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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'}

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

Yep, that's it @gerritholl.

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/apf.html

https://github.com/pyproj4/pyproj/blob/main/pyproj/crs/_cf1x8.py

So the question is: Should the CF->AreaDefinition conversion in pyresample accept this and fallback to use the crs_wkt? Or should the CF writer in Satpy fail if it detects that the CRS is not CF (1.7?) compliant?

@djhoese
Copy link
Member

djhoese commented Oct 7, 2022

Is all this extra information actually necessary, though? Isn't that redundant compared to a WKT that should fully describe the area?

@gerritholl the question is what is CF compliant and what do we support? I think @snowman2 has argued with the CF folks for a while now that WKT is how they should be describing their CRSes and not having to manually add new projections with new parameter/attribute names for each new revision of the CF standard (we don't need yet another way of describing CRSes, PROJ has enough as it is). I think the pyresample library could be updated to accept/understand this case and return an AreaDefinition with a CRS purely from the crs_wkt information.

@gerritholl
Copy link
Collaborator Author

From https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#appendix-grid-mappings:

The crs_wkt attribute is intended to act as a supplement to other single-property CF grid mapping attributes (as described in Appendix F); it is not intended to replace those attributes. If data producers omit the single-property grid mapping attributes in favour of the crs_wkt attribute, software which cannot interpret crs_wkt will be unable to use the grid_mapping information. Therefore the CRS should be described as thoroughly as possible with the single-property grid mapping attributes as well as by crs_wkt.

It would seem that the grid mapping returned by pyproj is not CF-compliant for CRS 4087, which is why the round-trip fails between the cf writer and the satpy_cf_nc reader. This could be resolved by making the satpy_cf_nc reader able to work with the crs_wkt attribute alone, assuming the information contained therein is a complete description.

For other CRSes, such as CRS 4326, pyproj does return a grid mapping including other keys.

@snowman2
Copy link

snowman2 commented Oct 7, 2022

It would seem that the grid mapping returned by pyproj is not CF-compliant for CRS 4087

It does a best effort. If the CF grid mapping does not exist, the crs_wkt is the only attribute that can be provided. It won't work in some of the CF software, but it is compatible with most geospatial software. See: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html

@gerritholl
Copy link
Collaborator Author

gerritholl commented Oct 7, 2022

The incomplete CF grid mapping is not the full story. Even if the CF grid mapping has more attributes, satpy may fail to read it. Example with EPSG 4326:

$ cat read-nc-area.py
from satpy import Scene
from satpy.utils import debug_on; debug_on()
sc = Scene(filenames=["/media/nas/x21308/scratch/wcm/tmp/Meteosat-9-seviri-20221006120009-20221006121241.nc"], reader=["satpy_cf_nc"])
sc.load(["IR_108"])
print(sc["IR_108"].attrs["area"])
$ python test-nc-area.py
/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 17:11:35 : satpy.readers.yaml_reader] Reading ('/data/gholl/checkouts/satpy/satpy/etc/readers/abi_l1b.yaml',)
[DEBUG: 2022-10-07 17:11:35 : satpy.readers.yaml_reader] Assigning to abi_l1b: ['/media/nas/x21308/scratch/wcm/in/OR_ABI-L1b-RadF-M6C13_G16_s20222791200206_e20222791209526_c20222791209590.nc']
[DEBUG: 2022-10-07 17:11:35 : satpy.composites.config_loader] Looking for composites config file abi.yaml
/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 17:11:35 : satpy.composites.config_loader] Looking for composites config file visir.yaml
[DEBUG: 2022-10-07 17:11:35 : satpy.readers.abi_l1b] Reading in get_dataset C13.
[DEBUG: 2022-10-07 17:11:35 : satpy.scene] Resampling DataID(name='C13', wavelength=WavelengthRange(min=10.1, central=10.35, max=10.6, unit='µm'), resolution=2000, calibration=<calibration.brightness_temperature>, modifiers=())
[INFO: 2022-10-07 17:11:36 : satpy.resample] Using default KDTree resampler
[DEBUG: 2022-10-07 17:11:36 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2022-10-07 17:11:36 : satpy.resample] Resampling Rad
[DEBUG: 2022-10-07 17:11:36 : satpy.writers] Reading ['/data/gholl/checkouts/satpy/satpy/etc/writers/cf.yaml']
[INFO: 2022-10-07 17:11:36 : satpy.writers.cf_writer] Saving datasets to NetCDF4/CF.
/data/gholl/checkouts/satpy/satpy/writers/cf_writer.py:571: FutureWarning: The default behaviour of the CF writer will soon change to not compress data by default.
  warnings.warn("The default behaviour of the CF writer will soon change to not compress data by default.",
/data/gholl/checkouts/satpy/satpy/writers/cf_writer.py:197: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  if LooseVersion(pyproj.__version__) < LooseVersion('2.4.1'):
[WARNING: 2022-10-07 17:11:36 : satpy.writers.cf_writer] No time dimension in datasets, skipping time bounds creation.
/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in cos
  return func(*(_execute_task(a, cache) for a in args))
/data/gholl/mambaforge/envs/py310/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in sin
  return func(*(_execute_task(a, cache) for a in args))
{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'}
/data/gholl/checkouts/satpy/satpy/_config.py:125: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in entry_points().get(name, []):
[DEBUG: 2022-10-07 17:11:42 : satpy.readers.yaml_reader] Reading ('/data/gholl/checkouts/satpy/satpy/etc/readers/satpy_cf_nc.yaml',)
[DEBUG: 2022-10-07 17:11:42 : satpy.readers.yaml_reader] Assigning to satpy_cf_nc: ['/tmp/GOES-16-abi-20221006120020-20221006120952.nc']
/data/gholl/checkouts/satpy/satpy/readers/satpy_cf_nc.py:241: DeprecationWarning: The truth value of an empty array is ambiguous. Returning False, but in future this will result in an error. Use `array.size > 0` to check that an array is not empty.
  if 'modifiers' in ds_info and not ds_info['modifiers']:
[DEBUG: 2022-10-07 17:11:42 : satpy.readers.satpy_cf_nc] Getting data for: C13
[ERROR: 2022-10-07 17:11:42 : satpy.readers.satpy_cf_nc] No AreaDefinition to load from nc file. Falling back to SwathDefinition.
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.
[DEBUG: 2022-10-07 17:11:42 : satpy.readers.yaml_reader] No coordinates found for DataID(name='C13', wavelength=WavelengthRange(min=10.1, central=10.35, max=10.6, unit='µm'), resolution=2000, calibration=<calibration.brightness_temperature>, modifiers=())
Traceback (most recent call last):
  File "/data/gholl/checkouts/protocode/mwe/test-nc-area.py", line 22, in <module>
    print(sc2["C13"].attrs["area"])
KeyError: 'area'
$ ncdump -h /media/nas/x21308/scratch/wcm/tmp/Meteosat-9-seviri-20221006120009-20221006121241.nc
netcdf Meteosat-9-seviri-20221006120009-20221006121241 {
dimensions:
        y = 500 ;
        x = 1000 ;
variables:
        int64 world4326 ;
                world4326:crs_wkt = "GEOGCRS[\"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]],CS[ellipsoidal,2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"Horizontal component of 3D system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]" ;
                world4326:geographic_crs_name = "WGS 84" ;
                world4326:grid_mapping_name = "latitude_longitude" ;
                world4326:inverse_flattening = 298.257223563 ;
                world4326:long_name = "world4326" ;
                world4326:longitude_of_prime_meridian = 0. ;
                world4326:prime_meridian_name = "Greenwich" ;
                world4326:reference_ellipsoid_name = "WGS 84" ;
                world4326:semi_major_axis = 6378137. ;
                world4326:semi_minor_axis = 6356752.31424518 ;
        double y(y) ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "m" ;
        double x(x) ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "m" ;
        short IR_108(y, x) ;
                IR_108:_FillValue = 0s ;
                IR_108:calibration = "brightness_temperature" ;
                IR_108:end_time = "2022-10-06 12:12:41.110000" ;
                IR_108:georef_offset_corrected = "true" ;
                IR_108:grid_mapping = "world4326" ;
                IR_108:modifiers = "" ;
                IR_108:nominal_end_time = "2022-10-06 12:15:09.596170" ;
                IR_108:nominal_start_time = "2022-10-06 12:00:09.961253" ;
                IR_108:orbital_parameters = "{\"projection_longitude\": 45.5, \"projection_latitude\": 0.0, \"projection_altitude\": 35785831.0, \"satellite_nominal_longitude\": 45.5, \"satellite_nominal_latitude\": 0.0, \"satellite_actual_longitude\": 45.7482511326184, \"satellite_actual_latitude\": 0.5197476531151144, \"satellite_actual_altitude\": 35777087.07270167}" ;
                IR_108:platform_name = "Meteosat-9" ;
                IR_108:reader = "seviri_l1b_hrit" ;
                IR_108:resolution = 3000.403165817 ;
                IR_108:sensor = "seviri" ;
                IR_108:standard_name = "toa_brightness_temperature" ;
                IR_108:start_time = "2022-10-06 12:00:09.961000" ;
                IR_108:units = "K" ;
                string IR_108:wavelength = "10.8 µm (9.8-11.8 µm)" ;
                IR_108:scale_factor = 0.1 ;

// global attributes:
                :history = "Created by pytroll/satpy on 2022-10-07 15:07:22.749624" ;
                :Conventions = "CF-1.7" ;
}

@gerritholl
Copy link
Collaborator Author

For EPSG 4326 the problem is different. The standard_name on y is projection_y_coordinate in m, which is wrong for EPSG 4326, as y should be longitude with units degrees. I'll open a separate issue for that.

@gerritholl
Copy link
Collaborator Author

Digging into this some more, it seems that EPSG 4087 (equirectangular cylindrical projection) isn't supported by the CF Conventions, and therefore, writing a CF conform file with this projection is not possible. https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#appendix-grid-mappings

@gerritholl
Copy link
Collaborator Author

Should the cf writer raise a warning or an exception when the user attempts to write data that cannot be written in a CF-conform way?

@djhoese
Copy link
Member

djhoese commented Oct 13, 2022

I'd say warning if we don't want any new keyword arguments. We shouldn't stop anyone from producing a file unless they have a keyword argument to allow a less-strict CF check.

There is a strict argument to the area2cf helper function which controls whether or not lon/lats should be written. Maybe that could be used for the same idea (the writer warns on strict=False and "bad" CRS, otherwise exception if strict=True and "bad" CRS).

https://satpy.readthedocs.io/en/stable/_modules/satpy/writers/cf_writer.html#area2cf

@gerritholl gerritholl changed the title satpy_cf_nc reader fails to read area written by cf writer satpy_cf_nc reader fails to read area written by cf writer for projections not supported by CF conventions Jan 10, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants