Merge pull request #5 from nsidc/due-cryo-119-netcdf-coords

Due cryo 119 netcdf coords

Author: andypbarrett

_quarto.yml

@@ -28,7 +28,13 @@ website:
       - section: "How do I..."
         href: nsidc-data-cookbook/how-to-guides/overview.qmd
         contents:
-          - section: "reproject data" 
+          - section: "reproject data?"
+          - section: "work with NetCDF files?"
+            contents:
+            - text: "get the bounding box of a netcdf file"
+              href: nsidc-data-cookbook/how-to-guides/netcdf_cf.qmd
+            - text: "get the latitude and longitudes for grid cells"
+              href: nsidc-data-cookbook/how-to-guides/get_latitude_and_longitude.qmd
       - section: "Tutorials"
         contents:
       - section: "Reference Guides"

nsidc-data-cookbook/how-to-guides/get_latitude_and_longitude.qmd

@@ -0,0 +1,47 @@
+---
+title: "How to get latitude and longitude for the grid cells of a NetCDF file"
+author: Andrew P. Barrett
+date: last-modified
+---
+
+## Problem
+My data are in a NetCDF file but there are no latitude and longitude coordinates.  How do I get the latitude and longitude for each grid cell center?
+
+## Solution
+```
+import numpy as np
+
+import rioxarray
+import xarray as xr
+from pyproj import CRS, Transformer
+
+
+# Load dataset using decode_coords='all'
+ds = xr.open_dataset('example_data/NSIDC0081_SEAICE_PS_N25km_20230627_v2.0.nc',
+                     decode_coords='all')
+
+# Instantiate source pyproj CRS from ds.rio.crs object
+source_crs = CRS(ds.rio.crs.to_wkt())
+# Instantiate destination pyproj CRS using EPSG code for WGS84 
+destination_crs = CRS(4326)
+
+# Instantiate pyproj transformer instance using source and destination CRS
+transform = Transformer.from_crs(source_crs, destination_crs)
+
+# Create 2D arrays of x and y coordinates
+x2d, y2d = np.meshgrid(ds.x, ds.y)
+
+# Calculate latitude and longitudes for each grid cell
+lat, lon = transform.transform(x2d, y2d)
+```
+
+## Discussion
+The workflow above assumes that the NetCDF file is CF-compliant and has a `grid_mapping` variable and projected coordinates defined.  See the `rioxarray` [CRS management documentation](https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#) if this is not the case.
+
+If the dataset `rio.crs` is not set, the `source_crs` could be defined using the EPSG code or other projection information directly, without having to set the dataset crs.
+
+For all datasets on projected grids, latitude and longitudes for grid cells will be 2-dimensional arrays.  x- and y-coordinates for these datasets will be vectors.  So 2-dimensional arrays of x and y coordinates need to be created.  `np.meshgrid` is designed to do this.
+
+If you want the latitudes and longitudes of the corners of grid cells, x and y coordinates can be incremented by half the grid cells width and height, and the resoluting projected coordinates passed to `transform.transformer`.
+
+It is worth considering if you really need to have latitude and longitude arrays for your dataset.  One reason that CF-conventions were changed, was that 2D arrays of geographic coordinates can increase the file size, especially for high resolution datasets.  Projected coordinates as vectors take up much less space.  Plotting tools such as `cartopy` can make maps without latitude and longitude.  If you need to extract cell values for specific geographic coordinates, for example for a weather station or buoy, it is quicker and easier to transform these coordinates into the projected coordinate system of the dataset.  The `xarray.Dataset.sel` or `xarray.DataArray.sel` methods can then be used to select the data for that location. 

nsidc-data-cookbook/how-to-guides/netcdf_cf.qmd

@@ -0,0 +1,61 @@
+---
+title: "How do I get the bounding box of a NetCDF file in latitude and longitude?"
+author: Andrew P. Barrett
+date: last-modified
+---
+
+## Problem
+My data are in a NetCDF file but there are no latitude and longitude coordinates.  How do I find the bounding box in latitude and longitude?
+
+## Solution
+### Using xarray and rioxarray
+```
+import rioxarray
+import xarray as xr
+
+
+# Load dataset using decode_coords='all'
+ds = xr.open_dataset('example_data/NSIDC0081_SEAICE_PS_N25km_20230627_v2.0.nc',
+                     decode_coords='all')
+
+# Get the bounds of the data grid in WGS84 (EPSG code 4362)
+ds.rio.transform_bounds(4326)
+```
+```
+(-180.0, 30.98056405144958, 180.0, 90.0)
+```
+
+## Discussion
+
+As of [CF v1.8](https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#grid-mappings-and-projections), latitude and longitude coordinate variables are no longer required in CF-compliant NetCDF files, as long as projected horizontal spatial coordinates (e.g. `x` and `y`, or `easting` and `northing`) and a `grid_mapping` variable is provided.  The `grid_mapping` variable defines the coordinate reference system (CRS) of the projected horizontal coordinates.
+
+The easiest, and frankly best, way to read and work with NetCDF files is to use `xarray`.  `rioxarray` extends the `xarray` package to use CRS and make geospatial tasks, such as reprojecting and regridding, easier.  By setting the keyword `decode-coords='all'`, `rioxarray` searches the `xarray.DataArray` or `xarray.Dataset` for the CRS.  `rioxarray` also determines the `transform` or `geotransform`, which defines the image CRS that transforms cell coordinates (column and row) into projected coordinates (x, y).  The `transform` is calculated from the coordinates of the data.  See the `rioxarray` [documentaton](https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#) for more details.
+
+Both these spatial coordinates can be accessed as follows:
+
+```
+ds.rio.crs
+```
+```
+CRS.from_wkt('PROJCS["NSIDC Sea Ice Polar Stereographic North",GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",SPHEROID["Hughes 1980",6378273,298.279411123061,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","6054"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4054"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3411"]]')
+```
+
+```
+ds.rio.transform()
+```
+```
+Affine(25000.0, 0.0, -3850000.0,
+       0.0, -25000.0, 5850000.0)
+```
+
+If this information is not found, all is not lost.  [`rio.write_crs`](https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Setting-the-CRS) and [`rio.write_transform`](https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Setting-the-transform-of-the-dataset) can be used to set the CRS and the transform of the dataset.
+
+Once these variables are set, the horizontal spatial bounds of the dataset can be queried as described above.  If all you want is the bounds in projected coordinates, these can be accessed using `ds.rio.bounds()`.
+
+```
+ds.rio.bounds()
+```
+```
+(-3850000.0, -5350000.0, 3750000.0, 5850000.0)
+```
+

nsidc-data-cookbook/how-to-guides/overview.qmd

@@ -25,4 +25,6 @@ _Maybe add a short description of each section here_
 - Searching for data; 
 - Accessing data;
 - Reprojecting and resampling data;
+- Working with CF-compliant NetCDF files
+  - [How do I get the bounding box of a NetCDF file in latitude and longitude?](netcdf_cf.qmd)