diff --git a/.gitignore b/.gitignore index 3be58bb..fb812f6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] +*.tif # C extensions *.so diff --git a/README.md b/README.md index 00ec601..fa8c4d4 100644 --- a/README.md +++ b/README.md @@ -41,3 +41,15 @@ pip install rasterix hatch env run --env test.py3.13 run-pytest # Run the tests without coverage reports hatch env run --env test.py3.13 run-coverage-html # Run the tests with an html coverage report ``` + +#### Using hatch + +Start a shell +``` +hatch env run --env test.py3.13 ipython +``` + +Expose kernel for notebooks in JupyterLab +``` +hatch env run --env test.py3.13 "python -m ipykernel install --user --name=rasterix" +``` diff --git a/notebooks/xproj_integration.ipynb b/notebooks/xproj_integration.ipynb new file mode 100644 index 0000000..abe00c9 --- /dev/null +++ b/notebooks/xproj_integration.ipynb @@ -0,0 +1,623 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4e6cc2f1-ea8f-4a2a-b5c5-38c54ab0ff42", + "metadata": {}, + "source": [ + "# Coordinate Reference System (CRS) Management\n", + "\n", + "Rasterix uses [XProj](https://xproj.readthedocs.io/en/latest/) to store and manage the CRS of a dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "b9faca03-d650-4e60-9d4f-3205e432e855", + "metadata": {}, + "outputs": [], + "source": [ + "import rasterix\n", + "import xarray as xr\n", + "import xproj\n", + "\n", + "# GDAL settings for remote data\n", + "import os\n", + "os.environ['VSICURL_PC_URL_SIGNING']='YES'\n", + "os.environ['GDAL_DISABLE_READDIR_ON_OPEN'] = 'EMPTY_DIR'\n", + "\n", + "xr.set_options(display_expand_indexes=True);" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "1674789d-3646-4d16-bfb8-2ec73eeed00a", + "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: 332, x: 316)> Size: 420kB\n",
+       "[104912 values with dtype=float32]\n",
+       "Coordinates:\n",
+       "  * band         (band) int64 8B 1\n",
+       "    spatial_ref  int64 8B ...\n",
+       "  * x            (x) float64 3kB -3.938e+06 -3.912e+06 ... 3.912e+06 3.938e+06\n",
+       "  * y            (y) float64 3kB 4.338e+06 4.312e+06 ... -3.912e+06 -3.938e+06\n",
+       "  * crs          int64 8B 0\n",
+       "Indexes:\n",
+       "  ┌ x        RasterIndex\n",
+       "  └ y\n",
+       "    crs      CRSIndex (crs=EPSG:3412)\n",
+       "Attributes:\n",
+       "    AREA_OR_POINT:  Area
" + ], + "text/plain": [ + " Size: 420kB\n", + "[104912 values with dtype=float32]\n", + "Coordinates:\n", + " * band (band) int64 8B 1\n", + " spatial_ref int64 8B ...\n", + " * x (x) float64 3kB -3.938e+06 -3.912e+06 ... 3.912e+06 3.938e+06\n", + " * y (y) float64 3kB 4.338e+06 4.312e+06 ... -3.912e+06 -3.938e+06\n", + " * crs int64 8B 0\n", + "Indexes:\n", + " ┌ x RasterIndex\n", + " └ y\n", + " crs CRSIndex (crs=EPSG:3412)\n", + "Attributes:\n", + " AREA_OR_POINT: Area" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "url = \"https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif\"\n", + "da = xr.open_dataarray(url, engine=\"rasterio\")\n", + "\n", + "# How to get this to automatically pick up CRS that is detected by rioxarray?\n", + "# Can this be done automatically by assign_index?\n", + "# Should assign_index happen automatically? or would this be done with engine='rasterix')\n", + "da = rasterix.assign_index(da).proj.assign_crs(crs=da.rio.crs)\n", + "da" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "ebbc212e-8e14-45fb-8789-b81dcb0e65d6", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/jovyan/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/IPython/core/interactiveshell.py:3153: SerializationWarning: saving variable None with floating point data as an integer dtype without any _FillValue to use for NaNs\n", + " result = runner(coro)\n" + ] + } + ], + "source": [ + "# Make sure we can save this back to a geotiff \n", + "# Hmmm, this warning isn't clear to me\n", + "da.rio.to_raster('/tmp/S_20240101_concentration_v3.0.tif')" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "81525346-2ef2-4294-9471-ba0c2adb0e36", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "File ‘S_20240101_concentration_v3.0.tif’ already there; not retrieving.\n", + "\n", + "-rw-r--r-- 1 jovyan jovyan 606187 Jan 2 2024 S_20240101_concentration_v3.0.tif\n" + ] + } + ], + "source": [ + "!wget -nc https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif\n", + "!ls -l S_20240101_concentration_v3.0.tif" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "5ca2dac1-943e-4cf2-92aa-0eda31a57fec", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-rw-r--r-- 1 jovyan jovyan 213021 May 13 13:51 /tmp/S_20240101_concentration_v3.0.tif\n" + ] + } + ], + "source": [ + "# Are the tiffs the same?\n", + "!ls -l /tmp/S_20240101_concentration_v3.0.tif" + ] + }, + { + "cell_type": "markdown", + "id": "ff4b111f-06ac-4e91-830b-4014a4f042ea", + "metadata": {}, + "source": [ + "## Two rasters with same GeoTransform but different CRS\n", + "\n", + "Here is a case for needing to account for both affine and CRS. [Sentinel-2](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) is delivered on the MGRS grid, which divides UTM Zones into smaller units. Adjacent zones have the same affines, and therefore appear to have the same coordinates. But the CRS differentiates them.\n", + "\n", + "We want Xarray to raise an error in this case before attempting any computations." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "176c3837-05de-4d18-9337-e22445b41ab0", + "metadata": {}, + "outputs": [ + { + "ename": "MergeError", + "evalue": "conflicting values/indexes on objects to be combined for coordinate 'crs'\nfirst index: CRSIndex\n\nName: WGS 84 / UTM zone 10N\nAxis Info [cartesian]:\n- [east]: Easting (metre)\n- [north]: Northing (metre)\nArea of Use:\n- undefined\nCoordinate Operation:\n- name: UTM zone 10N\n- method: Transverse Mercator\nDatum: World Geodetic System 1984\n- Ellipsoid: WGS 84\n- Prime Meridian: Greenwich\n\nsecond index: CRSIndex\n\nName: WGS 84 / UTM zone 11N\nAxis Info [cartesian]:\n- [east]: Easting (metre)\n- [north]: Northing (metre)\nArea of Use:\n- undefined\nCoordinate Operation:\n- name: UTM zone 11N\n- method: Transverse Mercator\nDatum: World Geodetic System 1984\n- Ellipsoid: WGS 84\n- Prime Meridian: Greenwich\n\nfirst variable: Size: 8B\narray(0)\nsecond variable: Size: 8B\narray(0)\n", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mMergeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[22]\u001b[39m\u001b[32m, line 17\u001b[39m\n\u001b[32m 14\u001b[39m url = \u001b[33m'\u001b[39m\u001b[33mhttps://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/11/T/MM/2025/04/13/S2B_MSIL2A_20250413T184919_N0511_R113_T11TMM_20250413T224733.SAFE/GRANULE/L2A_T11TMM_A042324_20250413T185100/IMG_DATA/R10m/T11TMM_20250413T184919_B04_10m.tif\u001b[39m\u001b[33m'\u001b[39m\n\u001b[32m 15\u001b[39m da2 = load_raster(url)\n\u001b[32m---> \u001b[39m\u001b[32m17\u001b[39m \u001b[43mda1\u001b[49m\u001b[43m \u001b[49m\u001b[43m+\u001b[49m\u001b[43m \u001b[49m\u001b[43mda2\u001b[49m \n", + "\u001b[36mFile \u001b[39m\u001b[32m~/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/xarray/core/_typed_ops.py:528\u001b[39m, in \u001b[36mDataArrayOpsMixin.__add__\u001b[39m\u001b[34m(self, other)\u001b[39m\n\u001b[32m 527\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34m__add__\u001b[39m(\u001b[38;5;28mself\u001b[39m, other: DaCompatible) -> Self | Dataset | DataTree:\n\u001b[32m--> \u001b[39m\u001b[32m528\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43m_binary_op\u001b[49m\u001b[43m(\u001b[49m\u001b[43mother\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43moperator\u001b[49m\u001b[43m.\u001b[49m\u001b[43madd\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/xarray/core/dataarray.py:4853\u001b[39m, in \u001b[36mDataArray._binary_op\u001b[39m\u001b[34m(self, other, f, reflexive)\u001b[39m\n\u001b[32m 4846\u001b[39m other_coords = \u001b[38;5;28mgetattr\u001b[39m(other, \u001b[33m\"\u001b[39m\u001b[33mcoords\u001b[39m\u001b[33m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[32m 4848\u001b[39m variable = (\n\u001b[32m 4849\u001b[39m f(\u001b[38;5;28mself\u001b[39m.variable, other_variable_or_arraylike)\n\u001b[32m 4850\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m reflexive\n\u001b[32m 4851\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m f(other_variable_or_arraylike, \u001b[38;5;28mself\u001b[39m.variable)\n\u001b[32m 4852\u001b[39m )\n\u001b[32m-> \u001b[39m\u001b[32m4853\u001b[39m coords, indexes = \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m.\u001b[49m\u001b[43mcoords\u001b[49m\u001b[43m.\u001b[49m\u001b[43m_merge_raw\u001b[49m\u001b[43m(\u001b[49m\u001b[43mother_coords\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreflexive\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 4854\u001b[39m name = result_name([\u001b[38;5;28mself\u001b[39m, other])\n\u001b[32m 4856\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m._replace(variable, coords, name, indexes=indexes)\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/xarray/core/coordinates.py:505\u001b[39m, in \u001b[36mCoordinates._merge_raw\u001b[39m\u001b[34m(self, other, reflexive)\u001b[39m\n\u001b[32m 503\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[32m 504\u001b[39m coord_list = [\u001b[38;5;28mself\u001b[39m, other] \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m reflexive \u001b[38;5;28;01melse\u001b[39;00m [other, \u001b[38;5;28mself\u001b[39m]\n\u001b[32m--> \u001b[39m\u001b[32m505\u001b[39m variables, indexes = \u001b[43mmerge_coordinates_without_align\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcoord_list\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 506\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m variables, indexes\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/xarray/structure/merge.py:418\u001b[39m, in \u001b[36mmerge_coordinates_without_align\u001b[39m\u001b[34m(objects, prioritized, exclude_dims, combine_attrs)\u001b[39m\n\u001b[32m 414\u001b[39m filtered = collected\n\u001b[32m 416\u001b[39m \u001b[38;5;66;03m# TODO: indexes should probably be filtered in collected elements\u001b[39;00m\n\u001b[32m 417\u001b[39m \u001b[38;5;66;03m# before merging them\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m418\u001b[39m merged_coords, merged_indexes = \u001b[43mmerge_collected\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 419\u001b[39m \u001b[43m \u001b[49m\u001b[43mfiltered\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mprioritized\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcombine_attrs\u001b[49m\u001b[43m=\u001b[49m\u001b[43mcombine_attrs\u001b[49m\n\u001b[32m 420\u001b[39m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 421\u001b[39m merged_indexes = filter_indexes_from_coords(merged_indexes, \u001b[38;5;28mset\u001b[39m(merged_coords))\n\u001b[32m 423\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m merged_coords, merged_indexes\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/.local/share/hatch/env/virtual/rasterix/UNYVAJKM/test.py3.13/lib/python3.13/site-packages/xarray/structure/merge.py:274\u001b[39m, in \u001b[36mmerge_collected\u001b[39m\u001b[34m(grouped, prioritized, compat, combine_attrs, equals)\u001b[39m\n\u001b[32m 270\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m other_var, other_index \u001b[38;5;129;01min\u001b[39;00m indexed_elements[\u001b[32m1\u001b[39m:]:\n\u001b[32m 271\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m indexes_equal(\n\u001b[32m 272\u001b[39m index, other_index, variable, other_var, index_cmp_cache\n\u001b[32m 273\u001b[39m ):\n\u001b[32m--> \u001b[39m\u001b[32m274\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m MergeError(\n\u001b[32m 275\u001b[39m \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mconflicting values/indexes on objects to be combined for coordinate \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mname\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m 276\u001b[39m \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mfirst index: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mindex\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33msecond index: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mother_index\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m 277\u001b[39m \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33mfirst variable: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mvariable\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33msecond variable: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mother_var\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m 278\u001b[39m )\n\u001b[32m 279\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m compat == \u001b[33m\"\u001b[39m\u001b[33midentical\u001b[39m\u001b[33m\"\u001b[39m:\n\u001b[32m 280\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m other_variable, _ \u001b[38;5;129;01min\u001b[39;00m indexed_elements[\u001b[32m1\u001b[39m:]:\n", + "\u001b[31mMergeError\u001b[39m: conflicting values/indexes on objects to be combined for coordinate 'crs'\nfirst index: CRSIndex\n\nName: WGS 84 / UTM zone 10N\nAxis Info [cartesian]:\n- [east]: Easting (metre)\n- [north]: Northing (metre)\nArea of Use:\n- undefined\nCoordinate Operation:\n- name: UTM zone 10N\n- method: Transverse Mercator\nDatum: World Geodetic System 1984\n- Ellipsoid: WGS 84\n- Prime Meridian: Greenwich\n\nsecond index: CRSIndex\n\nName: WGS 84 / UTM zone 11N\nAxis Info [cartesian]:\n- [east]: Easting (metre)\n- [north]: Northing (metre)\nArea of Use:\n- undefined\nCoordinate Operation:\n- name: UTM zone 11N\n- method: Transverse Mercator\nDatum: World Geodetic System 1984\n- Ellipsoid: WGS 84\n- Prime Meridian: Greenwich\n\nfirst variable: Size: 8B\narray(0)\nsecond variable: Size: 8B\narray(0)\n" + ] + } + ], + "source": [ + "def load_raster(href, overview_level=3, masked=True):\n", + " da = xr.open_dataarray(href, \n", + " engine='rasterio', \n", + " masked=masked,\n", + " open_kwargs=dict(overview_level=overview_level),\n", + " chunks='auto').squeeze()\n", + " \n", + " return rasterix.assign_index(da).proj.assign_crs(crs=da.rio.crs)\n", + "\n", + "\n", + "url = 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/10/T/DS/2025/04/14/S2C_MSIL2A_20250414T190931_N0511_R056_T10TDS_20250415T001314.SAFE/GRANULE/L2A_T10TDS_A003172_20250414T191846/IMG_DATA/R10m/T10TDS_20250414T190931_B04_10m.tif'\n", + "da1 = load_raster(url)\n", + "\n", + "url = 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/11/T/MM/2025/04/13/S2B_MSIL2A_20250413T184919_N0511_R113_T11TMM_20250413T224733.SAFE/GRANULE/L2A_T11TMM_A042324_20250413T185100/IMG_DATA/R10m/T11TMM_20250413T184919_B04_10m.tif'\n", + "da2 = load_raster(url)\n", + " \n", + "da1 + da2 " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7e25c4d-d56d-4b2e-8e9c-343c72e91cb8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "rasterix", + "language": "python", + "name": "rasterix" + }, + "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.13.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 7ef820d..4c1f641 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,6 +31,7 @@ dependencies = [ "pandas>=2", "numpy>=2", "xarray>=2025", + "xproj", ] dynamic=["version"] [project.optional-dependencies] @@ -42,13 +43,16 @@ rasterize = [ ] exactextract = ["exactextract"] test = [ - "geodatasets", "dask-geopandas", + "geodatasets", + "exactextract", + "ipykernel", + "ipython", + "netCDF4", "odc-geo", "rasterio", "rioxarray", - "exactextract", - "netCDF4" + "ruff", ] [tool.hatch] @@ -82,8 +86,9 @@ run-coverage = "pytest --cov-config=pyproject.toml --cov=pkg --cov-report xml -- run-coverage-html = "pytest --cov-config=pyproject.toml --cov=pkg --cov-report html --cov=src" run-pytest = "run-coverage --no-cov" run-verbose = "run-coverage --verbose" -run-mypy = "mypy src" list-env = "pip list" +install-kernel = "python -m ipykernel install --user --name=rasterix" +lint = "ruff check" [tool.ruff.lint] # E402: module level import not at top of file diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py index 25bfe3c..cf4ff4b 100644 --- a/src/rasterix/raster_index.py +++ b/src/rasterix/raster_index.py @@ -18,6 +18,9 @@ T_Xarray = TypeVar("T_Xarray", "DataArray", "Dataset") +import xproj +import warnings + def assign_index(obj: T_Xarray, *, x_dim: str | None = None, y_dim: str | None = None) -> T_Xarray: if x_dim is None or y_dim is None: @@ -284,7 +287,7 @@ def _filter_dim_indexers(index: WrappedIndex, indexers: Mapping) -> Mapping: return {dim: indexers[dim] for dim in dims if dim in indexers} -class RasterIndex(Index): +class RasterIndex(Index, xproj.ProjIndexMixin): """Xarray index for raster coordinates. RasterIndex is itself a wrapper around one or more Xarray indexes associated @@ -322,6 +325,21 @@ def __init__(self, indexes: Mapping[WrappedIndexCoords, WrappedIndex]): self._wrapped_indexes = dict(indexes) + self._crs = None + + def _proj_get_crs(self): + return self._crs + + def _proj_set_crs(self, spatial_ref, crs): + # note: `spatial_ref` is not used here + print(f"set CRS of index {self!r} to crs={crs}!") + + self._crs = crs + return self + + def _repr_inline_(self, max_width=70): + return f"{type(self).__name__} (crs={self._crs})" + @classmethod def from_transform( cls, affine: Affine, width: int, height: int, x_dim: str = "x", y_dim: str = "y" @@ -385,6 +403,13 @@ def isel(self, indexers: Mapping[Any, int | slice | np.ndarray | Variable]) -> R return None def sel(self, labels: dict[Any, Any], method=None, tolerance=None) -> IndexSelResult: + + if self._crs is not None: + warnings.warn( + f"make sure that indexer labels have CRS {self._crs}!", + UserWarning, + ) + results = [] for coord_names, index in self._wrapped_indexes.items(): diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py index 045ae3a..478c649 100644 --- a/tests/test_raster_index.py +++ b/tests/test_raster_index.py @@ -2,17 +2,34 @@ import rioxarray # noqa import xarray as xr from affine import Affine +import pytest + +from xarray.structure.merge import MergeError from rasterix import RasterIndex, assign_index +def _open_test_raster(): + source = "/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif" + return xr.open_dataarray(source, engine="rasterio") + def test_rectilinear(): - source = "/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif" - da_no_raster_index = xr.open_dataarray(source, engine="rasterio") + da_no_raster_index = _open_test_raster() da_raster_index = assign_index(da_no_raster_index) assert da_raster_index.equals(da_no_raster_index) +def test_different_crs_same_geotransform(): + da = _open_test_raster() + da1 = assign_index(da) + da2 = da1.copy() + da1 = da1.proj.assign_crs(spatial_ref='EPSG:32610') + da2 = da2.proj.assign_crs(spatial_ref='EPSG:32611') + with pytest.raises( + MergeError, match="conflicting values/indexes on objects to be combined for coordinate" + ): + da1 + da2 + # TODO: parameterize over # 1. y points up; # 2. y points down