search

pytroll / satpy

Python package for earth-observing satellite data processing
http://satpy.readthedocs.org/en/latest/
GNU General Public License v3.0
1.07k stars 295 forks source link

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

Open gerritholl opened 2 years ago

gerritholl commented 2 years ago

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:

gerritholl commented 2 years ago

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

gerritholl commented 2 years ago

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

gerritholl commented 2 years ago

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 commented 2 years ago

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

gerritholl commented 2 years ago

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 commented 2 years ago

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

gerritholl commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

@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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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