diff --git a/.gitignore b/.gitignore index 0951de11..2fd6749f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +*.bak + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 239a9732..b745e730 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -10,7 +10,7 @@ repos: hooks: - id: black - repo: https://github.com/timothycrosley/isort - rev: 5.10.1 + rev: 5.12.0 hooks: - id: isort args: [rioxarray/, test/, docs/] diff --git a/docs/examples/crs_index.ipynb b/docs/examples/crs_index.ipynb new file mode 100644 index 00000000..1985d6fd --- /dev/null +++ b/docs/examples/crs_index.ipynb @@ -0,0 +1,9170 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# CRS Index" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default rioxarray stores CRS information as a scalar coordinate in Xarray. This works well because the important metadata persists through all Xarray operations and when it comes time to save an output, we must retrieve this important metadata.\n", + "\n", + "As of Xarray v2022.09.0 it is possible to create custom Indexes (typically a `PandasIndex` is used behind the scenes that determines how dimensional coordinates are queried and aligned). With the new Xindexes interface, we can attach CRS as persistent metadata to a custom `CRSIndex`. Effectively we enhance the default `PandasIndex` by adding associated metadata (CRS) and new functionality like checking CRS compatibility across objects.\n", + "\n", + "This notebook showcases experimental CRSIndex functionality" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import rioxarray\n", + "import geopandas as gpd\n", + "import xvec\n", + "\n", + "xr.set_options(\n", + " display_expand_data=False,\n", + " display_expand_indexes=True,\n", + ");\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'band_data' (band: 1, y: 500, x: 500)>\n",
+       "[250000 values with dtype=float32]\n",
+       "Coordinates:\n",
+       "  * band         (band) int64 1\n",
+       "  * x            (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n",
+       "  * y            (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n",
+       "    spatial_ref  int64 ...\n",
+       "Attributes:\n",
+       "    AREA_OR_POINT:       Area\n",
+       "    OVR_RESAMPLING_ALG:  NEAREST
" + ], + "text/plain": [ + "\n", + "[250000 values with dtype=float32]\n", + "Coordinates:\n", + " * band (band) int64 1\n", + " * x (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n", + " * y (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n", + " spatial_ref int64 ...\n", + "Attributes:\n", + " AREA_OR_POINT: Area\n", + " OVR_RESAMPLING_ALG: NEAREST" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# DataArray\n", + "DA = xr.open_dataarray(\"../../test/test_data/input/cog.tif\", engine=\"rasterio\")\n", + "DA" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note `x` and `y` dimensions are backed by a corresponding separate `PandasIndex` by default. `spatial_ref` is a non-dimensional scalar coordinate and therefore does not appear in the `Indexes`" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "CRS.from_wkt('PROJCS[\"Albers\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378140,298.256999999996,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Albers_Conic_Equal_Area\"],PARAMETER[\"latitude_of_center\",23],PARAMETER[\"longitude_of_center\",-96],PARAMETER[\"standard_parallel_1\",29.5],PARAMETER[\"standard_parallel_2\",45.5],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]')" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DA.rio.crs" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'band_data' (band: 1, y: 500, x: 500)>\n",
+       "[250000 values with dtype=float32]\n",
+       "Coordinates:\n",
+       "  * band         (band) int64 1\n",
+       "    spatial_ref  int64 ...\n",
+       "  * x            (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n",
+       "  * y            (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=PROJCS["Albers",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS  ...)\n",
+       "    y            CRSIndex (crs=PROJCS["Albers",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS  ...)\n",
+       "Attributes:\n",
+       "    AREA_OR_POINT:       Area\n",
+       "    OVR_RESAMPLING_ALG:  NEAREST
" + ], + "text/plain": [ + "\n", + "[250000 values with dtype=float32]\n", + "Coordinates:\n", + " * band (band) int64 1\n", + " spatial_ref int64 ...\n", + " * x (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n", + " * y (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n", + "Indexes:\n", + " x CRSIndex (crs=PROJCS[\"Albers\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS ...)\n", + " y CRSIndex (crs=PROJCS[\"Albers\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS ...)\n", + "Attributes:\n", + " AREA_OR_POINT: Area\n", + " OVR_RESAMPLING_ALG: NEAREST" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# https://corteva.github.io/rioxarray/stable/rioxarray.html#rioxarray-set-options\n", + "with rioxarray.set_options(use_crs_index=True):\n", + " DAcrs = xr.open_dataarray(\"../../test/test_data/input/cog.tif\", engine=\"rasterio\")\n", + "DAcrs" + ] + }, + { + "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 'band_data' (band: 1, y: 500, x: 500)>\n",
+       "[250000 values with dtype=float32]\n",
+       "Coordinates:\n",
+       "  * band         (band) int64 1\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n",
+       "  * y            (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=PROJCS["Albers",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS  ...)\n",
+       "    y            CRSIndex (crs=PROJCS["Albers",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS  ...)\n",
+       "Attributes:\n",
+       "    AREA_OR_POINT:       Area\n",
+       "    OVR_RESAMPLING_ALG:  NEAREST
" + ], + "text/plain": [ + "\n", + "[250000 values with dtype=float32]\n", + "Coordinates:\n", + " * band (band) int64 1\n", + " spatial_ref int64 0\n", + " * x (x) float64 1.635e+06 1.635e+06 ... 1.784e+06 1.784e+06\n", + " * y (y) float64 2.715e+06 2.714e+06 ... 2.565e+06 2.565e+06\n", + "Indexes:\n", + " x CRSIndex (crs=PROJCS[\"Albers\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS ...)\n", + " y CRSIndex (crs=PROJCS[\"Albers\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS ...)\n", + "Attributes:\n", + " AREA_OR_POINT: Area\n", + " OVR_RESAMPLING_ALG: NEAREST" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Or alternatively configure index after opening:\n", + "with rioxarray.set_options(use_crs_index=True):\n", + " DAcrs = DA.rio.write_crs(DA.rio.crs)\n", + "DAcrs " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "rasterio.crs.CRS" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(DAcrs.xindexes['x'].crs)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the spatial relationship between `x` and `y` is indicated by their shared `CRSIndex`, which includes an abbreviated representation (often an EPSG code). The CRS is actually a `rasterio.crs.CRS` object with attributes and methods that propoagates through Xarray operations." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
+       "    spatial_ref  int64 ...\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 ...\n",
+       "    green        (time, y, x) float64 ...\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 10, y: 10)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n", + " * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n", + " spatial_ref int64 ...\n", + "Data variables:\n", + " blue (time, y, x) float64 ...\n", + " green (time, y, x) float64 ...\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DS = xr.open_dataset(\"../../test/test_data/input/PLANET_SCOPE_3D.nc\", engine=\"rasterio\")\n", + "DS" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "CRS.from_epsg(32722)" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DS.rio.crs" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 ...\n",
+       "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 ...\n",
+       "    green        (time, y, x) float64 ...\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:32722)\n",
+       "    y            CRSIndex (crs=EPSG:32722)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 10, y: 10)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 ...\n", + " * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n", + " * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n", + "Data variables:\n", + " blue (time, y, x) float64 ...\n", + " green (time, y, x) float64 ...\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:32722)\n", + " y CRSIndex (crs=EPSG:32722)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "with rioxarray.set_options(use_crs_index=True):\n", + " DSutm = xr.open_dataset(\"../../test/test_data/input/PLANET_SCOPE_3D.nc\", engine=\"rasterio\")\n", + "\n", + "DSutm" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "CRS.from_epsg(32722)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# CRS information is retrieved via the `.rio` accessor (In this case UTM https://epsg.io/32722)\n", + "DSutm.rio.crs" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 ...\n",
+       "    green        (time, y, x) float64 ...\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:32722)\n",
+       "    y            CRSIndex (crs=EPSG:32722)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 10, y: 10)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n", + " * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n", + "Data variables:\n", + " blue (time, y, x) float64 ...\n", + " green (time, y, x) float64 ...\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:32722)\n", + " y CRSIndex (crs=EPSG:32722)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Or Convert from existing PandasIndex to CRSIndex\n", + "with rioxarray.set_options(use_crs_index=True):\n", + " DSutm = DS.rio.write_crs(DS.rio.crs)\n", + "DSutm" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "CRS.from_epsg(32722)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Access CRS info the same way\n", + "DSutm.rio.crs" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "CRS.from_epsg(32722)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Or access via xindexes\n", + "DSutm.xindexes['x'].crs" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the spatial relationship between `x` and `y` is indicated by their shared `CRSIndex`, which includes an abbreviated representation (often an EPSG code)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using CRSIndex" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (x: 10, y: 10, time: 2)\n",
+       "Coordinates:\n",
+       "  * x            (x) float64 -51.32 -51.32 -51.32 ... -51.32 -51.32 -51.32\n",
+       "  * y            (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n",
+       "    green        (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (x: 10, y: 10, time: 2)\n", + "Coordinates:\n", + " * x (x) float64 -51.32 -51.32 -51.32 ... -51.32 -51.32 -51.32\n", + " * y (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + "Data variables:\n", + " blue (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n", + " green (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DSll = DS.rio.reproject('EPSG:4326')\n", + "DSll" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
+       "    spatial_ref  int64 ...\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n",
+       "    green        (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 10, y: 10)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n", + " * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n", + " spatial_ref int64 ...\n", + "Data variables:\n", + " blue (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n", + " green (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DS" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 0, y: 0)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "  * x            (x) float64 \n",
+       "  * y            (y) float64 \n",
+       "    spatial_ref  int64 0\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 \n",
+       "    green        (time, y, x) float64
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 0, y: 0)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " * x (x) float64 \n", + " * y (y) float64 \n", + " spatial_ref int64 0\n", + "Data variables:\n", + " blue (time, y, x) float64 \n", + " green (time, y, x) float64 " + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DS + DSll" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With default PandasIndex Xarray will happily align coordinates assuming they belong to the same CRS. In this case, none of the coordinates match (Lat/Lon versus UTM)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (x: 10, y: 10, time: 2)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 -51.32 -51.32 -51.32 ... -51.32 -51.32 -51.32\n",
+       "  * y            (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n",
+       "    green        (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:4326)\n",
+       "    y            CRSIndex (crs=EPSG:4326)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (x: 10, y: 10, time: 2)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * x (x) float64 -51.32 -51.32 -51.32 ... -51.32 -51.32 -51.32\n", + " * y (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n", + "Data variables:\n", + " blue (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n", + " green (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:4326)\n", + " y CRSIndex (crs=EPSG:4326)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NOTE: if methods create a new DataArray behind the scenes (like rio.reproject) have to chain with write_crs():\n", + "with rioxarray.set_options(use_crs_index=True):\n", + " DSll = DSutm.rio.reproject('EPSG:4326').rio.write_crs('EPSG:4326')\n", + "DSll" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "CRS mismatch between the CRS of index geometries and the CRS of input geometries.\nIndex CRS: EPSG:32722\nInput CRS: EPSG:4326\n", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[19], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m DSutm \u001b[39m+\u001b[39;49m DSll\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/_typed_ops.py:16\u001b[0m, in \u001b[0;36mDatasetOpsMixin.__add__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__add__\u001b[39m(\u001b[39mself\u001b[39m, other):\n\u001b[0;32m---> 16\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_binary_op(other, operator\u001b[39m.\u001b[39;49madd)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:6604\u001b[0m, in \u001b[0;36mDataset._binary_op\u001b[0;34m(self, other, f, reflexive, join)\u001b[0m\n\u001b[1;32m 6602\u001b[0m align_type \u001b[39m=\u001b[39m OPTIONS[\u001b[39m\"\u001b[39m\u001b[39marithmetic_join\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39mif\u001b[39;00m join \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39melse\u001b[39;00m join\n\u001b[1;32m 6603\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(other, (DataArray, Dataset)):\n\u001b[0;32m-> 6604\u001b[0m \u001b[39mself\u001b[39m, other \u001b[39m=\u001b[39m align(\u001b[39mself\u001b[39;49m, other, join\u001b[39m=\u001b[39;49malign_type, copy\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m) \u001b[39m# type: ignore[assignment]\u001b[39;00m\n\u001b[1;32m 6605\u001b[0m g \u001b[39m=\u001b[39m f \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m reflexive \u001b[39melse\u001b[39;00m \u001b[39mlambda\u001b[39;00m x, y: f(y, x)\n\u001b[1;32m 6606\u001b[0m ds \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_calculate_binary_op(g, other, join\u001b[39m=\u001b[39malign_type)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:787\u001b[0m, in \u001b[0;36malign\u001b[0;34m(join, copy, indexes, exclude, fill_value, *objects)\u001b[0m\n\u001b[1;32m 591\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 592\u001b[0m \u001b[39mGiven any number of Dataset and/or DataArray objects, returns new\u001b[39;00m\n\u001b[1;32m 593\u001b[0m \u001b[39mobjects with aligned indexes and dimension sizes.\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 777\u001b[0m \n\u001b[1;32m 778\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 779\u001b[0m aligner \u001b[39m=\u001b[39m Aligner(\n\u001b[1;32m 780\u001b[0m objects,\n\u001b[1;32m 781\u001b[0m join\u001b[39m=\u001b[39mjoin,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 785\u001b[0m fill_value\u001b[39m=\u001b[39mfill_value,\n\u001b[1;32m 786\u001b[0m )\n\u001b[0;32m--> 787\u001b[0m aligner\u001b[39m.\u001b[39;49malign()\n\u001b[1;32m 788\u001b[0m \u001b[39mreturn\u001b[39;00m aligner\u001b[39m.\u001b[39mresults\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:572\u001b[0m, in \u001b[0;36mAligner.align\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 570\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfind_matching_unindexed_dims()\n\u001b[1;32m 571\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39massert_no_index_conflict()\n\u001b[0;32m--> 572\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49malign_indexes()\n\u001b[1;32m 573\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39massert_unindexed_dim_sizes_equal()\n\u001b[1;32m 575\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mjoin \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39moverride\u001b[39m\u001b[39m\"\u001b[39m:\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:424\u001b[0m, in \u001b[0;36mAligner.align_indexes\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 417\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 418\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mcannot align objects with join=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mexact\u001b[39m\u001b[39m'\u001b[39m\u001b[39m where \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 419\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mindex/labels/sizes are not equal along \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 420\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mthese coordinates (dimensions): \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 421\u001b[0m \u001b[39m+\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m, \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mjoin(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mname\u001b[39m!r}\u001b[39;00m\u001b[39m \u001b[39m\u001b[39m{\u001b[39;00mdims\u001b[39m!r}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m name, dims \u001b[39min\u001b[39;00m key[\u001b[39m0\u001b[39m])\n\u001b[1;32m 422\u001b[0m )\n\u001b[1;32m 423\u001b[0m joiner \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_get_index_joiner(index_cls)\n\u001b[0;32m--> 424\u001b[0m joined_index \u001b[39m=\u001b[39m joiner(matching_indexes)\n\u001b[1;32m 425\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mjoin \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39mleft\u001b[39m\u001b[39m\"\u001b[39m:\n\u001b[1;32m 426\u001b[0m joined_index_vars \u001b[39m=\u001b[39m matching_index_vars[\u001b[39m0\u001b[39m]\n", + "File \u001b[0;32m~/GitHub/rioxarray/rioxarray/_crs_index.py:150\u001b[0m, in \u001b[0;36mCRSIndex.join\u001b[0;34m(self, other, how)\u001b[0m\n\u001b[1;32m 148\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mjoin\u001b[39m(\u001b[39mself\u001b[39m, other, how\u001b[39m=\u001b[39m\u001b[39m\"\u001b[39m\u001b[39minner\u001b[39m\u001b[39m\"\u001b[39m):\n\u001b[1;32m 149\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_check_crs(other\u001b[39m.\u001b[39mcrs, allow_none\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m):\n\u001b[0;32m--> 150\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_crs_mismatch_raise(other\u001b[39m.\u001b[39;49mcrs)\n\u001b[1;32m 152\u001b[0m new_indexes \u001b[39m=\u001b[39m {\n\u001b[1;32m 153\u001b[0m k: v\u001b[39m.\u001b[39mjoin(other\u001b[39m.\u001b[39m_indexes[k], how\u001b[39m=\u001b[39mhow) \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_indexes\u001b[39m.\u001b[39mitems()\n\u001b[1;32m 154\u001b[0m }\n\u001b[1;32m 156\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mtype\u001b[39m(\u001b[39mself\u001b[39m)(new_indexes, crs\u001b[39m=\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcrs)\n", + "File \u001b[0;32m~/GitHub/rioxarray/rioxarray/_crs_index.py:137\u001b[0m, in \u001b[0;36mCRSIndex._crs_mismatch_raise\u001b[0;34m(self, other_crs, stacklevel)\u001b[0m\n\u001b[1;32m 129\u001b[0m msg \u001b[39m=\u001b[39m (\n\u001b[1;32m 130\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mCRS mismatch between the CRS of index geometries \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 131\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mand the CRS of input geometries.\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 132\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mIndex CRS: \u001b[39m\u001b[39m{\u001b[39;00msrs\u001b[39m}\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 133\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mInput CRS: \u001b[39m\u001b[39m{\u001b[39;00mother_srs\u001b[39m}\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 134\u001b[0m )\n\u001b[1;32m 136\u001b[0m \u001b[39mif\u001b[39;00m get_option(CRS_INDEX_ERROR):\n\u001b[0;32m--> 137\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(msg)\n\u001b[1;32m 138\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m 139\u001b[0m warnings\u001b[39m.\u001b[39mwarn(msg, \u001b[39mUserWarning\u001b[39;00m, stacklevel\u001b[39m=\u001b[39mstacklevel)\n", + "\u001b[0;31mValueError\u001b[0m: CRS mismatch between the CRS of index geometries and the CRS of input geometries.\nIndex CRS: EPSG:32722\nInput CRS: EPSG:4326\n" + ] + } + ], + "source": [ + "DSutm + DSll" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/scott/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:424: UserWarning: CRS mismatch between the CRS of index geometries and the CRS of input geometries.\n", + "Index CRS: EPSG:32722\n", + "Input CRS: EPSG:4326\n", + "\n", + " joined_index = joiner(matching_indexes)\n", + "/Users/scott/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:506: UserWarning: CRS mismatch between the CRS of index geometries and the CRS of input geometries.\n", + "Index CRS: EPSG:4326\n", + "Input CRS: EPSG:32722\n", + "\n", + " indexers = obj_idx.reindex_like(aligned_idx, **self.reindex_kwargs)\n" + ] + } + ], + "source": [ + "# NOTE: can have a CRS mismatch warning instead of raising an error\n", + "with rioxarray.set_options(crs_index_error=False):\n", + " DSutm + DSll" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "cannot re-index or align objects with conflicting indexes found for the following coordinates: 'x' (2 conflicting indexes), 'y' (2 conflicting indexes)\nConflicting indexes may occur when\n- they relate to different sets of coordinate and/or dimension names\n- they don't have the same type\n- they may be used to reindex data along common dimensions", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[21], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[39m# Xarray raises a ValueError for alignment of defaut PandasIndex with new CRSIndex... we might want to allow this?\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m DSutm \u001b[39m+\u001b[39;49m DS\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/_typed_ops.py:16\u001b[0m, in \u001b[0;36mDatasetOpsMixin.__add__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__add__\u001b[39m(\u001b[39mself\u001b[39m, other):\n\u001b[0;32m---> 16\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_binary_op(other, operator\u001b[39m.\u001b[39;49madd)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:6604\u001b[0m, in \u001b[0;36mDataset._binary_op\u001b[0;34m(self, other, f, reflexive, join)\u001b[0m\n\u001b[1;32m 6602\u001b[0m align_type \u001b[39m=\u001b[39m OPTIONS[\u001b[39m\"\u001b[39m\u001b[39marithmetic_join\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39mif\u001b[39;00m join \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39melse\u001b[39;00m join\n\u001b[1;32m 6603\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(other, (DataArray, Dataset)):\n\u001b[0;32m-> 6604\u001b[0m \u001b[39mself\u001b[39m, other \u001b[39m=\u001b[39m align(\u001b[39mself\u001b[39;49m, other, join\u001b[39m=\u001b[39;49malign_type, copy\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m) \u001b[39m# type: ignore[assignment]\u001b[39;00m\n\u001b[1;32m 6605\u001b[0m g \u001b[39m=\u001b[39m f \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m reflexive \u001b[39melse\u001b[39;00m \u001b[39mlambda\u001b[39;00m x, y: f(y, x)\n\u001b[1;32m 6606\u001b[0m ds \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_calculate_binary_op(g, other, join\u001b[39m=\u001b[39malign_type)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:787\u001b[0m, in \u001b[0;36malign\u001b[0;34m(join, copy, indexes, exclude, fill_value, *objects)\u001b[0m\n\u001b[1;32m 591\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 592\u001b[0m \u001b[39mGiven any number of Dataset and/or DataArray objects, returns new\u001b[39;00m\n\u001b[1;32m 593\u001b[0m \u001b[39mobjects with aligned indexes and dimension sizes.\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 777\u001b[0m \n\u001b[1;32m 778\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 779\u001b[0m aligner \u001b[39m=\u001b[39m Aligner(\n\u001b[1;32m 780\u001b[0m objects,\n\u001b[1;32m 781\u001b[0m join\u001b[39m=\u001b[39mjoin,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 785\u001b[0m fill_value\u001b[39m=\u001b[39mfill_value,\n\u001b[1;32m 786\u001b[0m )\n\u001b[0;32m--> 787\u001b[0m aligner\u001b[39m.\u001b[39;49malign()\n\u001b[1;32m 788\u001b[0m \u001b[39mreturn\u001b[39;00m aligner\u001b[39m.\u001b[39mresults\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:571\u001b[0m, in \u001b[0;36mAligner.align\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 569\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfind_matching_indexes()\n\u001b[1;32m 570\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfind_matching_unindexed_dims()\n\u001b[0;32m--> 571\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49massert_no_index_conflict()\n\u001b[1;32m 572\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39malign_indexes()\n\u001b[1;32m 573\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39massert_unindexed_dim_sizes_equal()\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:314\u001b[0m, in \u001b[0;36mAligner.assert_no_index_conflict\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 310\u001b[0m \u001b[39mif\u001b[39;00m dup:\n\u001b[1;32m 311\u001b[0m items_msg \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m, \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mjoin(\n\u001b[1;32m 312\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mk\u001b[39m!r}\u001b[39;00m\u001b[39m (\u001b[39m\u001b[39m{\u001b[39;00mv\u001b[39m}\u001b[39;00m\u001b[39m conflicting indexes)\u001b[39m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m dup\u001b[39m.\u001b[39mitems()\n\u001b[1;32m 313\u001b[0m )\n\u001b[0;32m--> 314\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 315\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mcannot re-index or align objects with conflicting indexes found for \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 316\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mthe following \u001b[39m\u001b[39m{\u001b[39;00mmsg\u001b[39m}\u001b[39;00m\u001b[39m: \u001b[39m\u001b[39m{\u001b[39;00mitems_msg\u001b[39m}\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 317\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mConflicting indexes may occur when\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 318\u001b[0m \u001b[39m\"\u001b[39m\u001b[39m- they relate to different sets of coordinate and/or dimension names\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 319\u001b[0m \u001b[39m\"\u001b[39m\u001b[39m- they don\u001b[39m\u001b[39m'\u001b[39m\u001b[39mt have the same type\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 320\u001b[0m \u001b[39m\"\u001b[39m\u001b[39m- they may be used to reindex data along common dimensions\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 321\u001b[0m )\n", + "\u001b[0;31mValueError\u001b[0m: cannot re-index or align objects with conflicting indexes found for the following coordinates: 'x' (2 conflicting indexes), 'y' (2 conflicting indexes)\nConflicting indexes may occur when\n- they relate to different sets of coordinate and/or dimension names\n- they don't have the same type\n- they may be used to reindex data along common dimensions" + ] + } + ], + "source": [ + "# Xarray raises a ValueError for alignment of defaut PandasIndex with new CRSIndex... we might want to allow this?\n", + "DSutm + DS" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " {'x': [-51.3173]} {'method': 'nearest', 'tolerance': None}\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'blue' (time: 2, y: 10, x: 1)>\n",
+       "5.078 2.736 5.267 2.344 5.195 5.714 6.407 ... 3.473 2.622 5.093 4.75 3.218 1.888\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 -51.32\n",
+       "  * y            (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:4326)\n",
+       "    y            CRSIndex (crs=EPSG:4326)\n",
+       "Attributes:\n",
+       "    NETCDF_DIM_EXTRA:        {time}\n",
+       "    NETCDF_DIM_time_DEF:     [2. 6.]\n",
+       "    NETCDF_DIM_time_VALUES:  [     0.       872712.659688]\n",
+       "    units:                   ('DN', 'DN')
" + ], + "text/plain": [ + "\n", + "5.078 2.736 5.267 2.344 5.195 5.714 6.407 ... 3.473 2.622 5.093 4.75 3.218 1.888\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * x (x) float64 -51.32\n", + " * y (y) float64 -17.32 -17.32 -17.32 ... -17.32 -17.32 -17.32\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:4326)\n", + " y CRSIndex (crs=EPSG:4326)\n", + "Attributes:\n", + " NETCDF_DIM_EXTRA: {time}\n", + " NETCDF_DIM_time_DEF: [2. 6.]\n", + " NETCDF_DIM_time_VALUES: [ 0. 872712.659688]\n", + " units: ('DN', 'DN')" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# vector selection retains CRS\n", + "DSll.blue.sel(x=[-51.3173], method='nearest')" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/scott/GitHub/rioxarray/rioxarray/rioxarray.py:522: UserWarning: MissingSpatialDimensionError: x or y dimensions required for CRSIndex. See 'rio.set_spatial_dims()'\n", + " warnings.warn(\"MissingSpatialDimensionError: x or y dimensions required for CRSIndex. See 'rio.set_spatial_dims()'\")\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray (x: 2)>\n",
+       "nan nan\n",
+       "Coordinates:\n",
+       "  * x            (x) float64 -51.32 -51.32\n",
+       "    spatial_ref  int64 0
" + ], + "text/plain": [ + "\n", + "nan nan\n", + "Coordinates:\n", + " * x (x) float64 -51.32 -51.32\n", + " spatial_ref int64 0" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "with rioxarray.set_options(use_crs_index=True):\n", + " xvector = xr.DataArray(coords=dict(x=[-51.31734, -51.31745]), \n", + " dims=[\"x\"]).rio.write_crs('EPSG:4326')\n", + "\n", + "xvector" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray (x: 2, y: 2)>\n",
+       "nan nan nan nan\n",
+       "Coordinates:\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 -51.32 -51.32\n",
+       "  * y            (y) float64 -17.32 -17.32\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:4326)\n",
+       "    y            CRSIndex (crs=EPSG:4326)
" + ], + "text/plain": [ + "\n", + "nan nan nan nan\n", + "Coordinates:\n", + " spatial_ref int64 0\n", + " * x (x) float64 -51.32 -51.32\n", + " * y (y) float64 -17.32 -17.32\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:4326)\n", + " y CRSIndex (crs=EPSG:4326)" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# We can ensure we sample or interpolate with points in the same CRS\n", + "with rioxarray.set_options(use_crs_index=True):\n", + " points = xr.DataArray(coords=dict(x=[-51.31734, -51.31745], y=[-17.3228, -17.3229]), \n", + " dims=[\"x\",\"y\"]).rio.write_crs('EPSG:4326')\n", + "points" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray (x: 2, y: 2)>\n",
+       "nan nan nan nan\n",
+       "Coordinates:\n",
+       "    spatial_ref  int64 0\n",
+       "  * x            (x) float64 -51.32 -51.32\n",
+       "  * y            (y) float64 -17.32 -17.32\n",
+       "Indexes:\n",
+       "    x            CRSIndex (crs=EPSG:4326)\n",
+       "    y            CRSIndex (crs=EPSG:4326)
" + ], + "text/plain": [ + "\n", + "nan nan nan nan\n", + "Coordinates:\n", + " spatial_ref int64 0\n", + " * x (x) float64 -51.32 -51.32\n", + " * y (y) float64 -17.32 -17.32\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:4326)\n", + " y CRSIndex (crs=EPSG:4326)" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "points" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "cannot remove index(es) 'y', which would corrupt the following index built from coordinates 'x', 'y':\n", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[26], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m points\u001b[39m.\u001b[39;49mdrop_indexes(\u001b[39m'\u001b[39;49m\u001b[39my\u001b[39;49m\u001b[39m'\u001b[39;49m) \u001b[39m# ValueError: cannot remove index(es) 'y', which would corrupt the following index built from coordinates 'x', 'y':\u001b[39;00m\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataarray.py:3024\u001b[0m, in \u001b[0;36mDataArray.drop_indexes\u001b[0;34m(self, coord_names, errors)\u001b[0m\n\u001b[1;32m 3002\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mdrop_indexes\u001b[39m(\n\u001b[1;32m 3003\u001b[0m \u001b[39mself\u001b[39m: T_DataArray,\n\u001b[1;32m 3004\u001b[0m coord_names: Hashable \u001b[39m|\u001b[39m Iterable[Hashable],\n\u001b[1;32m 3005\u001b[0m \u001b[39m*\u001b[39m,\n\u001b[1;32m 3006\u001b[0m errors: ErrorOptions \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39mraise\u001b[39m\u001b[39m\"\u001b[39m,\n\u001b[1;32m 3007\u001b[0m ) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m T_DataArray:\n\u001b[1;32m 3008\u001b[0m \u001b[39m \u001b[39m\u001b[39m\"\"\"Drop the indexes assigned to the given coordinates.\u001b[39;00m\n\u001b[1;32m 3009\u001b[0m \n\u001b[1;32m 3010\u001b[0m \u001b[39m Parameters\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 3022\u001b[0m \u001b[39m A new dataarray with dropped indexes.\u001b[39;00m\n\u001b[1;32m 3023\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 3024\u001b[0m ds \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_to_temp_dataset()\u001b[39m.\u001b[39;49mdrop_indexes(coord_names, errors\u001b[39m=\u001b[39;49merrors)\n\u001b[1;32m 3025\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_from_temp_dataset(ds)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:5163\u001b[0m, in \u001b[0;36mDataset.drop_indexes\u001b[0;34m(self, coord_names, errors)\u001b[0m\n\u001b[1;32m 5158\u001b[0m \u001b[39mif\u001b[39;00m unindexed_coords:\n\u001b[1;32m 5159\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 5160\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mthose coordinates do not have an index: \u001b[39m\u001b[39m{\u001b[39;00munindexed_coords\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 5161\u001b[0m )\n\u001b[0;32m-> 5163\u001b[0m assert_no_index_corrupted(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mxindexes, coord_names, action\u001b[39m=\u001b[39;49m\u001b[39m\"\u001b[39;49m\u001b[39mremove index(es)\u001b[39;49m\u001b[39m\"\u001b[39;49m)\n\u001b[1;32m 5165\u001b[0m variables \u001b[39m=\u001b[39m {}\n\u001b[1;32m 5166\u001b[0m \u001b[39mfor\u001b[39;00m name, var \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_variables\u001b[39m.\u001b[39mitems():\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/indexes.py:1526\u001b[0m, in \u001b[0;36massert_no_index_corrupted\u001b[0;34m(indexes, coord_names, action)\u001b[0m\n\u001b[1;32m 1524\u001b[0m common_names_str \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m, \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mjoin(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mk\u001b[39m!r}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m k \u001b[39min\u001b[39;00m common_names)\n\u001b[1;32m 1525\u001b[0m index_names_str \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m, \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mjoin(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mk\u001b[39m!r}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m k \u001b[39min\u001b[39;00m index_coords)\n\u001b[0;32m-> 1526\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 1527\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mcannot \u001b[39m\u001b[39m{\u001b[39;00maction\u001b[39m}\u001b[39;00m\u001b[39m \u001b[39m\u001b[39m{\u001b[39;00mcommon_names_str\u001b[39m}\u001b[39;00m\u001b[39m, which would corrupt \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1528\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mthe following index built from coordinates \u001b[39m\u001b[39m{\u001b[39;00mindex_names_str\u001b[39m}\u001b[39;00m\u001b[39m:\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1529\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mindex\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 1530\u001b[0m )\n", + "\u001b[0;31mValueError\u001b[0m: cannot remove index(es) 'y', which would corrupt the following index built from coordinates 'x', 'y':\n" + ] + } + ], + "source": [ + "points.drop_indexes('y') # ValueError: cannot remove index(es) 'y', which would corrupt the following index built from coordinates 'x', 'y':" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'x' (x: 2)>\n",
+       "-51.32 -51.32\n",
+       "Coordinates:\n",
+       "    spatial_ref  int64 0\n",
+       "    x            (x) float64 -51.32 -51.32
" + ], + "text/plain": [ + "\n", + "-51.32 -51.32\n", + "Coordinates:\n", + " spatial_ref int64 0\n", + " x (x) float64 -51.32 -51.32" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "points.x" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " {'x': \n", + "nan nan nan nan\n", + "Coordinates:\n", + " spatial_ref int64 0\n", + " * x (x) float64 -51.32 -51.32\n", + " * y (y) float64 -17.32 -17.32\n", + "Indexes:\n", + " x CRSIndex (crs=EPSG:4326)\n", + " y CRSIndex (crs=EPSG:4326)} {'method': 'nearest', 'tolerance': None}\n" + ] + }, + { + "ename": "IndexError", + "evalue": "Dimensions of indexers mismatch: (slice(None, None, None), slice(None, None, None), \narray([[9, 9],\n [9, 9]]))", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:847\u001b[0m, in \u001b[0;36mVariable._broadcast_indexes_vectorized\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 846\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 847\u001b[0m variables \u001b[39m=\u001b[39m _broadcast_compat_variables(\u001b[39m*\u001b[39;49mvariables)\n\u001b[1;32m 848\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m:\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:3121\u001b[0m, in \u001b[0;36m_broadcast_compat_variables\u001b[0;34m(*variables)\u001b[0m\n\u001b[1;32m 3116\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"Create broadcast compatible variables, with the same dimensions.\u001b[39;00m\n\u001b[1;32m 3117\u001b[0m \n\u001b[1;32m 3118\u001b[0m \u001b[39mUnlike the result of broadcast_variables(), some variables may have\u001b[39;00m\n\u001b[1;32m 3119\u001b[0m \u001b[39mdimensions of size 1 instead of the size of the broadcast dimension.\u001b[39;00m\n\u001b[1;32m 3120\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m-> 3121\u001b[0m dims \u001b[39m=\u001b[39m \u001b[39mtuple\u001b[39m(_unified_dims(variables))\n\u001b[1;32m 3122\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mtuple\u001b[39m(var\u001b[39m.\u001b[39mset_dims(dims) \u001b[39mif\u001b[39;00m var\u001b[39m.\u001b[39mdims \u001b[39m!=\u001b[39m dims \u001b[39melse\u001b[39;00m var \u001b[39mfor\u001b[39;00m var \u001b[39min\u001b[39;00m variables)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:3108\u001b[0m, in \u001b[0;36m_unified_dims\u001b[0;34m(variables)\u001b[0m\n\u001b[1;32m 3107\u001b[0m \u001b[39melif\u001b[39;00m all_dims[d] \u001b[39m!=\u001b[39m s:\n\u001b[0;32m-> 3108\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 3109\u001b[0m \u001b[39m\"\u001b[39m\u001b[39moperands cannot be broadcast together \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 3110\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mwith mismatched lengths for dimension \u001b[39m\u001b[39m{\u001b[39;00md\u001b[39m!r}\u001b[39;00m\u001b[39m: \u001b[39m\u001b[39m{\u001b[39;00m(all_dims[d],\u001b[39m \u001b[39ms)\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 3111\u001b[0m )\n\u001b[1;32m 3112\u001b[0m \u001b[39mreturn\u001b[39;00m all_dims\n", + "\u001b[0;31mValueError\u001b[0m: operands cannot be broadcast together with mismatched lengths for dimension 'y': (10, 2)", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[28], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[39m# So no CRS check against selection vectors currently:\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m DSll\u001b[39m.\u001b[39;49msel(x\u001b[39m=\u001b[39;49mpoints, method\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mnearest\u001b[39;49m\u001b[39m'\u001b[39;49m)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:2594\u001b[0m, in \u001b[0;36mDataset.sel\u001b[0;34m(self, indexers, method, tolerance, drop, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 2591\u001b[0m query_results\u001b[39m.\u001b[39mdrop_coords\u001b[39m.\u001b[39mappend(k)\n\u001b[1;32m 2592\u001b[0m query_results\u001b[39m.\u001b[39mvariables \u001b[39m=\u001b[39m no_scalar_variables\n\u001b[0;32m-> 2594\u001b[0m result \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49misel(indexers\u001b[39m=\u001b[39;49mquery_results\u001b[39m.\u001b[39;49mdim_indexers, drop\u001b[39m=\u001b[39;49mdrop)\n\u001b[1;32m 2595\u001b[0m \u001b[39mreturn\u001b[39;00m result\u001b[39m.\u001b[39m_overwrite_indexes(\u001b[39m*\u001b[39mquery_results\u001b[39m.\u001b[39mas_tuple()[\u001b[39m1\u001b[39m:])\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:2434\u001b[0m, in \u001b[0;36mDataset.isel\u001b[0;34m(self, indexers, drop, missing_dims, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 2432\u001b[0m indexers \u001b[39m=\u001b[39m either_dict_or_kwargs(indexers, indexers_kwargs, \u001b[39m\"\u001b[39m\u001b[39misel\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 2433\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39many\u001b[39m(is_fancy_indexer(idx) \u001b[39mfor\u001b[39;00m idx \u001b[39min\u001b[39;00m indexers\u001b[39m.\u001b[39mvalues()):\n\u001b[0;32m-> 2434\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_isel_fancy(indexers, drop\u001b[39m=\u001b[39;49mdrop, missing_dims\u001b[39m=\u001b[39;49mmissing_dims)\n\u001b[1;32m 2436\u001b[0m \u001b[39m# Much faster algorithm for when all indexers are ints, slices, one-dimensional\u001b[39;00m\n\u001b[1;32m 2437\u001b[0m \u001b[39m# lists, or zero or one-dimensional np.ndarray's\u001b[39;00m\n\u001b[1;32m 2438\u001b[0m indexers \u001b[39m=\u001b[39m drop_dims_from_indexers(indexers, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdims, missing_dims)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/dataset.py:2490\u001b[0m, in \u001b[0;36mDataset._isel_fancy\u001b[0;34m(self, indexers, drop, missing_dims)\u001b[0m\n\u001b[1;32m 2486\u001b[0m var_indexers \u001b[39m=\u001b[39m {\n\u001b[1;32m 2487\u001b[0m k: v \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m valid_indexers\u001b[39m.\u001b[39mitems() \u001b[39mif\u001b[39;00m k \u001b[39min\u001b[39;00m var\u001b[39m.\u001b[39mdims\n\u001b[1;32m 2488\u001b[0m }\n\u001b[1;32m 2489\u001b[0m \u001b[39mif\u001b[39;00m var_indexers:\n\u001b[0;32m-> 2490\u001b[0m new_var \u001b[39m=\u001b[39m var\u001b[39m.\u001b[39;49misel(indexers\u001b[39m=\u001b[39;49mvar_indexers)\n\u001b[1;32m 2491\u001b[0m \u001b[39m# drop scalar coordinates\u001b[39;00m\n\u001b[1;32m 2492\u001b[0m \u001b[39m# https://github.com/pydata/xarray/issues/6554\u001b[39;00m\n\u001b[1;32m 2493\u001b[0m \u001b[39mif\u001b[39;00m name \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcoords \u001b[39mand\u001b[39;00m drop \u001b[39mand\u001b[39;00m new_var\u001b[39m.\u001b[39mndim \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m:\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:1336\u001b[0m, in \u001b[0;36mVariable.isel\u001b[0;34m(self, indexers, missing_dims, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 1333\u001b[0m indexers \u001b[39m=\u001b[39m drop_dims_from_indexers(indexers, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdims, missing_dims)\n\u001b[1;32m 1335\u001b[0m key \u001b[39m=\u001b[39m \u001b[39mtuple\u001b[39m(indexers\u001b[39m.\u001b[39mget(dim, \u001b[39mslice\u001b[39m(\u001b[39mNone\u001b[39;00m)) \u001b[39mfor\u001b[39;00m dim \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdims)\n\u001b[0;32m-> 1336\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m[key]\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:879\u001b[0m, in \u001b[0;36mVariable.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 866\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__getitem__\u001b[39m(\u001b[39mself\u001b[39m: T_Variable, key) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m T_Variable:\n\u001b[1;32m 867\u001b[0m \u001b[39m \u001b[39m\u001b[39m\"\"\"Return a new Variable object whose contents are consistent with\u001b[39;00m\n\u001b[1;32m 868\u001b[0m \u001b[39m getting the provided key from the underlying data.\u001b[39;00m\n\u001b[1;32m 869\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 877\u001b[0m \u001b[39m array `x.values` directly.\u001b[39;00m\n\u001b[1;32m 878\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 879\u001b[0m dims, indexer, new_order \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_broadcast_indexes(key)\n\u001b[1;32m 880\u001b[0m data \u001b[39m=\u001b[39m as_indexable(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_data)[indexer]\n\u001b[1;32m 881\u001b[0m \u001b[39mif\u001b[39;00m new_order:\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:725\u001b[0m, in \u001b[0;36mVariable._broadcast_indexes\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 723\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(k, Variable):\n\u001b[1;32m 724\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(k\u001b[39m.\u001b[39mdims) \u001b[39m>\u001b[39m \u001b[39m1\u001b[39m:\n\u001b[0;32m--> 725\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_broadcast_indexes_vectorized(key)\n\u001b[1;32m 726\u001b[0m dims\u001b[39m.\u001b[39mappend(k\u001b[39m.\u001b[39mdims[\u001b[39m0\u001b[39m])\n\u001b[1;32m 727\u001b[0m \u001b[39melif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39misinstance\u001b[39m(k, integer_types):\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/variable.py:849\u001b[0m, in \u001b[0;36mVariable._broadcast_indexes_vectorized\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 847\u001b[0m variables \u001b[39m=\u001b[39m _broadcast_compat_variables(\u001b[39m*\u001b[39mvariables)\n\u001b[1;32m 848\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m:\n\u001b[0;32m--> 849\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mIndexError\u001b[39;00m(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mDimensions of indexers mismatch: \u001b[39m\u001b[39m{\u001b[39;00mkey\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 851\u001b[0m out_key \u001b[39m=\u001b[39m [variable\u001b[39m.\u001b[39mdata \u001b[39mfor\u001b[39;00m variable \u001b[39min\u001b[39;00m variables]\n\u001b[1;32m 852\u001b[0m out_dims \u001b[39m=\u001b[39m \u001b[39mtuple\u001b[39m(out_dims_set)\n", + "\u001b[0;31mIndexError\u001b[0m: Dimensions of indexers mismatch: (slice(None, None, None), slice(None, None, None), \narray([[9, 9],\n [9, 9]]))" + ] + } + ], + "source": [ + "# So no CRS check against selection vectors currently:\n", + "DSll.sel(x=points, method='nearest')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integrations\n", + "If Python libraries store CRS information in predictable / standardized ways we can do neat things combining libraries. This will require some coordination..." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
stationgeometry
0siteAPOINT (-51.31734 -17.32280)
1siteBPOINT (-51.31745 -17.32290)
\n", + "
" + ], + "text/plain": [ + " station geometry\n", + "0 siteA POINT (-51.31734 -17.32280)\n", + "1 siteB POINT (-51.31745 -17.32290)" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gf = gpd.GeoDataFrame(dict(station=['siteA', 'siteB']), \n", + " geometry=gpd.points_from_xy(x=[-51.31734, -51.31745], y=[-17.3228, -17.3229]), \n", + " crs='EPSG:4326')\n", + "gf" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'station' (station: 2)>\n",
+       "'siteA' 'siteB'\n",
+       "Coordinates:\n",
+       "  * station  (station) object POINT (-51.31734 -17.3228) POINT (-51.31745 -17...\n",
+       "Indexes:\n",
+       "    station  GeometryIndex (crs=EPSG:4326)
" + ], + "text/plain": [ + "\n", + "'siteA' 'siteB'\n", + "Coordinates:\n", + " * station (station) object POINT (-51.31734 -17.3228) POINT (-51.31745 -17...\n", + "Indexes:\n", + " station GeometryIndex (crs=EPSG:4326)" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Turn points into an Xarray GeometryIndex with xvec (https://xvec.readthedocs.io)\n", + "da = xr.DataArray(gf.station, coords=[gf.geometry], dims=[\"station\"])\n", + "da = da.xvec.set_geom_indexes(\"station\", crs=gf.crs)\n", + "da" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
+       "    spatial_ref  int64 ...\n",
+       "Data variables:\n",
+       "    blue         (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n",
+       "    green        (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, x: 10, y: 10)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " * x (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n", + " * y (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n", + " spatial_ref int64 ...\n", + "Data variables:\n", + " blue (time, y, x) float64 6.611 5.581 0.3996 ... 3.491 5.056 3.368\n", + " green (time, y, x) float64 7.921 66.15 30.1 ... 21.76 27.29 18.41\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DS" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " {'x': \n", + "-51.32 -51.32\n", + "Dimensions without coordinates: geometry, 'y': \n", + "-17.32 -17.32\n", + "Dimensions without coordinates: geometry} {'method': 'nearest', 'tolerance': None}\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, geometry: 2)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * geometry     (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n",
+       "Data variables:\n",
+       "    blue         (time, geometry) float64 5.48 4.536 0.5377 6.746\n",
+       "    green        (time, geometry) float64 57.79 23.8 20.33 6.438\n",
+       "Indexes:\n",
+       "    geometry     GeometryIndex (crs=EPSG:4326)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, geometry: 2)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * geometry (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n", + "Data variables:\n", + " blue (time, geometry) float64 5.48 4.536 0.5377 6.746\n", + " green (time, geometry) float64 57.79 23.8 20.33 6.438\n", + "Indexes:\n", + " geometry GeometryIndex (crs=EPSG:4326)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NOTE: modify xvec.extract_points to automatically get x_dim, y_dim? \n", + "DSll.xvec.extract_points(da.station, DSll.rio.x_dim, DSll.rio.y_dim)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " {'x': \n", + "-51.32 -51.32\n", + "Dimensions without coordinates: geometry, 'y': \n", + "-17.32 -17.32\n", + "Dimensions without coordinates: geometry} {'method': 'nearest', 'tolerance': None}\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, geometry: 2)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * geometry     (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n",
+       "Data variables:\n",
+       "    blue         (time, geometry) float64 5.48 4.536 0.5377 6.746\n",
+       "    green        (time, geometry) float64 57.79 23.8 20.33 6.438\n",
+       "Indexes:\n",
+       "    geometry     GeometryIndex (crs=EPSG:4326)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, geometry: 2)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * geometry (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n", + "Data variables:\n", + " blue (time, geometry) float64 5.48 4.536 0.5377 6.746\n", + " green (time, geometry) float64 57.79 23.8 20.33 6.438\n", + "Indexes:\n", + " geometry GeometryIndex (crs=EPSG:4326)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NOTE: could modify extract_points() to work just with `gf` as input?\n", + "DSll.xvec.extract_points(gf.geometry, x_coords='x', y_coords='y')" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " {'x': \n", + "-51.32 -51.32\n", + "Dimensions without coordinates: geometry, 'y': \n", + "-17.32 -17.32\n", + "Dimensions without coordinates: geometry} {'method': 'nearest', 'tolerance': None}\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 2, geometry: 2)\n",
+       "Coordinates:\n",
+       "  * time         (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n",
+       "    spatial_ref  int64 0\n",
+       "  * geometry     (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n",
+       "Data variables:\n",
+       "    blue         (time, geometry) float64 0.8441 0.8441 1.844 1.844\n",
+       "    green        (time, geometry) float64 37.09 37.09 44.67 44.67\n",
+       "Indexes:\n",
+       "    geometry     GeometryIndex (crs=EPSG:4326)\n",
+       "Attributes:\n",
+       "    coordinates:  spatial_ref
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 2, geometry: 2)\n", + "Coordinates:\n", + " * time (time) object 2016-12-19 10:27:29.687763 2016-12-29 12:52:42...\n", + " spatial_ref int64 0\n", + " * geometry (geometry) object POINT (-51.31734 -17.3228) POINT (-51.3174...\n", + "Data variables:\n", + " blue (time, geometry) float64 0.8441 0.8441 1.844 1.844\n", + " green (time, geometry) float64 37.09 37.09 44.67 44.67\n", + "Indexes:\n", + " geometry GeometryIndex (crs=EPSG:4326)\n", + "Attributes:\n", + " coordinates: spatial_ref" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NOTE: want this to actually fail with a CRS check (points in raster in UTM, but points in lon/lat)\n", + "DSutm.xvec.extract_points(da.station, x_coords='x', y_coords='y')" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "CRS mismatch between the CRS of index geometries and the CRS of input geometries.\nIndex CRS: EPSG:32722\nInput CRS: EPSG:4326\n", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[35], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m xr\u001b[39m.\u001b[39;49malign(DSutm, DSll)\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:787\u001b[0m, in \u001b[0;36malign\u001b[0;34m(join, copy, indexes, exclude, fill_value, *objects)\u001b[0m\n\u001b[1;32m 591\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 592\u001b[0m \u001b[39mGiven any number of Dataset and/or DataArray objects, returns new\u001b[39;00m\n\u001b[1;32m 593\u001b[0m \u001b[39mobjects with aligned indexes and dimension sizes.\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 777\u001b[0m \n\u001b[1;32m 778\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m 779\u001b[0m aligner \u001b[39m=\u001b[39m Aligner(\n\u001b[1;32m 780\u001b[0m objects,\n\u001b[1;32m 781\u001b[0m join\u001b[39m=\u001b[39mjoin,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 785\u001b[0m fill_value\u001b[39m=\u001b[39mfill_value,\n\u001b[1;32m 786\u001b[0m )\n\u001b[0;32m--> 787\u001b[0m aligner\u001b[39m.\u001b[39;49malign()\n\u001b[1;32m 788\u001b[0m \u001b[39mreturn\u001b[39;00m aligner\u001b[39m.\u001b[39mresults\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:572\u001b[0m, in \u001b[0;36mAligner.align\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 570\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfind_matching_unindexed_dims()\n\u001b[1;32m 571\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39massert_no_index_conflict()\n\u001b[0;32m--> 572\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49malign_indexes()\n\u001b[1;32m 573\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39massert_unindexed_dim_sizes_equal()\n\u001b[1;32m 575\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mjoin \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39moverride\u001b[39m\u001b[39m\"\u001b[39m:\n", + "File \u001b[0;32m~/miniconda3/envs/rioxarray/lib/python3.11/site-packages/xarray/core/alignment.py:424\u001b[0m, in \u001b[0;36mAligner.align_indexes\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 417\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 418\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mcannot align objects with join=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mexact\u001b[39m\u001b[39m'\u001b[39m\u001b[39m where \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 419\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mindex/labels/sizes are not equal along \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 420\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mthese coordinates (dimensions): \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 421\u001b[0m \u001b[39m+\u001b[39m \u001b[39m\"\u001b[39m\u001b[39m, \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mjoin(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mname\u001b[39m!r}\u001b[39;00m\u001b[39m \u001b[39m\u001b[39m{\u001b[39;00mdims\u001b[39m!r}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m name, dims \u001b[39min\u001b[39;00m key[\u001b[39m0\u001b[39m])\n\u001b[1;32m 422\u001b[0m )\n\u001b[1;32m 423\u001b[0m joiner \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_get_index_joiner(index_cls)\n\u001b[0;32m--> 424\u001b[0m joined_index \u001b[39m=\u001b[39m joiner(matching_indexes)\n\u001b[1;32m 425\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mjoin \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39mleft\u001b[39m\u001b[39m\"\u001b[39m:\n\u001b[1;32m 426\u001b[0m joined_index_vars \u001b[39m=\u001b[39m matching_index_vars[\u001b[39m0\u001b[39m]\n", + "File \u001b[0;32m~/GitHub/rioxarray/rioxarray/_crs_index.py:150\u001b[0m, in \u001b[0;36mCRSIndex.join\u001b[0;34m(self, other, how)\u001b[0m\n\u001b[1;32m 148\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mjoin\u001b[39m(\u001b[39mself\u001b[39m, other, how\u001b[39m=\u001b[39m\u001b[39m\"\u001b[39m\u001b[39minner\u001b[39m\u001b[39m\"\u001b[39m):\n\u001b[1;32m 149\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_check_crs(other\u001b[39m.\u001b[39mcrs, allow_none\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m):\n\u001b[0;32m--> 150\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_crs_mismatch_raise(other\u001b[39m.\u001b[39;49mcrs)\n\u001b[1;32m 152\u001b[0m new_indexes \u001b[39m=\u001b[39m {\n\u001b[1;32m 153\u001b[0m k: v\u001b[39m.\u001b[39mjoin(other\u001b[39m.\u001b[39m_indexes[k], how\u001b[39m=\u001b[39mhow) \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_indexes\u001b[39m.\u001b[39mitems()\n\u001b[1;32m 154\u001b[0m }\n\u001b[1;32m 156\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mtype\u001b[39m(\u001b[39mself\u001b[39m)(new_indexes, crs\u001b[39m=\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcrs)\n", + "File \u001b[0;32m~/GitHub/rioxarray/rioxarray/_crs_index.py:137\u001b[0m, in \u001b[0;36mCRSIndex._crs_mismatch_raise\u001b[0;34m(self, other_crs, stacklevel)\u001b[0m\n\u001b[1;32m 129\u001b[0m msg \u001b[39m=\u001b[39m (\n\u001b[1;32m 130\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mCRS mismatch between the CRS of index geometries \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 131\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mand the CRS of input geometries.\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 132\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mIndex CRS: \u001b[39m\u001b[39m{\u001b[39;00msrs\u001b[39m}\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 133\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mInput CRS: \u001b[39m\u001b[39m{\u001b[39;00mother_srs\u001b[39m}\u001b[39;00m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 134\u001b[0m )\n\u001b[1;32m 136\u001b[0m \u001b[39mif\u001b[39;00m get_option(CRS_INDEX_ERROR):\n\u001b[0;32m--> 137\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(msg)\n\u001b[1;32m 138\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m 139\u001b[0m warnings\u001b[39m.\u001b[39mwarn(msg, \u001b[39mUserWarning\u001b[39;00m, stacklevel\u001b[39m=\u001b[39mstacklevel)\n", + "\u001b[0;31mValueError\u001b[0m: CRS mismatch between the CRS of index geometries and the CRS of input geometries.\nIndex CRS: EPSG:32722\nInput CRS: EPSG:4326\n" + ] + } + ], + "source": [ + "xr.align(DSutm, DSll)" + ] + } + ], + "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.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/rioxarray/_crs_index.py b/rioxarray/_crs_index.py new file mode 100644 index 00000000..695aaeb0 --- /dev/null +++ b/rioxarray/_crs_index.py @@ -0,0 +1,165 @@ +""" +Experimental implementation of xarray custom index + +https://github.com/corteva/rioxarray/issues/588 +https://github.com/xarray-contrib/xvec/blob/main/xvec/index.py +""" +from __future__ import annotations + +import warnings +from collections.abc import Hashable, Mapping +from typing import Any + +import pandas as pd +from rasterio.crs import CRS +from xarray import Variable, get_options +from xarray.core.indexes import Index, PandasIndex, get_indexer_nd +from xarray.core.indexing import IndexSelResult, merge_sel_results + +from rioxarray._options import CRS_INDEX_ERROR, get_option + + +# From XVEC +def _format_crs(crs: CRS, max_width: int = 50) -> str: + if crs is not None: + srs = crs.to_string() + else: + srs = "None" + + return srs if len(srs) <= max_width else " ".join([srs[:max_width], "..."]) + + +class CRSIndex(Index): + """Coordinate-Refernce System aware, Xarray compatible 2D index for rasters. + A 'MetaIndex' that adds CRS information to 2 related spatial dimensions (x,y) + The CRS defines coordinate system wherein the (x,y) points reside + """ + + def __init__(self, indexes: [PandasIndex, PandasIndex], crs: CRS): + self._indexes = indexes + self._crs = crs + + @property + def crs(self) -> CRS | None: + """Returns the coordinate reference system of the index as a + :class:`rasterio.crs.CRS` object. + """ + return self._crs + + @classmethod + def from_variables( + cls, + variables: Mapping[Any, Variable], + *, + options: Mapping[str, Any], + ): + + xy_indexes = { + k: PandasIndex.from_variables({k: v}, options=options) + for k, v in variables.items() + if k in ["x", "y"] + } + + return cls(xy_indexes, crs=options.get("crs")) + + def create_variables(self, variables=None): + idx_variables = {} + + for index in self._indexes.values(): + idx_variables.update(index.create_variables(variables)) + + return idx_variables + + def sel(self, labels, **kwargs) -> IndexSelResult: + print(self, labels, kwargs) + results = [] + for k, index in self._indexes.items(): + if k in labels: + results.append(index.sel({k: labels[k]}, **kwargs)) + + result = merge_sel_results(results) + return result + + def isel(self, indexers): + # Do a CRS check + # print(indexers) + results = {} + for k, index in self._indexes.items(): + if k in indexers: + result_idx = index.isel(indexers) + if result_idx is not None: + results[k] = result_idx + else: + results[k] = index + + return type(self)(results, crs=self.crs) + + def to_pandas_index(self) -> pd.Index: + return self._index.index + + def _repr_inline_(self, max_width: int): + # TODO: remove when fixed in XArray, Open Issue? + if max_width is None: + max_width = get_options()["display_width"] + srs = _format_crs(self.crs, max_width=max_width) + return f"{self.__class__.__name__} (crs={srs})" + + def _check_crs(self, other_crs: CRS | None, allow_none: bool = False) -> bool: + """Check if the index's projection is the same than the given one. + If allow_none is True, empty CRS is treated as the same. + """ + # TODO: ignore_axis_order https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.equals + if allow_none: + if self.crs is None or other_crs is None: + return True + if not self.crs == other_crs: + return False + return True + + def _crs_mismatch_raise(self, other_crs: CRS | None, stacklevel: int = 3): + """Raise a CRS mismatch error or warning with the information + on the assigned CRS. + """ + srs = _format_crs(self.crs, max_width=50) + other_srs = _format_crs(other_crs, max_width=50) + + msg = ( + "CRS mismatch between the CRS of index geometries " + "and the CRS of input geometries.\n" + f"Index CRS: {srs}\n" + f"Input CRS: {other_srs}\n" + ) + + if get_option(CRS_INDEX_ERROR): + raise ValueError(msg) + else: + warnings.warn(msg, UserWarning, stacklevel=stacklevel) + + def equals(self, other: Index) -> bool: + if not isinstance(other, CRSIndex): + return False + if not self._check_crs(other.crs, allow_none=True): + return False + + def join(self, other, how="inner"): + if not self._check_crs(other.crs, allow_none=True): + self._crs_mismatch_raise(other.crs) + + new_indexes = { + k: v.join(other._indexes[k], how=how) for k, v in self._indexes.items() + } + + return type(self)(new_indexes, crs=self.crs) + + def reindex_like(self, other, method=None, tolerance=None) -> dict[Hashable, Any]: + if not self._check_crs(other.crs, allow_none=True): + self._crs_mismatch_raise(other.crs) + + new_index = { + k: get_indexer_nd( + self._indexes[k].index, other._indexes[k].index, method, tolerance + ) + for k in self._indexes.keys() + } + + return new_index diff --git a/rioxarray/_io.py b/rioxarray/_io.py index f874c8b4..eaddc21f 100644 --- a/rioxarray/_io.py +++ b/rioxarray/_io.py @@ -793,6 +793,7 @@ def _ds_close(): dataset = dataset.pop() else: dataset = Dataset(attrs=global_tags) + return dataset @@ -1256,10 +1257,11 @@ def open_rasterio( result.rio.write_transform(riods.transform, inplace=True) rio_crs = riods.crs or result.rio.crs if rio_crs: - result.rio.write_crs(rio_crs, inplace=True) + # unable to get inplace=True to work with open_dataarray or open_dataset engine='rasterio' + # result.rio.write_crs(rio_crs, inplace=True) + result = result.rio.write_crs(rio_crs) if has_gcps: result.rio.write_gcps(*riods.gcps, inplace=True) - if chunks is not None: result = _prepare_dask(result, riods, filename, chunks) diff --git a/rioxarray/_options.py b/rioxarray/_options.py index 1dc55ffa..f4133fa1 100644 --- a/rioxarray/_options.py +++ b/rioxarray/_options.py @@ -10,10 +10,14 @@ EXPORT_GRID_MAPPING = "export_grid_mapping" SKIP_MISSING_SPATIAL_DIMS = "skip_missing_spatial_dims" +USE_CRS_INDEX = "use_crs_index" +CRS_INDEX_ERROR = "crs_index_error" OPTIONS = { EXPORT_GRID_MAPPING: True, SKIP_MISSING_SPATIAL_DIMS: False, + USE_CRS_INDEX: False, + CRS_INDEX_ERROR: True, } OPTION_NAMES = set(OPTIONS) @@ -60,6 +64,10 @@ class set_options: # pylint: disable=invalid-name If True, it will not perform spatial operations on variables within a :class:`xarray.Dataset` if the spatial dimensions are not found. + use_crs_index: bool, default=False + If True, will store CRS information as an xarray custom Index + crs_index_error: bool, default=True + If True raise an exception when CRSIndex are different, otherwise raise a warning. Usage as a context manager:: diff --git a/rioxarray/merge.py b/rioxarray/merge.py index 13d8899b..0a2dd877 100644 --- a/rioxarray/merge.py +++ b/rioxarray/merge.py @@ -204,7 +204,8 @@ def merge_arrays( xda.rio.write_nodata( nodata if nodata is not None else representative_array.rio.nodata, inplace=True ) - xda.rio.write_crs(representative_array.rio.crs, inplace=True) + # xda.rio.write_crs(representative_array.rio.crs, inplace=True) + xda = xda.rio.write_crs(representative_array.rio.crs) xda.rio.write_transform(merged_transform, inplace=True) return xda @@ -279,5 +280,6 @@ def merge_datasets( ), attrs=representative_ds.attrs, ) - xds.rio.write_crs(merged_data[data_var].rio.crs, inplace=True) + # xds.rio.write_crs(merged_data[data_var].rio.crs, inplace=True) + xds = xds.rio.write_crs(merged_data[data_var].rio.crs) return xds diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py index 30c7c1bd..97e31248 100644 --- a/rioxarray/rioxarray.py +++ b/rioxarray/rioxarray.py @@ -17,7 +17,8 @@ from rasterio.control import GroundControlPoint from rasterio.crs import CRS -from rioxarray._options import EXPORT_GRID_MAPPING, get_option +from rioxarray._crs_index import CRSIndex +from rioxarray._options import EXPORT_GRID_MAPPING, USE_CRS_INDEX, get_option from rioxarray.crs import crs_from_user_input from rioxarray.exceptions import ( DimensionError, @@ -203,7 +204,8 @@ def _get_spatial_dims( def _has_spatial_dims( - obj: Union[xarray.Dataset, xarray.DataArray], var: Union[Any, Hashable] + obj: Union[xarray.Dataset, xarray.DataArray], + var: Union[Any, Hashable], ) -> bool: """ Check to see if the variable in the Dataset has spatial dimensions @@ -511,10 +513,52 @@ def write_crs( # remove old crs if exists data_obj.attrs.pop("crs", None) - return data_obj.rio.write_grid_mapping( + if get_option(USE_CRS_INDEX): + no_existing_index = CRSIndex not in [ + type(i) for i in data_obj.xindexes.values() + ] + + try: + has_spatial_dims = hasattr(data_obj.rio, "x_dim") and hasattr( + data_obj.rio, "y_dim" + ) + except MissingSpatialDimensionError: + warnings.warn( + "MissingSpatialDimensionError: x or y dimensions required for CRSIndex. See 'rio.set_spatial_dims()'" + ) + has_spatial_dims = False + pass + + if no_existing_index and has_spatial_dims: + # when loading netcdf with multiple variables, not all coords parsed + has_coords = ( + data_obj.rio.x_dim and data_obj.rio.y_dim in data_obj.coords + ) + if has_coords: + xdim, ydim = data_obj.rio.x_dim, data_obj.rio.y_dim + + # Must use either old .indexes or new .xindexes + # avoid ValueError: those coordinates already have an index: {'y', 'x'} + data_obj = data_obj.drop_indexes([xdim, ydim]) + + data_obj = data_obj.set_xindex( + ( + xdim, + ydim, + ), + CRSIndex, + crs=data_obj.rio.crs, + ) + + # Fix to appease some tests: test_nonstandard_dims* + data_obj.rio.set_spatial_dims(xdim, ydim, inplace=True) + + result = data_obj.rio.write_grid_mapping( grid_mapping_name=grid_mapping_name, inplace=True ) + return result + def estimate_utm_crs(self, datum_name: str = "WGS 84") -> rasterio.crs.CRS: """Returns the estimated UTM CRS based on the bounds of the dataset. diff --git a/rioxarray/xarray_plugin.py b/rioxarray/xarray_plugin.py index 01f6c4d1..e6936c1a 100644 --- a/rioxarray/xarray_plugin.py +++ b/rioxarray/xarray_plugin.py @@ -82,6 +82,7 @@ def open_dataset( ) if drop_variables is not None: rds = rds.drop_vars(drop_variables) + return rds def guess_can_open(self, filename_or_obj): # pylint: disable=arguments-renamed diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index 76e9b34b..b3d8b65d 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/test/integration/test_integration_rioxarray.py @@ -41,6 +41,9 @@ open_rasterio_engine, ) +# Set global option to test CRSIndex across all tests +# rioxarray.set_options(use_crs_index=True) + try: SCIPY_LT_17 = version.parse(importlib.metadata.version("scipy")) < version.parse( "1.7.0" @@ -2128,6 +2131,20 @@ def test_crs_is_removed(): assert "crs" not in test_ds.attrs +def test_write_crs_index(): + test_da = xarray.DataArray( + numpy.zeros((5, 5)), + dims=("y", "x"), + coords={"y": numpy.arange(1, 6), "x": numpy.arange(2, 7)}, + ) + with rioxarray.set_options(use_crs_index=True): + test_da = test_da.rio.write_crs(4326) + assert hasattr(test_da, "indexes") is False + assert hasattr(test_da, "xindexes") is True + assert isinstance(test_da.xindexes["x"], rioxarray._crs_index.CRSIndex) + assert test_da.xindexes["y"].crs == test_da.rio.crs + + def test_write_crs_cf(): test_da = xarray.DataArray(1) test_da = test_da.rio.write_crs(4326)