|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "id": "0", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "# Coordinate Reference System (CRS) Management\n", |
| 9 | + "\n", |
| 10 | + "Rasterix uses [XProj](https://xproj.readthedocs.io/en/latest/) to store and manage the CRS of a dataset" |
| 11 | + ] |
| 12 | + }, |
| 13 | + { |
| 14 | + "cell_type": "code", |
| 15 | + "execution_count": null, |
| 16 | + "id": "1", |
| 17 | + "metadata": {}, |
| 18 | + "outputs": [], |
| 19 | + "source": [ |
| 20 | + "# GDAL settings for remote data\n", |
| 21 | + "import os\n", |
| 22 | + "\n", |
| 23 | + "import xarray as xr\n", |
| 24 | + "import xproj\n", |
| 25 | + "\n", |
| 26 | + "import rasterix\n", |
| 27 | + "\n", |
| 28 | + "os.environ[\"VSICURL_PC_URL_SIGNING\"] = \"YES\"\n", |
| 29 | + "os.environ[\"GDAL_DISABLE_READDIR_ON_OPEN\"] = \"EMPTY_DIR\"\n", |
| 30 | + "\n", |
| 31 | + "xr.set_options(display_expand_indexes=True);" |
| 32 | + ] |
| 33 | + }, |
| 34 | + { |
| 35 | + "cell_type": "code", |
| 36 | + "execution_count": null, |
| 37 | + "id": "2", |
| 38 | + "metadata": {}, |
| 39 | + "outputs": [], |
| 40 | + "source": [ |
| 41 | + "url = \"https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif\"\n", |
| 42 | + "da = xr.open_dataarray(url, engine=\"rasterio\")\n", |
| 43 | + "\n", |
| 44 | + "da = rasterix.assign_index(da).set_xindex(\"spatial_ref\", xproj.CRSIndex)\n", |
| 45 | + "\n", |
| 46 | + "da" |
| 47 | + ] |
| 48 | + }, |
| 49 | + { |
| 50 | + "cell_type": "code", |
| 51 | + "execution_count": null, |
| 52 | + "id": "3", |
| 53 | + "metadata": {}, |
| 54 | + "outputs": [], |
| 55 | + "source": [ |
| 56 | + "# Make sure we can save this back to a geotiff\n", |
| 57 | + "# Hmmm, this warning isn't clear to me\n", |
| 58 | + "da.rio.to_raster(\"/tmp/S_20240101_concentration_v3.0.tif\")" |
| 59 | + ] |
| 60 | + }, |
| 61 | + { |
| 62 | + "cell_type": "code", |
| 63 | + "execution_count": null, |
| 64 | + "id": "4", |
| 65 | + "metadata": {}, |
| 66 | + "outputs": [], |
| 67 | + "source": [ |
| 68 | + "!wget -nc https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif\n", |
| 69 | + "!ls -l S_20240101_concentration_v3.0.tif" |
| 70 | + ] |
| 71 | + }, |
| 72 | + { |
| 73 | + "cell_type": "code", |
| 74 | + "execution_count": null, |
| 75 | + "id": "5", |
| 76 | + "metadata": {}, |
| 77 | + "outputs": [], |
| 78 | + "source": [ |
| 79 | + "# Are the tiffs the same?\n", |
| 80 | + "!ls -l /tmp/S_20240101_concentration_v3.0.tif" |
| 81 | + ] |
| 82 | + }, |
| 83 | + { |
| 84 | + "cell_type": "markdown", |
| 85 | + "id": "6", |
| 86 | + "metadata": {}, |
| 87 | + "source": [ |
| 88 | + "## Two rasters with same GeoTransform but different CRS\n", |
| 89 | + "\n", |
| 90 | + "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", |
| 91 | + "\n", |
| 92 | + "We want Xarray to raise an error in this case before attempting any computations." |
| 93 | + ] |
| 94 | + }, |
| 95 | + { |
| 96 | + "cell_type": "code", |
| 97 | + "execution_count": null, |
| 98 | + "id": "7", |
| 99 | + "metadata": {}, |
| 100 | + "outputs": [], |
| 101 | + "source": [ |
| 102 | + "def load_raster(href, overview_level=3, masked=True):\n", |
| 103 | + " da = xr.open_dataarray(\n", |
| 104 | + " href, engine=\"rasterio\", masked=masked, open_kwargs=dict(overview_level=overview_level), chunks=\"auto\"\n", |
| 105 | + " ).squeeze()\n", |
| 106 | + "\n", |
| 107 | + " return rasterix.assign_index(da).set_xindex(\"spatial_ref\", xproj.CRSIndex)\n", |
| 108 | + "\n", |
| 109 | + "\n", |
| 110 | + "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", |
| 111 | + "da1 = load_raster(url)\n", |
| 112 | + "\n", |
| 113 | + "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", |
| 114 | + "da2 = load_raster(url)\n", |
| 115 | + "\n", |
| 116 | + "da1 + da2" |
| 117 | + ] |
| 118 | + }, |
| 119 | + { |
| 120 | + "cell_type": "code", |
| 121 | + "execution_count": null, |
| 122 | + "id": "8", |
| 123 | + "metadata": {}, |
| 124 | + "outputs": [], |
| 125 | + "source": [] |
| 126 | + } |
| 127 | + ], |
| 128 | + "metadata": { |
| 129 | + "language_info": { |
| 130 | + "codemirror_mode": { |
| 131 | + "name": "ipython", |
| 132 | + "version": 3 |
| 133 | + }, |
| 134 | + "file_extension": ".py", |
| 135 | + "mimetype": "text/x-python", |
| 136 | + "name": "python", |
| 137 | + "nbconvert_exporter": "python", |
| 138 | + "pygments_lexer": "ipython3" |
| 139 | + } |
| 140 | + }, |
| 141 | + "nbformat": 4, |
| 142 | + "nbformat_minor": 5 |
| 143 | +} |
0 commit comments