Export pyresample cartopy projection info from area_def to netCDF file?

[Amy Huff - NOAA Affiliate]

Hi developers,

I am using pyresample to regrid TEMPO convolved L1b radiances granule data to make a true color RGB composite. I create an AreaDefinition and get the projection info for Cartopy (see below).  Ultimately, I use the tempo_area.to_cartopy_crs() projection info to plot the RGB composite data array on a map.

I'd like to export (save) the tempo_area.to_cartopy_crs() projection info, ideally in a netCDF4 (.nc) file via an xarray dataset. Is this possible?  My objective is to have a stand-alone .nc file containing the RGB data array and the projection info, so an end user can open up the .nc file and plot the RGB composite on a map.  I can make the .nc file with the RGB data array, but I can't figure out how to include the projection info, or if it's even possible to export it, despite internet searches.

I am able to print the properties of the projection vars(crs) and the projection string crs.proj4_init  but I am not sure what to do with this info.

Thanks,

Amy

    # Create AreaDefinition (resolution ~3km)

    tempo_area = create_area_def('tempo_area', {'proj': 'eqc', 'lon_0': 0},

                                 area_extent=[west_lon_map, south_lat_map,

                                              east_lon_map, north_lat_map],

                                 units='degrees', resolution=(4.3E-2, 1.8E-2))

    # Get uniform grid projection for Cartopy

    crs = tempo_area.to_cartopy_crs()

[David Hoese]

Hi Amy,

The `tempo_area.crs` property is a pyproj CRS object and has a "to_cf" method that you can use to get a dictionary of the grid mapping attributes that a CF standard NetCDF file should have in its grid mapping variable.

https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.to_cf

The cartopy CRS, in modern versions of cartopy, should be a subclass of pyproj's CRS class if I recall correctly so it may also have a `to_cf` method. There are other libraries like Satpy that have utilities for converting an Xarray Dataset or DataArray into a CF NetCDF file but that might be overkill for your situation it sounds like. So I think your steps would be:

1. Get the CF attributes dict using "to_cf".

2. Create a single scalar variable for your grid mapping and call it anything (ex. "temp_projection"). I think people typically create this as an int32 with a value of 0. It's just one of those weird things in CF and doesn't really matter as the value of the variable is never used.

3. Add all of the CF dict items to the grid mapping variable as attributes.

4. Add a "grid_mapping" attribute to your RGB variable and its value should be set to the same name as your grid mapping variable (ex. "tempo_projection").

Good luck.

Dave

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/pytroll/a1e335f9-6df9-4c9c-8e34-e31913f13a8cn%40googlegroups.com.

[Amy Huff - NOAA Affiliate]

Hi Dave,

Thanks for your help.  I can't quite get it to work, however.  I can do all of the steps you suggested, with no problem.  But I can't successfully pull out the projection info from my new .nc file and use it to define a Cartopy projection.  

Below is a copy of the metadata from my Dataset.  Here is how I added the "tempo_projection" info:

tempo_area, crs = create_tempo_uniform_grid(west_lon_map,
                                            east_lon_map,
                                            south_lat_map,
                                            north_lat_map)

# Get the cartopy projection dictionary
cf_dict = crs.to_cf()

# Create a new variable for the grid mapping
ds['tempo_projection'] = 0
ds['tempo_projection'].attrs = cf_dict

# Link the 'rgb_scan' data variable to the grid mapping variable
ds['rgb_scan'].attrs['grid_mapping'] = 'tempo_projection'

Subsequently, I tried several different ways to pull out the projection info, including:

CRS.from_wkt(ds['tempo_projection'].attrs['crs_wkt'])      (using pyproj)

ccrs.Projection.from_cf(ds['tempo_projection'])      (using cartopy)

ccrs.Projection.from_wkt(ds['tempo_projection'].attrs['crs_wkt'])      (using cartopy)

None of these worked.  The error messages basically amounted to "the projection isn't valid."  So maybe I missed a step somewhere?

I'd appreciate any advice you may have.

Thanks,

Amy

xarray.Dataset

• Dimensions:
  • granule: 9
  • rows: 2167
  • columns: 1628
  • color: 3
• Coordinates: (0)
• Data variables:
  • rgb_scan
    (granule, rows, columns, color)
    float32
    nan nan nan nan ... nan nan nan nan
    long_name :true color red-green-blue composite arraygrid_mapping :tempo_projectiontempo_bounds :(-14471533.803125564, -6679169.447596413, 2003750.8342789244, 6345210.975216594)
  • coords
    ()
    object
    None
  • tempo_projection
    ()
    int64
    0
    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]]]]
• Indexes: (0)
• Attributes: (12)
  

[David Hoese]

Sorry Amy, I had missed that part of your original email.

My guess is that cartopy's CRS object, which is now a subclass of pyproj's CRS class, doesn't handle the `from_*` methods very well. I was able to get a CRS object without an error by first creating a pyproj CRS object then passing it to `cartopy.crs.Projection(pyproj_crs_obj)`. So if you have WKT then `pyproj_crs_obj = pyproj.crs.CRS.from_wkt(wkt)`. I haven't tested this with actual plotting though.

Also note that for a normal CF-compatible projection your `to_cf()` dictionary should have a lot more than just `crs_wkt`.

Hopefully this makes sense. Good luck.

Dave

To view this discussion, visit https://groups.google.com/d/msgid/pytroll/724b2f90-82c0-4e5b-be4b-ac3d3937aff8n%40googlegroups.com.