Add example notebooks testing performance with large datasets

examples/Grand_Canyon_DEM.ipynb

@@ -0,0 +1,232 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Grand Canyon DEM - Large Dataset Test"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Requires the \"play\" dependency group (`uv sync --group play`)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import geopandas as gpd\n",
+    "import leafmap\n",
+    "import rioxarray\n",
+    "import shapely\n",
+    "from rasterio.warp import transform_bounds\n",
+    "from rio_vrt import build_vrt\n",
+    "\n",
+    "from jupyter_xarray_tiler.titiler import add_data_array"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Download the data!\n",
+    "\n",
+    "⚠️ Manual step!\n",
+    "\n",
+    "Download `DEME_Zone1-Zone15_2021.zip` from <https://www.sciencebase.gov/catalog/item/62ab6770d34e74f0d80eb3af>, and unzip its contents into the `data/` subdirectory.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dem_dir = Path() / \"data\" / \"DEME_Zone1-Zone15_2021\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load all TIFFs\n",
+    "tif_files = sorted(dem_dir.glob(\"*.tif\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "datasets = [rioxarray.open_rasterio(tif) for tif in tif_files]\n",
+    "bboxes = [shapely.box(*ds.rio.bounds()) for ds in datasets]\n",
+    "bboxes_gdf = gpd.GeoDataFrame(geometry=bboxes, crs=datasets[0].rio.crs)\n",
+    "bboxes_gdf.explore()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Use all the TIFs through a VRT w/ Dask\n",
+    "\n",
+    "⚠️ Currently only works at higher zoom levels. When zoomed out, TiTiler 500s with body:\n",
+    "\n",
+    "```json\n",
+    "{\"detail\":\"Maximum array limit 1000000000 reached, trying to put DataArray of (1, 32311, 8408) in memory.\"}\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create VRT (virtual mosaic) and open with lazy loading\n",
+    "vrt_path = dem_dir / \"mosaic.vrt\"\n",
+    "build_vrt(str(vrt_path), [str(f) for f in tif_files])\n",
+    "dem = rioxarray.open_rasterio(vrt_path, chunks=\"auto\")\n",
+    "dem"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Calculate map center and zoom from data extent\n",
+    "bounds = dem.rio.bounds()\n",
+    "lon_min, lat_min, lon_max, lat_max = transform_bounds(dem.rio.crs, \"EPSG:4326\", *bounds)\n",
+    "\n",
+    "center = [(lat_min + lat_max) / 2, (lon_min + lon_max) / 2]\n",
+    "max_diff = max(lat_max - lat_min, lon_max - lon_min)\n",
+    "zoom = math.floor(math.log2(360 / max_diff)) - 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Add to TiTiler\n",
+    "url = await add_data_array(dem, rescale=(0, 1000), colormap_name=\"terrain\")\n",
+    "\n",
+    "# Create map centered on data\n",
+    "m = leafmap.Map(center=center, zoom=zoom)\n",
+    "m"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Add to map\n",
+    "m.add_tile_layer(url=url, name=\"DEM\", attribution=\"\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Use all the TIFs through `open_mfdataset` w/ Dask"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# TODO"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Single TIF test  (without Dask)\n",
+    "\n",
+    "Loading a single TIF works correctly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load the first tif into memory without dask\n",
+    "single_tif = rioxarray.open_rasterio(tif_files[0])\n",
+    "\n",
+    "# Calculate map bounds\n",
+    "bounds = single_tif.rio.bounds()\n",
+    "lon_min, lat_min, lon_max, lat_max = transform_bounds(\n",
+    "    single_tif.rio.crs, \"EPSG:4326\", *bounds\n",
+    ")\n",
+    "center = [(lat_min + lat_max) / 2, (lon_min + lon_max) / 2]\n",
+    "zoom = math.floor(math.log2(360 / max(lat_max - lat_min, lon_max - lon_min))) - 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# add to TiTiler\n",
+    "url = await add_data_array(single_tif, rescale=(0, 1000), colormap_name=\"terrain\")\n",
+    "\n",
+    "# Display map\n",
+    "m = leafmap.Map(center=center, zoom=zoom)\n",
+    "m"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Add layer to map\n",
+    "m.add_tile_layer(url=url, name=\"Single TIF\", attribution=\"\")"
+   ]
+  }
+ ],
+ "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.14.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}