Add lightweight GeoTIFF/COG reader and writer

docs/source/user_guide/multispectral.ipynb

@@ -41,18 +41,7 @@
    },
    "outputs": [],
    "source": [
-    "import datashader as ds\n",
-    "from datashader.colors import Elevation\n",
-    "import datashader.transfer_functions as tf\n",
-    "from datashader.transfer_functions import shade\n",
-    "from datashader.transfer_functions import stack\n",
-    "from datashader.transfer_functions import dynspread\n",
-    "from datashader.transfer_functions import set_background\n",
-    "from datashader.transfer_functions import Images, Image\n",
-    "from datashader.utils import orient_array\n",
-    "import numpy as np\n",
-    "import xarray as xr\n",
-    "import rioxarray"
+    "import datashader as ds\nfrom datashader.colors import Elevation\nimport datashader.transfer_functions as tf\nfrom datashader.transfer_functions import shade\nfrom datashader.transfer_functions import stack\nfrom datashader.transfer_functions import dynspread\nfrom datashader.transfer_functions import set_background\nfrom datashader.transfer_functions import Images, Image\nfrom datashader.utils import orient_array\nimport numpy as np\nimport xarray as xr\nfrom xrspatial.geotiff import read_geotiff"
    ]
   },
   {
@@ -143,23 +132,7 @@
     }
    ],
    "source": [
-    "SCENE_ID = \"LC80030172015001LGN00\"\n",
-    "EXTS = {\n",
-    "    \"blue\": \"B2\",\n",
-    "    \"green\": \"B3\",\n",
-    "    \"red\": \"B4\",\n",
-    "    \"nir\": \"B5\",\n",
-    "}\n",
-    "\n",
-    "cvs = ds.Canvas(plot_width=1024, plot_height=1024)\n",
-    "layers = {}\n",
-    "for name, ext in EXTS.items():\n",
-    "    layer = rioxarray.open_rasterio(f\"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff\").load()[0]\n",
-    "    layer.name = name\n",
-    "    layer = cvs.raster(layer, agg=\"mean\")\n",
-    "    layer.data = orient_array(layer)\n",
-    "    layers[name] = layer\n",
-    "layers"
+    "SCENE_ID = \"LC80030172015001LGN00\"\nEXTS = {\n    \"blue\": \"B2\",\n    \"green\": \"B3\",\n    \"red\": \"B4\",\n    \"nir\": \"B5\",\n}\n\ncvs = ds.Canvas(plot_width=1024, plot_height=1024)\nlayers = {}\nfor name, ext in EXTS.items():\n    layer = read_geotiff(f\"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff\", band=0)\n    layer.name = name\n    layer = cvs.raster(layer, agg=\"mean\")\n    layer.data = orient_array(layer)\n    layers[name] = layer\nlayers"
    ]
   },
   {
@@ -362,7 +335,7 @@
        "}\n",
        "\n",
        ".xr-group-name::before {\n",
-       "  content: \"📁\";\n",
+       "  content: \"\ud83d\udcc1\";\n",
        "  padding-right: 0.3em;\n",
        "}\n",
        "\n",
@@ -425,7 +398,7 @@
        "\n",
        ".xr-section-summary-in + label:before {\n",
        "  display: inline-block;\n",
-       "  content: \"►\";\n",
+       "  content: \"\u25ba\";\n",
        "  font-size: 11px;\n",
        "  width: 15px;\n",
        "  text-align: center;\n",
@@ -436,7 +409,7 @@
        "}\n",
        "\n",
        ".xr-section-summary-in:checked + label:before {\n",
-       "  content: \"▼\";\n",
+       "  content: \"\u25bc\";\n",
        "}\n",
        "\n",
        ".xr-section-summary-in:checked + label > span {\n",

examples/user_guide/25_GLCM_Texture.ipynb

@@ -264,7 +264,7 @@
    "id": "ec79xdunce9",
    "metadata": {},
    "source": [
-    "### Step 1 — Download a Sentinel-2 NIR band\n",
+    "### Step 1 \u2014 Download a Sentinel-2 NIR band\n",
     "\n",
     "We read a 500 x 500 pixel window (5 km x 5 km at 10 m resolution) straight from a\n",
     "Cloud-Optimized GeoTIFF hosted on AWS. The scene is\n",
@@ -282,47 +282,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "import rioxarray\n",
-    "\n",
-    "os.environ['AWS_NO_SIGN_REQUEST'] = 'YES'\n",
-    "os.environ['GDAL_DISABLE_READDIR_ON_OPEN'] = 'EMPTY_DIR'\n",
-    "\n",
-    "COG_URL = (\n",
-    "    'https://sentinel-cogs.s3.us-west-2.amazonaws.com/'\n",
-    "    'sentinel-s2-l2a-cogs/10/S/EG/2023/9/'\n",
-    "    'S2B_10SEG_20230921_0_L2A/B08.tif'\n",
-    ")\n",
-    "\n",
-    "try:\n",
-    "    nir_da = rioxarray.open_rasterio(COG_URL).isel(band=0, y=slice(2100, 2600), x=slice(5300, 5800))\n",
-    "    nir = nir_da.load().values.astype(np.float64)\n",
-    "    print(f'Downloaded NIR band: {nir.shape}, range {nir.min():.0f} to {nir.max():.0f}')\n",
-    "except Exception as exc:\n",
-    "    print(f'Remote read failed ({exc}), using synthetic fallback')\n",
-    "    rng_sat = np.random.default_rng(99)\n",
-    "    nir = np.zeros((500, 500), dtype=np.float64)\n",
-    "    nir[:, 250:] = rng_sat.normal(80, 10, (500, 250)).clip(20, 200)\n",
-    "    nir[:, :250] = rng_sat.normal(1800, 400, (500, 250)).clip(300, 4000)\n",
-    "\n",
-    "satellite = xr.DataArray(nir, dims=['y', 'x'],\n",
-    "                         coords={'y': np.arange(nir.shape[0], dtype=float),\n",
-    "                                 'x': np.arange(nir.shape[1], dtype=float)})\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(7, 7))\n",
-    "satellite.plot.imshow(ax=ax, cmap='gray', vmax=float(np.percentile(nir, 98)),\n",
-    "                      add_colorbar=False)\n",
-    "ax.set_title('Sentinel-2 NIR band')\n",
-    "ax.set_axis_off()\n",
-    "plt.tight_layout()"
+    "import os\nfrom xrspatial.geotiff import read_geotiff\n\n\nCOG_URL = (\n    'https://sentinel-cogs.s3.us-west-2.amazonaws.com/'\n    'sentinel-s2-l2a-cogs/10/S/EG/2023/9/'\n    'S2B_10SEG_20230921_0_L2A/B08.tif'\n)\n\ntry:\n    nir_da = read_geotiff(COG_URL, band=0, window=(2100, 5300, 2600, 5800))\n    nir = nir_da.values.astype(np.float64)\n    print(f'Downloaded NIR band: {nir.shape}, range {nir.min():.0f} to {nir.max():.0f}')\nexcept Exception as exc:\n    print(f'Remote read failed ({exc}), using synthetic fallback')\n    rng_sat = np.random.default_rng(99)\n    nir = np.zeros((500, 500), dtype=np.float64)\n    nir[:, 250:] = rng_sat.normal(80, 10, (500, 250)).clip(20, 200)\n    nir[:, :250] = rng_sat.normal(1800, 400, (500, 250)).clip(300, 4000)\n\nsatellite = xr.DataArray(nir, dims=['y', 'x'],\n                         coords={'y': np.arange(nir.shape[0], dtype=float),\n                                 'x': np.arange(nir.shape[1], dtype=float)})\n\nfig, ax = plt.subplots(figsize=(7, 7))\nsatellite.plot.imshow(ax=ax, cmap='gray', vmax=float(np.percentile(nir, 98)),\n                      add_colorbar=False)\nax.set_title('Sentinel-2 NIR band')\nax.set_axis_off()\nplt.tight_layout()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "joxz7n8olpc",
    "metadata": {},
    "source": [
-    "### Step 2 — Compute GLCM texture features\n",
+    "### Step 2 \u2014 Compute GLCM texture features\n",
     "\n",
     "We pick four metrics that tend to separate water (uniform, high energy, high homogeneity) from land (rough, high contrast):\n",
     "\n",
@@ -485,4 +453,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}

examples/viewshed_gpu.ipynb

@@ -34,7 +34,9 @@
     }
    },
    "outputs": [],
-   "source": "import pandas\nimport matplotlib.pyplot as plt\nimport geopandas as gpd\n\nimport xarray as xr\nimport numpy as np\nimport cupy\nimport rioxarray\n\nimport xrspatial"
+   "source": [
+    "import pandas\nimport matplotlib.pyplot as plt\nimport geopandas as gpd\n\nimport xarray as xr\nimport numpy as np\nimport cupy\nfrom xrspatial.geotiff import read_geotiff\n\nimport xrspatial"
+   ]
   },
   {
    "cell_type": "markdown",
@@ -64,15 +66,7 @@
    },
    "outputs": [],
    "source": [
-    "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n",
-    "\n",
-    "raster = rioxarray.open_rasterio(file_name).sel(band=1).drop_vars('band')\n",
-    "raster.name = 'Colorado Elevation Raster'\n",
-    "\n",
-    "xmin, xmax = raster.x.data.min(), raster.x.data.max()\n",
-    "ymin, ymax = raster.y.data.min(), raster.y.data.max()\n",
-    "\n",
-    "xmin, xmax, ymin, ymax"
+    "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n\nraster = read_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
    ]
   },
   {

examples/xarray-spatial_classification-methods.ipynb

@@ -46,7 +46,9 @@
     }
    },
    "outputs": [],
