Rename GeoTIFF API to xarray conventions (#1047)

README.md

@@ -139,24 +139,30 @@ Native GeoTIFF and Cloud Optimized GeoTIFF reader/writer. No GDAL required.
 
 | Name | Description | NumPy | Dask | CuPy GPU | Dask+CuPy GPU | Cloud |
 |:-----|:------------|:-----:|:----:|:--------:|:-------------:|:-----:|
-| [read_geotiff](xrspatial/geotiff/__init__.py) | Read GeoTIFF / COG / VRT | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
-| [write_geotiff](xrspatial/geotiff/__init__.py) | Write DataArray as GeoTIFF / COG | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
+| [open_geotiff](xrspatial/geotiff/__init__.py) | Read GeoTIFF / COG / VRT | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
+| [to_geotiff](xrspatial/geotiff/__init__.py) | Write DataArray as GeoTIFF / COG | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
 | [write_vrt](xrspatial/geotiff/__init__.py) | Generate VRT mosaic from GeoTIFFs | ✅️ | | | | |
 
-`read_geotiff` and `write_geotiff` auto-dispatch to the correct backend:
+`open_geotiff` and `to_geotiff` auto-dispatch to the correct backend:
 
 ```python
-read_geotiff('dem.tif')                              # NumPy
-read_geotiff('dem.tif', chunks=512)                  # Dask
-read_geotiff('dem.tif', gpu=True)                    # CuPy (nvCOMP + GDS)
-read_geotiff('dem.tif', gpu=True, chunks=512)        # Dask + CuPy
-read_geotiff('https://example.com/cog.tif')          # HTTP COG
-read_geotiff('s3://bucket/dem.tif')                  # Cloud (S3/GCS/Azure)
-read_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected)
-
-write_geotiff(cupy_array, 'out.tif')                 # auto-detects GPU
-write_geotiff(data, 'out.tif', gpu=True)             # force GPU compress
+from xrspatial.geotiff import open_geotiff, to_geotiff
+
+open_geotiff('dem.tif')                              # NumPy
+open_geotiff('dem.tif', chunks=512)                  # Dask
+open_geotiff('dem.tif', gpu=True)                    # CuPy (nvCOMP + GDS)
+open_geotiff('dem.tif', gpu=True, chunks=512)        # Dask + CuPy
+open_geotiff('https://example.com/cog.tif')          # HTTP COG
+open_geotiff('s3://bucket/dem.tif')                  # Cloud (S3/GCS/Azure)
+open_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected)
+
+to_geotiff(cupy_array, 'out.tif')                    # auto-detects GPU
+to_geotiff(data, 'out.tif', gpu=True)                # force GPU compress
 write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
+
+# Accessor methods
+da.xrs.to_geotiff('out.tif', compression='lzw')     # write from DataArray
+ds.xrs.open_geotiff('large_dem.tif')                 # read windowed to Dataset extent
 ```
 
 **Compression codecs:** Deflate, LZW (Numba JIT), ZSTD, PackBits, JPEG (Pillow), JPEG 2000 (glymur), uncompressed
@@ -461,10 +467,10 @@ Importing `xrspatial` registers an `.xrs` accessor on DataArrays and Datasets, g
 
 ```python
 import xrspatial
-from xrspatial.geotiff import read_geotiff
+from xrspatial.geotiff import open_geotiff
 
 # Read a GeoTIFF (no GDAL required)
-elevation = read_geotiff('dem.tif')
+elevation = open_geotiff('dem.tif')
 
 # Surface analysis — call operations directly on the DataArray
 slope = elevation.xrs.slope()

docs/source/user_guide/multispectral.ipynb

@@ -41,7 +41,7 @@
    },
    "outputs": [],
    "source": [
-    "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"
+    "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 open_geotiff"
    ]
   },
   {
@@ -132,7 +132,7 @@
     }
    ],
    "source": [
-    "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"
+    "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 = open_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"
    ]
   },
   {

examples/user_guide/25_GLCM_Texture.ipynb

@@ -282,7 +282,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "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()"
+    "import os\nfrom xrspatial.geotiff import open_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 = open_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()"
    ]
   },
   {

examples/user_guide/35_JPEG2000_Compression.ipynb

@@ -9,7 +9,7 @@
   {
    "cell_type": "code",
    "id": "kamu534xsm",
-   "source": "import numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport tempfile\nimport os\n\nfrom xrspatial.geotiff import read_geotiff, write_geotiff",
+   "source": "import numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport tempfile\nimport os\n\nfrom xrspatial.geotiff import open_geotiff, to_geotiff",
    "metadata": {},
    "execution_count": null,
    "outputs": []
@@ -31,27 +31,27 @@
   {
    "cell_type": "markdown",
    "id": "8tsuyr3jbay",
-   "source": "## Write with JPEG 2000 (lossless)\n\nPass `compression='jpeg2000'` to `write_geotiff`. The default is lossless encoding.",
+   "source": "## Write with JPEG 2000 (lossless)\n\nPass `compression='jpeg2000'` to `to_geotiff`. The default is lossless encoding.",
    "metadata": {}
   },
   {
    "cell_type": "code",
    "id": "ystjp6v30d",
-   "source": "tmpdir = tempfile.mkdtemp(prefix='j2k_demo_')\n\n# Write with JPEG 2000 compression\nj2k_path = os.path.join(tmpdir, 'elevation_j2k.tif')\nwrite_geotiff(da, j2k_path, compression='jpeg2000')\n\n# Compare file sizes with deflate\ndeflate_path = os.path.join(tmpdir, 'elevation_deflate.tif')\nwrite_geotiff(da, deflate_path, compression='deflate')\n\nnone_path = os.path.join(tmpdir, 'elevation_none.tif')\nwrite_geotiff(da, none_path, compression='none')\n\nj2k_size = os.path.getsize(j2k_path)\ndeflate_size = os.path.getsize(deflate_path)\nnone_size = os.path.getsize(none_path)\n\nprint(f\"Uncompressed:  {none_size:>8,} bytes\")\nprint(f\"Deflate:       {deflate_size:>8,} bytes  ({deflate_size/none_size:.1%} of original)\")\nprint(f\"JPEG 2000:     {j2k_size:>8,} bytes  ({j2k_size/none_size:.1%} of original)\")",
+   "source": "tmpdir = tempfile.mkdtemp(prefix='j2k_demo_')\n\n# Write with JPEG 2000 compression\nj2k_path = os.path.join(tmpdir, 'elevation_j2k.tif')\nto_geotiff(da, j2k_path, compression='jpeg2000')\n\n# Compare file sizes with deflate\ndeflate_path = os.path.join(tmpdir, 'elevation_deflate.tif')\nto_geotiff(da, deflate_path, compression='deflate')\n\nnone_path = os.path.join(tmpdir, 'elevation_none.tif')\nto_geotiff(da, none_path, compression='none')\n\nj2k_size = os.path.getsize(j2k_path)\ndeflate_size = os.path.getsize(deflate_path)\nnone_size = os.path.getsize(none_path)\n\nprint(f\"Uncompressed:  {none_size:>8,} bytes\")\nprint(f\"Deflate:       {deflate_size:>8,} bytes  ({deflate_size/none_size:.1%} of original)\")\nprint(f\"JPEG 2000:     {j2k_size:>8,} bytes  ({j2k_size/none_size:.1%} of original)\")",
    "metadata": {},
    "execution_count": null,
    "outputs": []
   },
   {
    "cell_type": "markdown",
    "id": "89y9zun97nb",
-   "source": "## Read it back and verify lossless roundtrip\n\n`read_geotiff` auto-detects the compression from the TIFF header. No special arguments needed.",
+   "source": "## Read it back and verify lossless roundtrip\n\n`open_geotiff` auto-detects the compression from the TIFF header. No special arguments needed.",
    "metadata": {}
   },
   {
    "cell_type": "code",
    "id": "8vf9ljxkx03",
-   "source": "# Read back and check lossless roundtrip\nda_read = read_geotiff(j2k_path)\n\nprint(f\"Shape: {da_read.shape}\")\nprint(f\"Dtype: {da_read.dtype}\")\nprint(f\"CRS:   {da_read.attrs.get('crs')}\")\nprint(f\"Exact match: {np.array_equal(da_read.values, terrain)}\")\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 5))\nda.plot(ax=axes[0], cmap='terrain')\naxes[0].set_title('Original')\nda_read.plot(ax=axes[1], cmap='terrain')\naxes[1].set_title('After JPEG 2000 roundtrip')\nplt.tight_layout()\nplt.show()",
+   "source": "# Read back and check lossless roundtrip\nda_read = open_geotiff(j2k_path)\n\nprint(f\"Shape: {da_read.shape}\")\nprint(f\"Dtype: {da_read.dtype}\")\nprint(f\"CRS:   {da_read.attrs.get('crs')}\")\nprint(f\"Exact match: {np.array_equal(da_read.values, terrain)}\")\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 5))\nda.plot(ax=axes[0], cmap='terrain')\naxes[0].set_title('Original')\nda_read.plot(ax=axes[1], cmap='terrain')\naxes[1].set_title('After JPEG 2000 roundtrip')\nplt.tight_layout()\nplt.show()",
    "metadata": {},
    "execution_count": null,
    "outputs": []
@@ -65,15 +65,15 @@
   {
    "cell_type": "code",
    "id": "mgv9xhsrcen",
-   "source": "# Create a 3-band uint8 image\nrgb = np.zeros((128, 128, 3), dtype=np.uint8)\nrgb[:, :, 0] = np.linspace(0, 255, 128).astype(np.uint8)[None, :]  # red gradient\nrgb[:, :, 1] = np.linspace(0, 255, 128).astype(np.uint8)[:, None]  # green gradient\nrgb[:, :, 2] = 128  # constant blue\n\nda_rgb = xr.DataArray(\n    rgb, dims=['y', 'x', 'band'],\n    coords={'y': np.arange(128), 'x': np.arange(128), 'band': [0, 1, 2]},\n)\n\nrgb_path = os.path.join(tmpdir, 'rgb_j2k.tif')\nwrite_geotiff(da_rgb, rgb_path, compression='jpeg2000')\n\nda_rgb_read = read_geotiff(rgb_path)\nprint(f\"RGB shape: {da_rgb_read.shape}, dtype: {da_rgb_read.dtype}\")\nprint(f\"Exact match: {np.array_equal(da_rgb_read.values, rgb)}\")\n\nfig, axes = plt.subplots(1, 2, figsize=(10, 4))\naxes[0].imshow(rgb)\naxes[0].set_title('Original RGB')\naxes[1].imshow(da_rgb_read.values)\naxes[1].set_title('After J2K roundtrip')\nplt.tight_layout()\nplt.show()",
+   "source": "# Create a 3-band uint8 image\nrgb = np.zeros((128, 128, 3), dtype=np.uint8)\nrgb[:, :, 0] = np.linspace(0, 255, 128).astype(np.uint8)[None, :]  # red gradient\nrgb[:, :, 1] = np.linspace(0, 255, 128).astype(np.uint8)[:, None]  # green gradient\nrgb[:, :, 2] = 128  # constant blue\n\nda_rgb = xr.DataArray(\n    rgb, dims=['y', 'x', 'band'],\n    coords={'y': np.arange(128), 'x': np.arange(128), 'band': [0, 1, 2]},\n)\n\nrgb_path = os.path.join(tmpdir, 'rgb_j2k.tif')\nto_geotiff(da_rgb, rgb_path, compression='jpeg2000')\n\nda_rgb_read = open_geotiff(rgb_path)\nprint(f\"RGB shape: {da_rgb_read.shape}, dtype: {da_rgb_read.dtype}\")\nprint(f\"Exact match: {np.array_equal(da_rgb_read.values, rgb)}\")\n\nfig, axes = plt.subplots(1, 2, figsize=(10, 4))\naxes[0].imshow(rgb)\naxes[0].set_title('Original RGB')\naxes[1].imshow(da_rgb_read.values)\naxes[1].set_title('After J2K roundtrip')\nplt.tight_layout()\nplt.show()",
    "metadata": {},
    "execution_count": null,
    "outputs": []
   },
   {
    "cell_type": "markdown",
    "id": "zzga5hc3a99",
-   "source": "## GPU acceleration\n\nOn systems with nvJPEG2000 installed (CUDA toolkit, RAPIDS environments), pass `gpu=True` to use GPU-accelerated J2K encode/decode. The API is the same -- it falls back to CPU automatically if the library isn't found.\n\n```python\n# GPU write (nvJPEG2000 if available, else CPU fallback)\nwrite_geotiff(cupy_data, \"output.tif\", compression=\"jpeg2000\", gpu=True)\n\n# GPU read (nvJPEG2000 decode if available)\nda = read_geotiff(\"satellite.tif\", gpu=True)\n```",
+   "source": "## GPU acceleration\n\nOn systems with nvJPEG2000 installed (CUDA toolkit, RAPIDS environments), pass `gpu=True` to use GPU-accelerated J2K encode/decode. The API is the same -- it falls back to CPU automatically if the library isn't found.\n\n```python\n# GPU write (nvJPEG2000 if available, else CPU fallback)\nto_geotiff(cupy_data, \"output.tif\", compression=\"jpeg2000\", gpu=True)\n\n# GPU read (nvJPEG2000 decode if available)\nda = open_geotiff(\"satellite.tif\", gpu=True)\n```",
    "metadata": {}
   },
   {

examples/viewshed_gpu.ipynb

@@ -35,7 +35,7 @@
    },
    "outputs": [],
    "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"
+    "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 open_geotiff\n\nimport xrspatial"
    ]
   },
   {
@@ -66,7 +66,7 @@
    },
    "outputs": [],
    "source": [
-    "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"
+    "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n\nraster = open_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

@@ -47,7 +47,7 @@
    },
    "outputs": [],
    "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"
+    "import xarray as xr\nfrom xrspatial.geotiff import open_geotiff\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = open_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"
    ]
   },
   {