{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Coordinate Reference System Management" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "xarray \"... is particularly tailored to working with netCDF files, which were the source of xarray’s data model...\" (http://xarray.pydata.org).\n", "\n", "For netCDF files, the GIS community uses CF conventions (http://cfconventions.org/).\n", "\n", "Additionally, GDAL also supports these attributes:\n", "\n", "- spatial_ref (Well Known Text)\n", "- GeoTransform (GeoTransform array)\n", "\n", "References:\n", "\n", "- Esri: https://pro.arcgis.com/en/pro-app/latest/help/data/multidimensional/spatial-reference-for-netcdf-data.htm\n", "- GDAL: https://gdal.org/drivers/raster/netcdf.html#georeference\n", "- pyproj: https://pyproj4.github.io/pyproj/stable/build_crs_cf.html\n", "\n", "Operations on xarray objects can cause data loss. Due to this, rioxarray writes and expects the spatial reference information to exist in the coordinates." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accessing the CRS object" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have opened a dataset and the Coordinate Reference System (CRS) can be determined, you can access it via the `rio.crs` accessor.\n", "\n", "#### Search order for the CRS (DataArray and Dataset):\n", "1. Look in attributes (`attrs`) of your data array for the `grid_mapping` coordinate name.\n", " Inside the `grid_mapping` coordinate first look for `spatial_ref` then `crs_wkt` and lastly the CF grid mapping attributes.\n", " This is in line with the Climate and Forecast (CF) conventions for storing the CRS as well as GDAL netCDF conventions.\n", "2. Look in the `crs` attribute and load in the CRS from there. This is for backwards compatibility with `xarray.open_rasterio`, which is deprecated since version 0.20.0. We recommend using `rioxarray.open_rasterio` instead.\n", "\n", "The value for the `crs` is anything accepted by `rasterio.crs.CRS.from_user_input()`\n", "\n", "#### Search order for the CRS for Dataset:\n", "If the CRS is not found using the search methods above, it also searches the `data_vars` and uses the\n", "first valid CRS found.\n", "\n", "#### decode_coords=\"all\"\n", "\n", "If you use one of xarray's open methods such as ``xarray.open_dataset`` to load netCDF files\n", "with the default engine, it is recommended to use `decode_coords=\"all\"`. This will load the grid mapping\n", "variable into coordinates for compatibility with rioxarray.\n", "\n", "#### API Documentation\n", "\n", "- [rio.write_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_crs)\n", "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)\n", "- [rio.estimate_utm_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.estimate_utm_crs)\n", "- [rio.set_spatial_dims()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.set_spatial_dims)\n", "- [rio.write_coordinate_system()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_coordinate_system)\n", "- [rio.write_transform()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_transform)\n", "- [rio.transform()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.transform)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rioxarray # activate the rio accessor\n", "import xarray\n", "from affine import Affine" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "rds = xarray.open_dataset(\"../../test/test_data/input/PLANET_SCOPE_3D.nc\", decode_coords=\"all\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'units': 'DN', 'nodata': 0.0}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rds.green.attrs" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'spatial_ref' ()>\n",
       "array(0)\n",
       "Coordinates:\n",
       "    spatial_ref  int64 0\n",
       "Attributes:\n",
       "    spatial_ref:  PROJCS["WGS 84 / UTM zone 22S",GEOGCS["WGS 84",DATUM["WGS_1...
" ], "text/plain": [ "\n", "array(0)\n", "Coordinates:\n", " spatial_ref int64 0\n", "Attributes:\n", " spatial_ref: PROJCS[\"WGS 84 / UTM zone 22S\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1..." ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rds.green.spatial_ref" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "CRS.from_epsg(32722)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rds.green.rio.crs" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Setting the CRS\n", "\n", "Use the `rio.write_crs` method to set the CRS on your `xarray.Dataset` or `xarray.DataArray`.\n", "This modifies the `xarray.Dataset` or `xarray.DataArray` and sets the CRS in a CF compliant manner.\n", "\n", "- [rio.write_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_crs)\n", "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'spatial_ref' ()>\n",
       "array(0)\n",
       "Coordinates:\n",
       "    spatial_ref  int64 0\n",
       "Attributes:\n",
       "    crs_wkt:                      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...\n",
       "    semi_major_axis:              6378137.0\n",
       "    semi_minor_axis:              6356752.314245179\n",
       "    inverse_flattening:           298.257223563\n",
       "    reference_ellipsoid_name:     WGS 84\n",
       "    longitude_of_prime_meridian:  0.0\n",
       "    prime_meridian_name:          Greenwich\n",
       "    geographic_crs_name:          WGS 84\n",
       "    grid_mapping_name:            latitude_longitude\n",
       "    spatial_ref:                  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...
" ], "text/plain": [ "\n", "array(0)\n", "Coordinates:\n", " spatial_ref int64 0\n", "Attributes:\n", " crs_wkt: GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"...\n", " semi_major_axis: 6378137.0\n", " semi_minor_axis: 6356752.314245179\n", " inverse_flattening: 298.257223563\n", " reference_ellipsoid_name: WGS 84\n", " longitude_of_prime_meridian: 0.0\n", " prime_meridian_name: Greenwich\n", " geographic_crs_name: WGS 84\n", " grid_mapping_name: latitude_longitude\n", " spatial_ref: GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"..." ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xda = xarray.DataArray(1)\n", "xda.rio.write_crs(4326, inplace=True)\n", "xda.spatial_ref" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "CRS.from_epsg(4326)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xda.rio.crs" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Spatial dimensions\n", "\n", "Only 1-dimensional X and Y dimensions are supported.\n", "\n", "The expected X/Y dimension names searched for in the `coords` are:\n", "\n", "- x | y\n", "- longitude | latitude\n", "- Coordinates (`coords`) with the CF attributes in `attrs`:\n", " - axis: X | Y\n", " - standard_name: longitude | latitude or projection_x_coordinate | projection_y_coordinate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Option 1: Write the CF attributes for non-standard dimension names\n", "\n", "If you don't want to rename your dimensions/coordinates,\n", "you can write the CF attributes so the coordinates can be found.\n", "\n", "- [rio.set_spatial_dims()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.set_spatial_dims)\n", "- [rio.write_coordinate_system()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_coordinate_system)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rds.rio.write_crs(\n", " 4326\n", " inplace=True,\n", ").rio.set_spatial_dims(\n", " x_dim=\"lon\",\n", " y_dim=\"lat\"\n", " inplace=True,\n", ").rio.write_coordinate_system(inplace=True)" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "Option 2: Rename your coordinates\n", "\n", "[xarray.Dataset.rename](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.rename.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rds = rds.rename(lon=longitute, lat=latitude) " ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Setting the transform of the dataset\n", "\n", "The transform can be calculated from the coordinates of your data.\n", "This method is useful if your netCDF file does not have coordinates present.\n", "Use the `rio.write_transform` method to set the transform on your `xarray.Dataset` or `xarray.DataArray`.\n", "\n", "- [rio.write_transform()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_transform)\n", "- [rio.transform()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.transform)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'466266.0 3.0 0.0 8084700.0 0.0 -3.0'" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transform = Affine(3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0)\n", "xda.rio.write_transform(transform, inplace=True)\n", "xda.spatial_ref.GeoTransform" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Affine(3.0, 0.0, 466266.0,\n", " 0.0, -3.0, 8084700.0)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xda.rio.transform()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.4" } }, "nbformat": 4, "nbformat_minor": 4 }