-   "source": "import xarray as xr\nimport rioxarray\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = rioxarray.open_rasterio(file_name).sel(band=1).drop_vars('band')\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
+   "source": [
+    "import xarray as xr\nfrom xrspatial.geotiff import read_geotiff\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = read_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
+   ]
   },
   {
    "cell_type": "code",

setup.cfg

@@ -23,6 +23,8 @@ install_requires =
     scipy
     xarray
     numpy
+    matplotlib
+    zstandard
 packages = find:
 python_requires = >=3.12
 setup_requires = setuptools_scm

xrspatial/accessor.py

@@ -21,6 +21,33 @@ class XrsSpatialDataArrayAccessor:
     def __init__(self, obj):
         self._obj = obj
 
+    # ---- Plot ----
+
+    def plot(self, **kwargs):
+        """Plot the DataArray, using an embedded TIFF colormap if present.
+
+        For palette/indexed-color GeoTIFFs (read via ``read_geotiff``),
+        the TIFF's color table is applied automatically with correct
+        normalization.  For all other DataArrays, falls through to the
+        standard ``da.plot()``.
+
+        Usage::
+
+            da = read_geotiff('landcover.tif')
+            da.xrs.plot()  # palette colors used automatically
+        """
+        import numpy as np
+        cmap = self._obj.attrs.get('cmap')
+        if cmap is not None and 'cmap' not in kwargs:
+            from matplotlib.colors import BoundaryNorm
+            n_colors = len(cmap.colors)
+            boundaries = np.arange(n_colors + 1) - 0.5
+            norm = BoundaryNorm(boundaries, n_colors)
+            kwargs.setdefault('cmap', cmap)
+            kwargs.setdefault('norm', norm)
+            kwargs.setdefault('add_colorbar', True)
+        return self._obj.plot(**kwargs)
+
     # ---- Surface ----
 
     def slope(self, **kwargs):