diff --git a/analyses/tree_density_methodology.ipynb b/analyses/tree_density_methodology.ipynb index 424e220..e07daf7 100644 --- a/analyses/tree_density_methodology.ipynb +++ b/analyses/tree_density_methodology.ipynb @@ -5,9 +5,11 @@ "id": "title", "metadata": {}, "source": [ - "# Street tree density visual check\n", + "# Polygon-intersection tree density prototype\n", "\n", - "This notebook picks three postcodes with a fixed random seed and shows only the final map plus the tree-density percentage for each postcode. The percentage is read from the production output generated by `pipeline.transform.tree_density`; the map uses the same TOW type filter, postcode centroid, 50m buffer, and TOW polygon centroid rule used by that pipeline." + "This notebook prototypes the next algorithm: expand each postcode boundary by `N` metres, intersect that expanded polygon with Forest Research Trees Outside Woodland polygons, and sum the actual intersecting canopy area.\n", + "\n", + "The map shows three seeded random postcodes side by side. Each panel shows the original postcode boundary, the expanded postcode boundary, the intersecting foliage polygons, and the resulting percentage." ] }, { @@ -16,10 +18,10 @@ "id": "imports-config", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T20:42:01.528837Z", - "iopub.status.busy": "2026-05-12T20:42:01.528729Z", - "iopub.status.idle": "2026-05-12T20:42:01.986980Z", - "shell.execute_reply": "2026-05-12T20:42:01.986661Z" + "iopub.execute_input": "2026-05-12T20:47:25.026395Z", + "iopub.status.busy": "2026-05-12T20:47:25.026327Z", + "iopub.status.idle": "2026-05-12T20:47:26.779075Z", + "shell.execute_reply": "2026-05-12T20:47:26.778675Z" } }, "outputs": [], @@ -39,7 +41,7 @@ "import shapely\n", "from IPython.display import HTML\n", "from pyproj import Transformer\n", - "from shapely.geometry import Point, mapping, shape\n", + "from shapely.geometry import mapping, shape\n", "from shapely.ops import transform as shapely_transform\n", "\n", "ROOT = Path.cwd()\n", @@ -64,15 +66,15 @@ "ARCGIS = DATA_DIR / \"arcgis_data.parquet\"\n", "POSTCODE_TREE_DENSITY = DATA_DIR / \"tree_density_by_postcode.parquet\"\n", "POSTCODE_BOUNDARY_UNITS = DATA_DIR / \"postcode_boundaries\" / \"units\"\n", - "OUTPUT_HTML = ROOT / \"analyses\" / \"tree_density_three_postcodes_map.html\"\n", + "OUTPUT_HTML = ROOT / \"analyses\" / \"tree_density_polygon_intersection_map.html\"\n", "\n", "RANDOM_SEED = 20260512\n", "N_POSTCODES = 3\n", - "RADIUS_M = 50\n", + "BOUNDARY_EXPANSION_M = 50\n", + "CANDIDATE_SAMPLE_SIZE = 60\n", "\n", "tow_types = _parse_csv_arg(\",\".join(DEFAULT_TOW_TYPES))\n", - "density_col, area_col, count_col, height_col = _metric_columns(RADIUS_M)\n", - "buffer_area = math.pi * RADIUS_M**2\n" + "density_col, area_col, count_col, height_col = _metric_columns(50)\n" ] }, { @@ -81,18 +83,23 @@ "id": "helpers", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T20:42:01.988203Z", - "iopub.status.busy": "2026-05-12T20:42:01.988094Z", - "iopub.status.idle": "2026-05-12T20:42:02.015890Z", - "shell.execute_reply": "2026-05-12T20:42:02.015642Z" + "iopub.execute_input": "2026-05-12T20:47:26.780431Z", + "iopub.status.busy": "2026-05-12T20:47:26.780263Z", + "iopub.status.idle": "2026-05-12T20:47:26.846177Z", + "shell.execute_reply": "2026-05-12T20:47:26.845648Z" } }, "outputs": [], "source": [ + "_to_bng = Transformer.from_crs(\"EPSG:4326\", \"EPSG:27700\", always_xy=True)\n", "_to_wgs84 = Transformer.from_crs(\"EPSG:27700\", \"EPSG:4326\", always_xy=True)\n", "\n", "\n", - "def bng_geom_to_wgs84(geom):\n", + "def to_bng(geom):\n", + " return shapely_transform(_to_bng.transform, geom)\n", + "\n", + "\n", + "def to_wgs84(geom):\n", " return shapely_transform(_to_wgs84.transform, geom)\n", "\n", "\n", @@ -111,8 +118,15 @@ " raise ValueError(f\"Postcode boundary not found for {postcode}\")\n", "\n", "\n", - "def tow_features_for_postcode(dataset_path: str, layer_names: list[str], x: float, y: float):\n", - " bbox = (x - RADIUS_M, y - RADIUS_M, x + RADIUS_M, y + RADIUS_M)\n", + "def postcode_boundary_bng(postcode: str):\n", + " geom = to_bng(postcode_boundary_wgs84(postcode))\n", + " if not geom.is_valid:\n", + " geom = shapely.make_valid(geom)\n", + " return geom\n", + "\n", + "\n", + "def intersecting_tow_features(dataset_path: str, layer_names: list[str], target_bng):\n", + " bbox = target_bng.bounds\n", " rows = []\n", " where = _where_for_tow_types(tow_types)\n", "\n", @@ -128,11 +142,11 @@ " ) as (_meta, reader):\n", " for batch in reader:\n", " names = batch.schema.names\n", - " area = np.asarray(\n", + " tow_area = np.asarray(\n", " batch.column(names.index(\"TOW_Area_M\")).to_numpy(zero_copy_only=False),\n", " dtype=np.float64,\n", " )\n", - " height = np.asarray(\n", + " mean_height = np.asarray(\n", " batch.column(names.index(\"MEANHT\")).to_numpy(zero_copy_only=False),\n", " dtype=np.float64,\n", " )\n", @@ -143,41 +157,49 @@ " )\n", "\n", " geoms = shapely.from_wkb(geometry)\n", - " centroids = shapely.centroid(geoms)\n", - " centroid_x = shapely.get_x(centroids)\n", - " centroid_y = shapely.get_y(centroids)\n", - " valid = np.isfinite(area) & (area > 0) & np.isfinite(centroid_x) & np.isfinite(centroid_y)\n", - "\n", + " valid = np.isfinite(tow_area) & (tow_area > 0)\n", " for i in np.flatnonzero(valid):\n", - " distance_m = float(np.hypot(centroid_x[i] - x, centroid_y[i] - y))\n", + " geom = geoms[i]\n", + " if geom is None or shapely.is_empty(geom):\n", + " continue\n", + " if not shapely.is_valid(geom):\n", + " geom = shapely.make_valid(geom)\n", + " if not shapely.intersects(geom, target_bng):\n", + " continue\n", + "\n", + " intersection = shapely.intersection(geom, target_bng)\n", + " intersection_area = float(shapely.area(intersection))\n", + " if intersection_area <= 0:\n", + " continue\n", + "\n", " rows.append(\n", " {\n", " \"layer\": layer,\n", - " \"geometry_bng\": geoms[i],\n", - " \"centroid_bng\": Point(float(centroid_x[i]), float(centroid_y[i])),\n", + " \"geometry_bng\": geom,\n", + " \"intersection_bng\": intersection,\n", " \"woodland_type\": woodland_type[i],\n", - " \"area_m2\": float(area[i]),\n", - " \"mean_height_m\": None if not np.isfinite(height[i]) else float(height[i]),\n", - " \"distance_m\": distance_m,\n", - " \"counted\": distance_m <= RADIUS_M,\n", + " \"tow_area_m2\": float(tow_area[i]),\n", + " \"intersection_area_m2\": intersection_area,\n", + " \"mean_height_m\": None if not np.isfinite(mean_height[i]) else float(mean_height[i]),\n", " }\n", " )\n", " return rows\n", "\n", "\n", - "def build_postcode_map(dataset_path: str, layer_names: list[str], row: dict) -> str:\n", - " postcode = row[\"postcode\"]\n", - " point_bng = Point(row[\"x\"], row[\"y\"])\n", - " point_lon, point_lat = _to_wgs84.transform(point_bng.x, point_bng.y)\n", - " buffer_bng = point_bng.buffer(RADIUS_M, resolution=96)\n", - " boundary = postcode_boundary_wgs84(postcode)\n", - " features = tow_features_for_postcode(dataset_path, layer_names, point_bng.x, point_bng.y)\n", + "def build_postcode_card(dataset_path: str, layer_names: list[str], postcode: str) -> tuple[str, dict]:\n", + " original_bng = postcode_boundary_bng(postcode)\n", + " expanded_bng = original_bng.buffer(BOUNDARY_EXPANSION_M)\n", + " if not expanded_bng.is_valid:\n", + " expanded_bng = shapely.make_valid(expanded_bng)\n", "\n", - " counted_area = sum(feature[\"area_m2\"] for feature in features if feature[\"counted\"])\n", - " visual_density = round(min(counted_area / buffer_area * 100, 100), 1)\n", + " features = intersecting_tow_features(dataset_path, layer_names, expanded_bng)\n", + " intersecting_area = sum(feature[\"intersection_area_m2\"] for feature in features)\n", + " expanded_area = float(shapely.area(expanded_bng))\n", + " polygon_density = 0.0 if expanded_area <= 0 else min(intersecting_area / expanded_area * 100, 100)\n", "\n", + " center_lon, center_lat = to_wgs84(expanded_bng.centroid).coords[0]\n", " m = folium.Map(\n", - " location=[point_lat, point_lon],\n", + " location=[center_lat, center_lon],\n", " zoom_start=18,\n", " tiles=\"CartoDB positron\",\n", " control_scale=True,\n", @@ -186,18 +208,8 @@ " )\n", "\n", " folium.GeoJson(\n", - " mapping(boundary),\n", - " name=\"postcode boundary\",\n", - " style_function=lambda _feature: {\n", - " \"color\": \"#2563eb\",\n", - " \"weight\": 3,\n", - " \"fillColor\": \"#93c5fd\",\n", - " \"fillOpacity\": 0.10,\n", - " },\n", - " ).add_to(m)\n", - " folium.GeoJson(\n", - " mapping(bng_geom_to_wgs84(buffer_bng)),\n", - " name=\"50m buffer\",\n", + " mapping(to_wgs84(expanded_bng)),\n", + " name=f\"expanded postcode boundary (+{BOUNDARY_EXPANSION_M}m)\",\n", " style_function=lambda _feature: {\n", " \"color\": \"#f97316\",\n", " \"weight\": 3,\n", @@ -205,70 +217,73 @@ " \"fillOpacity\": 0.18,\n", " },\n", " ).add_to(m)\n", + " folium.GeoJson(\n", + " mapping(to_wgs84(original_bng)),\n", + " name=\"original postcode boundary\",\n", + " style_function=lambda _feature: {\n", + " \"color\": \"#2563eb\",\n", + " \"weight\": 3,\n", + " \"fillColor\": \"#93c5fd\",\n", + " \"fillOpacity\": 0.12,\n", + " },\n", + " ).add_to(m)\n", "\n", - " counted_group = folium.FeatureGroup(name=\"counted foliage\", show=True)\n", - " nearby_group = folium.FeatureGroup(name=\"nearby, not counted\", show=False)\n", - "\n", + " foliage_group = folium.FeatureGroup(name=\"foliage boundary\", show=True)\n", + " intersection_group = folium.FeatureGroup(name=\"intersection area\", show=True)\n", " for index, feature in enumerate(features, start=1):\n", - " group = counted_group if feature[\"counted\"] else nearby_group\n", - " geom_wgs84 = bng_geom_to_wgs84(feature[\"geometry_bng\"])\n", - " centroid_lon, centroid_lat = _to_wgs84.transform(feature[\"centroid_bng\"].x, feature[\"centroid_bng\"].y)\n", - " style = {\n", - " \"color\": \"#15803d\" if feature[\"counted\"] else \"#6b7280\",\n", - " \"weight\": 2 if feature[\"counted\"] else 1,\n", - " \"fillColor\": \"#22c55e\" if feature[\"counted\"] else \"#9ca3af\",\n", - " \"fillOpacity\": 0.45 if feature[\"counted\"] else 0.20,\n", - " }\n", " popup_html = (\n", " f\"TOW polygon {index}
\"\n", - " f\"Status: {'counted' if feature['counted'] else 'not counted'}
\"\n", " f\"Type: {feature['woodland_type']}
\"\n", - " f\"Area: {feature['area_m2']:.1f} sqm
\"\n", - " f\"Centroid distance: {feature['distance_m']:.1f} m\"\n", + " f\"Full TOW area: {feature['tow_area_m2']:.1f} sqm
\"\n", + " f\"Intersecting area: {feature['intersection_area_m2']:.1f} sqm\"\n", " )\n", " folium.GeoJson(\n", - " {\n", - " \"type\": \"Feature\",\n", - " \"geometry\": mapping(geom_wgs84),\n", - " \"properties\": {\"popup\": popup_html},\n", + " mapping(to_wgs84(feature[\"geometry_bng\"])),\n", + " name=\"foliage boundary\",\n", + " style_function=lambda _feature: {\n", + " \"color\": \"#166534\",\n", + " \"weight\": 1,\n", + " \"fillColor\": \"#22c55e\",\n", + " \"fillOpacity\": 0.16,\n", " },\n", - " style_function=lambda _feature, style=style: style,\n", - " tooltip=\"counted foliage\" if feature[\"counted\"] else \"nearby foliage, not counted\",\n", - " popup=folium.Popup(popup_html, max_width=280),\n", - " ).add_to(group)\n", - " folium.CircleMarker(\n", - " [centroid_lat, centroid_lon],\n", - " radius=3,\n", - " color=\"#14532d\" if feature[\"counted\"] else \"#4b5563\",\n", - " fill=True,\n", - " fill_opacity=0.95,\n", - " tooltip=f\"TOW centroid: {feature['distance_m']:.1f}m\",\n", - " ).add_to(group)\n", + " tooltip=\"foliage boundary\",\n", + " popup=folium.Popup(popup_html, max_width=300),\n", + " ).add_to(foliage_group)\n", + " folium.GeoJson(\n", + " mapping(to_wgs84(feature[\"intersection_bng\"])),\n", + " name=\"intersection area\",\n", + " style_function=lambda _feature: {\n", + " \"color\": \"#14532d\",\n", + " \"weight\": 2,\n", + " \"fillColor\": \"#16a34a\",\n", + " \"fillOpacity\": 0.58,\n", + " },\n", + " tooltip=\"counted intersection area\",\n", + " popup=folium.Popup(popup_html, max_width=300),\n", + " ).add_to(intersection_group)\n", "\n", - " counted_group.add_to(m)\n", - " nearby_group.add_to(m)\n", - " folium.CircleMarker(\n", - " [point_lat, point_lon],\n", - " radius=7,\n", - " color=\"#1d4ed8\",\n", - " fill=True,\n", - " fill_color=\"#1d4ed8\",\n", - " fill_opacity=1,\n", - " tooltip=f\"{postcode} centroid\",\n", - " ).add_to(m)\n", + " foliage_group.add_to(m)\n", + " intersection_group.add_to(m)\n", " folium.LayerControl(collapsed=True).add_to(m)\n", "\n", " title = escape(\n", - " f\"{postcode}: {row[density_col]:.1f}% | {int(row[count_col])} features | \"\n", - " f\"visual check {visual_density:.1f}%\"\n", + " f\"{postcode}: {polygon_density:.1f}% | \"\n", + " f\"intersection {intersecting_area:.0f} sqm | +{BOUNDARY_EXPANSION_M}m\"\n", " )\n", " iframe = escape(m.get_root().render(), quote=True)\n", - " return f\"\"\"\n", + " card = f\"\"\"\n", "
\n", "

{title}

\n", - " \n", + " \n", "
\n", " \"\"\"\n", + " return card, {\n", + " \"postcode\": postcode,\n", + " \"polygon_density_pct\": round(polygon_density, 1),\n", + " \"intersecting_area_sqm\": round(intersecting_area, 1),\n", + " \"expanded_boundary_area_sqm\": round(expanded_area, 1),\n", + " \"intersecting_tow_features\": len(features),\n", + " }\n", "\n", "\n", "def side_by_side_document(cards: list[str]) -> str:\n", @@ -278,7 +293,7 @@ " \n", " \n", " \n", - " Street tree density visual check\n", + " Polygon-intersection tree density prototype\n", " \n", + "shape: (3, 5)
postcodepolygon_density_pctintersecting_area_sqmexpanded_boundary_area_sqmintersecting_tow_features
strf64f64f64i64
"B61 7EF"3.720421.5554408.7182
"B94 6HR"4.32281.352697.133
"BB1 9JJ"11.52888.425147.639
" + ], + "text/plain": [ + "shape: (3, 5)\n", + "┌──────────┬─────────────────────┬─────────────────────┬─────────────────────┬─────────────────────┐\n", + "│ postcode ┆ polygon_density_pct ┆ intersecting_area_s ┆ expanded_boundary_a ┆ intersecting_tow_fe │\n", + "│ --- ┆ --- ┆ qm ┆ rea_sqm ┆ atures │\n", + "│ str ┆ f64 ┆ --- ┆ --- ┆ --- │\n", + "│ ┆ ┆ f64 ┆ f64 ┆ i64 │\n", + "╞══════════╪═════════════════════╪═════════════════════╪═════════════════════╪═════════════════════╡\n", + "│ B61 7EF ┆ 3.7 ┆ 20421.5 ┆ 554408.7 ┆ 182 │\n", + "│ B94 6HR ┆ 4.3 ┆ 2281.3 ┆ 52697.1 ┆ 33 │\n", + "│ BB1 9JJ ┆ 11.5 ┆ 2888.4 ┆ 25147.6 ┆ 39 │\n", + "└──────────┴─────────────────────┴─────────────────────┴─────────────────────┴─────────────────────┘" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, { "data": { "text/html": [ @@ -318,7 +362,7 @@ " \n", " \n", " \n", - " Street tree density visual check\n", + " Polygon-intersection tree density prototype\n", "