You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ji...@apache.org on 2026/05/06 06:38:41 UTC
(sedona) branch master updated: [GH-2700] Add 04-flood-snapshot notebook: SAR mask → flood polygon → affected buildings (#2905)
This is an automated email from the ASF dual-hosted git repository.
jiayuasu pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sedona.git
The following commit(s) were added to refs/heads/master by this push:
new edcccc21ff [GH-2700] Add 04-flood-snapshot notebook: SAR mask → flood polygon → affected buildings (#2905)
edcccc21ff is described below
commit edcccc21ff85b2d98d7bc547165f6ed26879d3ac
Author: Jia Yu <ji...@apache.org>
AuthorDate: Tue May 5 23:38:35 2026 -0700
[GH-2700] Add 04-flood-snapshot notebook: SAR mask → flood polygon → affected buildings (#2905)
---
docs/usecases/04-flood-snapshot.ipynb | 368 ++++++++++++++++++++++++++++++++++
1 file changed, 368 insertions(+)
diff --git a/docs/usecases/04-flood-snapshot.ipynb b/docs/usecases/04-flood-snapshot.ipynb
new file mode 100644
index 0000000000..f37abd7a46
--- /dev/null
+++ b/docs/usecases/04-flood-snapshot.ipynb
@@ -0,0 +1,368 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": "<!--\n Licensed to the Apache Software Foundation (ASF) under one\n or more contributor license agreements. See the NOTICE file\n distributed with this work for additional information\n regarding copyright ownership. The ASF licenses this file\n to you under the Apache License, Version 2.0 (the\n \"License\"); you may not use this file except in compliance\n with the License. You may obtain a copy of the License at\n\n http://www.apache.org/licenses/LICENSE-2.0\n\n Unl [...]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1. Connect to Spark through SedonaContext"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sedona.spark import SedonaContext\n",
+ "\n",
+ "config = SedonaContext.builder().master(\"spark://localhost:7077\").getOrCreate()\n",
+ "sedona = SedonaContext.create(config)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2. Synthesize a Sentinel-1-VV-style backscatter raster\n",
+ "\n",
+ "Sentinel-1 SAR over flat water returns very low backscatter (water specularly reflects radar away from the satellite); over land the surface scatters energy back. We simulate that contrast with a single-band `float32` raster: high values everywhere except a circular flood patch in the south-west quadrant which we drop to low values. Written as a tiled GeoTIFF so the Sedona raster reader accepts it."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "\n",
+ "import numpy as np\n",
+ "import rasterio\n",
+ "from rasterio.transform import from_bounds\n",
+ "\n",
+ "WORK = \"/tmp/flood-snapshot\"\n",
+ "os.makedirs(WORK, exist_ok=True)\n",
+ "\n",
+ "AOI = (-91.10, 41.50, -91.00, 41.60)\n",
+ "W = H = 128\n",
+ "transform = from_bounds(*AOI, W, H)\n",
+ "rng = np.random.default_rng(11)\n",
+ "\n",
+ "ys, xs = np.mgrid[0:H, 0:W]\n",
+ "flood_mask = ((xs - 32) ** 2 + (ys - 96) ** 2) < 28**2 # circle in SW corner\n",
+ "backscatter = 0.5 + 0.05 * rng.standard_normal((H, W)) # land\n",
+ "backscatter = np.where(\n",
+ " flood_mask, 0.05 + 0.02 * rng.standard_normal((H, W)), backscatter\n",
+ ")\n",
+ "backscatter = backscatter.clip(0, 1).astype(\"float32\")\n",
+ "\n",
+ "scene_path = f\"{WORK}/sentinel1_vv.tif\"\n",
+ "with rasterio.open(\n",
+ " scene_path,\n",
+ " \"w\",\n",
+ " driver=\"GTiff\",\n",
+ " tiled=True,\n",
+ " blockxsize=128,\n",
+ " blockysize=128,\n",
+ " height=H,\n",
+ " width=W,\n",
+ " count=1,\n",
+ " dtype=\"float32\",\n",
+ " crs=\"EPSG:4326\",\n",
+ " transform=transform,\n",
+ ") as dst:\n",
+ " dst.write(backscatter, 1)\n",
+ " dst.set_band_description(1, \"vv\")\n",
+ "\n",
+ "print(\n",
+ " f\"backscatter min={backscatter.min():.2f} max={backscatter.max():.2f} \"\n",
+ " f\"flood pixels={int(flood_mask.sum())}/{H*W}\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 3. Threshold the backscatter into a water-mask raster\n",
+ "\n",
+ "Single-raster `RS_MapAlgebra` with a comparison expression: every pixel whose VV backscatter is below `0.2` becomes 1 (water), every other pixel becomes 0 (dry). The result is a `uint8` mask raster the same size and CRS as the input."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scene = sedona.read.format(\"raster\").load(scene_path)\n",
+ "scene.createOrReplaceTempView(\"scene\")\n",
+ "\n",
+ "mask = sedona.sql(\"\"\"\n",
+ " SELECT x, y,\n",
+ " RS_MapAlgebra(\n",
+ " rast, 'B',\n",
+ " 'out[0] = rast[0] < 0.2 ? 1 : 0;'\n",
+ " ) AS rast\n",
+ " FROM scene\n",
+ "\"\"\")\n",
+ "mask.cache()\n",
+ "mask.createOrReplaceTempView(\"mask\")\n",
+ "mask.selectExpr(\n",
+ " \"x\",\n",
+ " \"y\",\n",
+ " \"RS_BandPixelType(rast, 1) as dtype\",\n",
+ " \"RS_Width(rast) as w\",\n",
+ " \"RS_Height(rast) as h\",\n",
+ ").show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 4. Vectorize the mask into a single dissolved flood polygon\n",
+ "\n",
+ "`RS_PixelAsPolygons(raster, band)` returns an array of `{geom, val, x, y}` structs \u2014 one per pixel. We `explode` that into rows, keep only the wet pixels (`val = 1`), and dissolve them into a single `MultiPolygon` via `ST_Union_Agg`. A small `ST_Buffer` (0.0005\u00b0, ~50 m) absorbs SAR speckle along the edge of the flood."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": "wet_pixels = sedona.sql(\"\"\"\n SELECT explode(RS_PixelAsPolygons(rast, 1)) AS pixel\n FROM mask\n\"\"\").selectExpr(\"pixel.geom AS geom\", \"pixel.value AS value\").where(\"value = 1\")\nwet_pixels.createOrReplaceTempView(\"wet_pixels\")\n\nflood = sedona.sql(\"\"\"\n SELECT ST_Buffer(ST_Union_Agg(geom), 0.0005) AS geom\n FROM wet_pixels\n\"\"\")\nflood_geom = flood.first()[\"geom\"]\nprint(\n f\"flood polygon: type={flood_geom.geom_type}, \"\n f\"area= [...]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": "## 4. Vectorize the mask into a single dissolved flood polygon\n\n`RS_PixelAsPolygons(raster, band)` returns an array of `{geom, value, x, y}` structs \u2014 one per pixel. We `explode` that into rows, keep only the wet pixels (`value = 1`), and dissolve them into a single `MultiPolygon` via `ST_Union_Agg`. A small `ST_Buffer` (0.0005\u00b0, ~50 m) absorbs SAR speckle along the edge of the flood."
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pyspark.sql import Row\n",
+ "\n",
+ "xmin, ymin, xmax, ymax = AOI\n",
+ "step_x = (xmax - xmin) / 4\n",
+ "step_y = (ymax - ymin) / 4\n",
+ "half = 0.00125\n",
+ "types_grid = [\n",
+ " [\"residential\", \"hospital\", \"residential\", \"school\"],\n",
+ " [\"residential\", \"residential\", \"fire_station\", \"residential\"],\n",
+ " [\"school\", \"residential\", \"residential\", \"hospital\"],\n",
+ " [\"residential\", \"fire_station\", \"residential\", \"residential\"],\n",
+ "]\n",
+ "\n",
+ "building_rows = []\n",
+ "for ix in range(4):\n",
+ " for iy in range(4):\n",
+ " cx = xmin + (ix + 0.5) * step_x\n",
+ " cy = ymin + (iy + 0.5) * step_y\n",
+ " wkt = (\n",
+ " f\"POLYGON(({cx-half} {cy-half}, {cx+half} {cy-half}, \"\n",
+ " f\"{cx+half} {cy+half}, {cx-half} {cy+half}, {cx-half} {cy-half}))\"\n",
+ " )\n",
+ " building_rows.append(Row(bid=f\"B{ix}{iy}\", type=types_grid[iy][ix], wkt=wkt))\n",
+ "\n",
+ "buildings = sedona.createDataFrame(building_rows).selectExpr(\n",
+ " \"bid\", \"type\", \"ST_SetSRID(ST_GeomFromText(wkt), 4326) as geom\"\n",
+ ")\n",
+ "buildings.createOrReplaceTempView(\"buildings\")\n",
+ "buildings.show(16)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 6. Find affected buildings + flag critical infrastructure\n",
+ "\n",
+ "`ST_Intersects(building, flood)` is the canonical \"affected by\" test. We add a `is_critical` flag for any affected building whose type is hospital / school / fire_station \u2014 those are the rows an emergency-management dashboard would surface first."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "affected = sedona.sql(\"\"\"\n",
+ " SELECT b.bid,\n",
+ " b.type,\n",
+ " b.type IN ('hospital', 'school', 'fire_station') AS is_critical,\n",
+ " b.geom\n",
+ " FROM buildings b, flood f\n",
+ " WHERE ST_Intersects(b.geom, f.geom)\n",
+ " ORDER BY is_critical DESC, b.bid\n",
+ "\"\"\")\n",
+ "affected.cache()\n",
+ "affected.createOrReplaceTempView(\"affected\")\n",
+ "affected.select(\"bid\", \"type\", \"is_critical\").show()\n",
+ "\n",
+ "summary = sedona.sql(\"\"\"\n",
+ " SELECT COUNT(*) AS affected_buildings,\n",
+ " SUM(CAST(is_critical AS INT)) AS affected_critical\n",
+ " FROM affected\n",
+ "\"\"\")\n",
+ "summary.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 7. Persist the affected buildings as GeoParquet 1.1\n",
+ "\n",
+ "Auto-populated covering-bbox metadata + projjson CRS (derived from SRID 4326) make the result file ready for downstream tools to filter on bbox without scanning every row. Round-trip read-back confirms the metadata lands."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import shutil\n",
+ "\n",
+ "out_path = f\"{WORK}/affected_buildings.parquet\"\n",
+ "if os.path.exists(out_path):\n",
+ " shutil.rmtree(out_path)\n",
+ "affected.coalesce(1).write.format(\"geoparquet\").save(out_path)\n",
+ "\n",
+ "sedona.read.format(\"geoparquet\").load(out_path).select(\n",
+ " \"bid\", \"type\", \"is_critical\"\n",
+ ").show(truncate=False)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 8. Visualize: backscatter + flood polygon + buildings color-coded\n",
+ "\n",
+ "Single matplotlib axes: VV backscatter as the basemap, the dissolved flood polygon overlaid in blue, building footprints filled by status \u2014 red for affected critical, orange for affected residential, white for unaffected."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.patches as mpatches\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "buildings_pdf = buildings.selectExpr(\"bid\", \"type\", \"ST_AsText(geom) as wkt\").toPandas()\n",
+ "affected_ids = set(r[\"bid\"] for r in affected.select(\"bid\").collect())\n",
+ "critical_affected_ids = set(\n",
+ " r[\"bid\"] for r in affected.where(\"is_critical\").select(\"bid\").collect()\n",
+ ")\n",
+ "\n",
+ "flood_wkt = sedona.sql(\"SELECT ST_AsText(geom) as wkt FROM flood\").first()[\"wkt\"]\n",
+ "\n",
+ "fig, ax = plt.subplots(1, 1, figsize=(7, 7))\n",
+ "extent = (AOI[0], AOI[2], AOI[1], AOI[3])\n",
+ "ax.imshow(backscatter, vmin=0, vmax=1, cmap=\"gray\", extent=extent)\n",
+ "\n",
+ "from shapely import wkt as shapely_wkt\n",
+ "\n",
+ "flood_shapely = shapely_wkt.loads(flood_wkt)\n",
+ "for poly in (\n",
+ " flood_shapely.geoms if hasattr(flood_shapely, \"geoms\") else [flood_shapely]\n",
+ "):\n",
+ " xs_p, ys_p = poly.exterior.xy\n",
+ " ax.fill(\n",
+ " xs_p, ys_p, facecolor=\"#268bd2\", edgecolor=\"#073642\", linewidth=1, alpha=0.5\n",
+ " )\n",
+ "\n",
+ "for _, row in buildings_pdf.iterrows():\n",
+ " bbox = row[\"wkt\"].split(\"((\")[1].split(\"))\")[0]\n",
+ " pts = [tuple(float(x) for x in p.strip().split()) for p in bbox.split(\",\")]\n",
+ " xs_b = [p[0] for p in pts]\n",
+ " ys_b = [p[1] for p in pts]\n",
+ " if row[\"bid\"] in critical_affected_ids:\n",
+ " face = \"#dc322f\" # red\n",
+ " elif row[\"bid\"] in affected_ids:\n",
+ " face = \"#cb4b16\" # orange\n",
+ " else:\n",
+ " face = \"#fdf6e3\" # cream\n",
+ " ax.fill(xs_b, ys_b, facecolor=face, edgecolor=\"black\", linewidth=0.6)\n",
+ " ax.text(\n",
+ " sum(xs_b[:-1]) / 4,\n",
+ " sum(ys_b[:-1]) / 4,\n",
+ " f\"{row['bid']}\\n{row['type'][:4]}\",\n",
+ " ha=\"center\",\n",
+ " va=\"center\",\n",
+ " fontsize=6,\n",
+ " )\n",
+ "\n",
+ "ax.set_title(\n",
+ " f\"Flood snapshot: {len(affected_ids)} buildings affected, \"\n",
+ " f\"{len(critical_affected_ids)} critical\"\n",
+ ")\n",
+ "ax.set_xlabel(\"longitude\")\n",
+ "ax.set_ylabel(\"latitude\")\n",
+ "ax.legend(\n",
+ " handles=[\n",
+ " mpatches.Patch(facecolor=\"#268bd2\", alpha=0.5, label=\"flood extent\"),\n",
+ " mpatches.Patch(facecolor=\"#dc322f\", label=\"affected critical\"),\n",
+ " mpatches.Patch(facecolor=\"#cb4b16\", label=\"affected residential\"),\n",
+ " mpatches.Patch(facecolor=\"#fdf6e3\", edgecolor=\"black\", label=\"unaffected\"),\n",
+ " ],\n",
+ " loc=\"upper right\",\n",
+ " fontsize=8,\n",
+ ")\n",
+ "fig.tight_layout()\n",
+ "fig"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## What just happened?\n",
+ "\n",
+ "An end-to-end disaster-response pipeline in a handful of SQL passes:\n",
+ "\n",
+ "1. Synthesized a SAR-VV-style backscatter raster.\n",
+ "2. Thresholded with single-raster `RS_MapAlgebra` to derive a binary water mask.\n",
+ "3. Vectorized the mask via `RS_PixelAsPolygons` + `explode` + `ST_Union_Agg` and absorbed speckle with `ST_Buffer` \u2014 the canonical raster\u2192vector dissolve in Sedona.\n",
+ "4. Built a Spark DataFrame of building footprints labelled with infrastructure type.\n",
+ "5. `ST_Intersects` spatial join surfaced the buildings inside the flood; a `type IN (...)` flag marked the critical ones.\n",
+ "6. Persisted to GeoParquet 1.1 (auto covering-bbox + projjson CRS), read back, plotted everything on a single matplotlib axes.\n",
+ "\n",
+ "Replace the synthesis cell in section 2 with `sedona.read.format(\"raster\").load(\"<sentinel-1 vv asset URL discovered via sedona.read.format('stac')>\")` and the rest of the pipeline runs unchanged on real Copernicus data."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "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.11"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}