{ "cells": [ { "cell_type": "markdown", "id": "529a0622", "metadata": {}, "source": [ "# GeoJSON Tiler for Risk Analyzer Workflows\n", "\n", "Cut a large GeoJSON polygon (an AOI) into smaller, **grid-aligned tiles** that can be\n", "submitted to the Risk Analyzer one at a time, then reassemble the per-tile results\n", "into a single map.\n", "\n", "**Why aligned tiles matter:** the grid is snapped to a fixed origin in a projected\n", "CRS, so tile edges land on exact multiples of the tile size. That means:\n", "\n", "1. Tiles from different runs (or different counties) line up with each other.\n", "2. The union of the clipped tiles exactly reproduces the original polygon — no gaps, no overlaps.\n", "3. Per-tile analysis outputs can be mosaicked back together cleanly.\n", "\n", "**Workflow**\n", "1. Configure inputs (Section 1)\n", "2. Load & inspect the AOI (Section 2)\n", "3. Build the aligned grid and clip it to the AOI (Section 3)\n", "4. Validate coverage (Section 4)\n", "5. Export one GeoJSON per tile for the Risk Analyzer (Section 5)\n", "6. Preview tiles on an interactive map (Section 6)\n", "7. After analysis: reassemble per-tile results onto one map (Section 7)\n" ] }, { "cell_type": "markdown", "id": "1886bb6a", "metadata": {}, "source": [ "## 1. Configuration\n", "\n", "Edit these values, then run everything below. Rules of thumb for `TILE_SIZE_KM`:\n", "- If the Risk Analyzer fails on a whole county (~50–70 km across), try 20–25 km tiles.\n", "- If tiles still time out on NASA data download, halve the tile size — because the grid\n", " origin is fixed, smaller tiles will nest inside the larger ones.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a3e37432", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "# --- Inputs -----------------------------------------------------------------\n", "INPUT_GEOJSON = \"Kerr_County.geojson\" # path to the AOI to cut\n", "AOI_NAME = \"kerr_county\" # short slug used in output filenames\n", "\n", "# --- Tiling -----------------------------------------------------------------\n", "TILE_SIZE_KM = 20 # tile edge length, in kilometers\n", "GRID_CRS = \"EPSG:5070\" # equal-area CONUS Albers; good default for Texas.\n", " # Any projected CRS in meters works (e.g. a UTM zone).\n", "GRID_ORIGIN = (0.0, 0.0) # fixed grid origin in GRID_CRS coordinates.\n", " # Leave at (0, 0) so every run/AOI shares one global grid.\n", "\n", "MIN_AREA_FRAC = 0.0 # drop clipped tiles smaller than this fraction of a full\n", " # tile (0.0 keeps every sliver -> exact coverage; use e.g.\n", " # 0.01 only if the Risk Analyzer rejects tiny polygons)\n", "\n", "# --- Outputs ----------------------------------------------------------------\n", "OUT_DIR = Path(f\"tiles_{AOI_NAME}_{TILE_SIZE_KM}km\") # per-tile geojsons\n", "RESULTS_DIR = Path(f\"results_{AOI_NAME}_{TILE_SIZE_KM}km\") # drop Risk Analyzer outputs here\n", "OUT_DIR.mkdir(exist_ok=True)\n", "RESULTS_DIR.mkdir(exist_ok=True)\n", "\n", "print(f\"Tiles will be written to: {OUT_DIR}/\")\n", "print(f\"Put Risk Analyzer results (geojson) in: {RESULTS_DIR}/\")" ] }, { "cell_type": "markdown", "id": "d13e3061", "metadata": {}, "source": [ "## 2. Load and inspect the AOI\n", "\n", "The AOI is loaded, assigned WGS84 if the file has no CRS (plain GeoJSON is WGS84 by\n", "spec), dissolved to a single geometry, and reprojected to the grid CRS for tiling.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cd92f804", "metadata": {}, "outputs": [], "source": [ "import geopandas as gpd\n", "import pandas as pd\n", "from shapely.geometry import box\n", "from shapely.ops import unary_union\n", "import math, json\n", "\n", "aoi = gpd.read_file(INPUT_GEOJSON)\n", "if aoi.crs is None:\n", " aoi = aoi.set_crs(\"EPSG:4326\") # GeoJSON spec default\n", "\n", "print(f\"Features: {len(aoi)} CRS: {aoi.crs}\")\n", "display(aoi.drop(columns=\"geometry\").head())\n", "\n", "# Dissolve to a single AOI geometry and project to the grid CRS\n", "aoi_geom_wgs = unary_union(aoi.geometry.values)\n", "aoi_proj = gpd.GeoSeries([aoi_geom_wgs], crs=\"EPSG:4326\").to_crs(GRID_CRS)\n", "aoi_geom = aoi_proj.iloc[0].buffer(0) # buffer(0) fixes minor invalidities\n", "\n", "area_km2 = aoi_geom.area / 1e6\n", "minx, miny, maxx, maxy = aoi_geom.bounds\n", "print(f\"AOI area: {area_km2:,.0f} km2 ({area_km2/2.58999:,.0f} sq mi)\")\n", "print(f\"AOI extent: {(maxx-minx)/1000:,.1f} km x {(maxy-miny)/1000:,.1f} km\")\n", "est = area_km2 / TILE_SIZE_KM**2\n", "print(f\"Expect roughly {est:,.0f}+ tiles at {TILE_SIZE_KM} km\")" ] }, { "cell_type": "markdown", "id": "6b324c53", "metadata": {}, "source": [ "## 3. Build the aligned grid and clip to the AOI\n", "\n", "Grid cells are placed at **exact multiples of the tile size from the fixed origin**\n", "(`floor(min / size) * size`), not at the AOI's bounding-box corner. This is the\n", "alignment guarantee: any AOI tiled with the same `GRID_CRS`, `GRID_ORIGIN`, and\n", "`TILE_SIZE_KM` falls on the same global grid. Tile IDs encode the grid column/row\n", "(`cXXXX_rYYYY`), so the same ID always refers to the same patch of ground.\n", "\n", "Each grid cell is then **clipped to the AOI boundary**, so edge tiles follow the\n", "county line and the tiles collectively reproduce the original polygon exactly.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c0a4f3f1", "metadata": {}, "outputs": [], "source": [ "tile_m = TILE_SIZE_KM * 1000.0\n", "ox, oy = GRID_ORIGIN\n", "\n", "# Snap the AOI bounds outward to the global grid\n", "col0 = math.floor((minx - ox) / tile_m)\n", "col1 = math.ceil ((maxx - ox) / tile_m)\n", "row0 = math.floor((miny - oy) / tile_m)\n", "row1 = math.ceil ((maxy - oy) / tile_m)\n", "\n", "records = []\n", "for col in range(col0, col1):\n", " for row in range(row0, row1):\n", " cell = box(ox + col*tile_m, oy + row*tile_m,\n", " ox + (col+1)*tile_m, oy + (row+1)*tile_m)\n", " if not cell.intersects(aoi_geom):\n", " continue\n", " clipped = cell.intersection(aoi_geom)\n", " if clipped.is_empty or clipped.area == 0:\n", " continue\n", " records.append({\n", " \"tile_id\": f\"{AOI_NAME}_c{col:+05d}_r{row:+05d}\",\n", " \"grid_col\": col,\n", " \"grid_row\": row,\n", " \"area_km2\": clipped.area / 1e6,\n", " \"frac_full\": clipped.area / cell.area,\n", " \"geometry\": clipped,\n", " })\n", "\n", "tiles = gpd.GeoDataFrame(records, crs=GRID_CRS)\n", "\n", "# Optionally drop slivers (keep at 0.0 for exact coverage)\n", "if MIN_AREA_FRAC > 0:\n", " before = len(tiles)\n", " tiles = tiles[tiles[\"frac_full\"] >= MIN_AREA_FRAC].copy()\n", " print(f\"Dropped {before - len(tiles)} sliver tiles (< {MIN_AREA_FRAC:.0%} of a full tile)\")\n", "\n", "tiles = tiles.sort_values([\"grid_row\", \"grid_col\"]).reset_index(drop=True)\n", "print(f\"{len(tiles)} tiles\")\n", "display(tiles.drop(columns=\"geometry\").head(10))" ] }, { "cell_type": "markdown", "id": "d801ab14", "metadata": {}, "source": [ "## 4. Validate coverage\n", "\n", "Confirms the tiles reassemble to the original AOI: total tile area should match the\n", "AOI area, and the symmetric difference between the tile union and the AOI should be\n", "~zero. If `MIN_AREA_FRAC > 0`, the dropped-sliver area shows up here — that's\n", "expected, and tells you exactly how much ground you're excluding.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9bcedbde", "metadata": {}, "outputs": [], "source": [ "tile_union = unary_union(tiles.geometry.values)\n", "sym_diff_km2 = tile_union.symmetric_difference(aoi_geom).area / 1e6\n", "\n", "print(f\"AOI area: {aoi_geom.area/1e6:,.3f} km2\")\n", "print(f\"Sum of tiles: {tiles['area_km2'].sum():,.3f} km2\")\n", "print(f\"Union vs AOI symmetric difference: {sym_diff_km2:.6f} km2\")\n", "\n", "if sym_diff_km2 < 1e-3:\n", " print(\"PASS - tiles exactly reproduce the AOI (no gaps, no overlaps)\")\n", "elif MIN_AREA_FRAC > 0:\n", " print(f\"NOTE - {sym_diff_km2:.3f} km2 excluded by the sliver filter\")\n", "else:\n", " print(\"WARN - unexpected mismatch; inspect the geometry\")" ] }, { "cell_type": "markdown", "id": "7a5c2d57", "metadata": {}, "source": [ "## 5. Export one GeoJSON per tile\n", "\n", "Each tile is written as its own WGS84 FeatureCollection, ready to load into the\n", "Risk Analyzer. A `tile_index.geojson` (all tiles in one file, with IDs) and a\n", "`tile_index.csv` are also written so you can track which tiles have been analyzed.\n", "\n", "**Suggested convention:** when the Risk Analyzer produces a result for\n", "`kerr_county_c-0012_r+0045.geojson`, save the output into the results folder with\n", "the same `tile_id` in the filename — Section 7 uses that to stitch results together.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c286f9ee", "metadata": {}, "outputs": [], "source": [ "tiles_wgs = tiles.to_crs(\"EPSG:4326\")\n", "\n", "for _, t in tiles_wgs.iterrows():\n", " single = gpd.GeoDataFrame([t], crs=\"EPSG:4326\")\n", " single.to_file(OUT_DIR / f\"{t.tile_id}.geojson\", driver=\"GeoJSON\")\n", "\n", "tiles_wgs.to_file(OUT_DIR / \"tile_index.geojson\", driver=\"GeoJSON\")\n", "tiles_wgs.drop(columns=\"geometry\").to_csv(OUT_DIR / \"tile_index.csv\", index=False)\n", "\n", "print(f\"Wrote {len(tiles_wgs)} tile files + tile_index.geojson + tile_index.csv to {OUT_DIR}/\")\n", "sorted(p.name for p in OUT_DIR.glob(\"*.geojson\"))[:8]" ] }, { "cell_type": "markdown", "id": "8c680b4b", "metadata": {}, "source": [ "## 6. Preview the tiles on an interactive map\n", "\n", "Hover a tile for its ID and area. Use this to sanity-check the cut before running\n", "analyses, and to decide whether the tile size is right.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eda0186e", "metadata": {}, "outputs": [], "source": [ "import folium\n", "\n", "center = [tiles_wgs.geometry.union_all().centroid.y,\n", " tiles_wgs.geometry.union_all().centroid.x]\n", "m = folium.Map(location=center, zoom_start=9, tiles=\"cartodbpositron\")\n", "\n", "folium.GeoJson(\n", " aoi.to_crs(\"EPSG:4326\"),\n", " name=\"Original AOI\",\n", " style_function=lambda f: {\"color\": \"#333\", \"weight\": 3, \"fill\": False, \"dashArray\": \"6\"},\n", ").add_to(m)\n", "\n", "folium.GeoJson(\n", " tiles_wgs,\n", " name=\"Tiles\",\n", " style_function=lambda f: {\"color\": \"#1f6feb\", \"weight\": 1.5,\n", " \"fillColor\": \"#1f6feb\", \"fillOpacity\": 0.12},\n", " highlight_function=lambda f: {\"weight\": 3, \"fillOpacity\": 0.35},\n", " tooltip=folium.GeoJsonTooltip(fields=[\"tile_id\", \"area_km2\"],\n", " aliases=[\"Tile\", \"Area (km2)\"]),\n", ").add_to(m)\n", "\n", "folium.LayerControl().add_to(m)\n", "m" ] }, { "cell_type": "markdown", "id": "4c88de77", "metadata": {}, "source": [ "## 7. Reassemble per-tile results\n", "\n", "After running each tile through the Risk Analyzer, drop the output GeoJSONs into the\n", "results folder (`RESULTS_DIR`). This section loads every result file, tags each\n", "feature with its source tile, concatenates them, and maps the combined result.\n", "\n", "- If the Risk Analyzer returns **polygons with a risk attribute**, set `RESULT_VALUE_FIELD`\n", " to that attribute name to get a choropleth.\n", "- If results are just geometries (or you haven't run any tiles yet), the cell falls\n", " back to a demo mode so you can see how the reassembly will look.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "95989d2a", "metadata": {}, "outputs": [], "source": [ "RESULT_VALUE_FIELD = None # e.g. \"risk_score\" once you know the output schema\n", "\n", "result_files = sorted(RESULTS_DIR.glob(\"*.geojson\"))\n", "print(f\"Found {len(result_files)} result file(s) in {RESULTS_DIR}/\")\n", "\n", "if result_files:\n", " parts = []\n", " for f in result_files:\n", " g = gpd.read_file(f)\n", " if g.crs is None:\n", " g = g.set_crs(\"EPSG:4326\")\n", " g = g.to_crs(\"EPSG:4326\")\n", " # Recover the tile_id from the filename (matches the naming convention)\n", " g[\"source_tile\"] = f.stem\n", " parts.append(g)\n", " combined = gpd.GeoDataFrame(pd.concat(parts, ignore_index=True), crs=\"EPSG:4326\")\n", " print(f\"Combined result: {len(combined)} features from {len(result_files)} tiles\")\n", "else:\n", " # ---- DEMO MODE: pretend each tile came back with a risk score ----------\n", " import numpy as np\n", " rng = np.random.default_rng(42)\n", " combined = tiles_wgs.copy()\n", " combined[\"source_tile\"] = combined[\"tile_id\"]\n", " combined[\"risk_score\"] = rng.uniform(0, 1, len(combined)).round(3)\n", " RESULT_VALUE_FIELD = \"risk_score\"\n", " print(\"DEMO MODE: no results found yet - showing synthetic per-tile risk scores\")\n", "\n", "# Optional: merge/dissolve analysis coverage back into one footprint\n", "merged_footprint = unary_union(combined.geometry.values)\n", "print(f\"Merged analysis footprint: {gpd.GeoSeries([merged_footprint], crs='EPSG:4326').to_crs(GRID_CRS).area.iloc[0]/1e6:,.0f} km2\")" ] }, { "cell_type": "code", "execution_count": null, "id": "9f90950f", "metadata": {}, "outputs": [], "source": [ "m2 = folium.Map(location=center, zoom_start=9, tiles=\"cartodbpositron\")\n", "\n", "folium.GeoJson(\n", " aoi.to_crs(\"EPSG:4326\"),\n", " name=\"Original AOI\",\n", " style_function=lambda f: {\"color\": \"#333\", \"weight\": 3, \"fill\": False, \"dashArray\": \"6\"},\n", ").add_to(m2)\n", "\n", "if RESULT_VALUE_FIELD and RESULT_VALUE_FIELD in combined.columns:\n", " combined.explore(\n", " m=m2,\n", " column=RESULT_VALUE_FIELD,\n", " cmap=\"YlOrRd\",\n", " name=\"Analysis results\",\n", " tooltip=[\"source_tile\", RESULT_VALUE_FIELD],\n", " style_kwds={\"weight\": 1, \"color\": \"#888\"},\n", " legend=True,\n", " )\n", "else:\n", " folium.GeoJson(\n", " combined,\n", " name=\"Analysis results\",\n", " style_function=lambda f: {\"color\": \"#b91c1c\", \"weight\": 1,\n", " \"fillColor\": \"#b91c1c\", \"fillOpacity\": 0.25},\n", " tooltip=folium.GeoJsonTooltip(fields=[\"source_tile\"], aliases=[\"Tile\"]),\n", " ).add_to(m2)\n", "\n", "folium.LayerControl().add_to(m2)\n", "m2" ] }, { "cell_type": "code", "execution_count": null, "id": "e9f0c19d", "metadata": {}, "outputs": [], "source": [ "# Export the reassembled result as a single GeoJSON\n", "merged_path = f\"{AOI_NAME}_{TILE_SIZE_KM}km_combined_results.geojson\"\n", "combined.to_file(merged_path, driver=\"GeoJSON\")\n", "print(f\"Wrote combined results to {merged_path}\")" ] }, { "cell_type": "markdown", "id": "0a7e521f", "metadata": {}, "source": [ "---\n", "### Notes & tips\n", "\n", "- **Alignment across runs:** keep `GRID_CRS`, `GRID_ORIGIN`, and `TILE_SIZE_KM`\n", " constant across every AOI you cut and the tiles will all share one global grid.\n", " Halving the tile size keeps the smaller tiles nested inside the larger ones.\n", "- **Tile IDs are global coordinates:** `c-0012_r+0045` means grid column −12, row 45\n", " on the shared grid — the same ID always maps to the same square of ground.\n", "- **Rectangular vs clipped tiles:** this notebook clips tiles to the AOI boundary so\n", " no data outside the county is requested. If the Risk Analyzer prefers rectangles,\n", " replace `cell.intersection(aoi_geom)` with `cell` in Section 3 (tiles will then\n", " overlap the county line, and the coverage check in Section 4 will report the excess).\n", "- **If a specific tile still fails** (NASA download too large), re-run just that area\n", " at a smaller `TILE_SIZE_KM` — nesting means the sub-tiles will fit exactly inside it.\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }