{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 6 — Programmatic Spatial Data Analysis\n",
    "\n",
    "**Companion notebook** for the Lesson 6 slides. Uses `../data/geo/nc_geo.gpkg` and\n",
    "the ABT from Lesson 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check(label, got, expected, tol=1e-6):\n",
    "    \"\"\"Compare your answer with the expected value; prints ✅ or ❌ (never raises).\"\"\"\n",
    "    import math\n",
    "    try:\n",
    "        if isinstance(expected, float) or isinstance(got, float):\n",
    "            ok = math.isclose(float(got), float(expected), rel_tol=tol, abs_tol=tol)\n",
    "        elif isinstance(expected, (list, tuple)) and expected and isinstance(expected[0], str):\n",
    "            ok = list(got) == list(expected)\n",
    "        elif hasattr(expected, \"__len__\") and not isinstance(expected, str):\n",
    "            import numpy as np\n",
    "            ok = np.allclose(np.asarray(got, dtype=float), np.asarray(expected, dtype=float), rtol=tol, atol=tol)\n",
    "        else:\n",
    "            ok = got == expected\n",
    "    except Exception:\n",
    "        ok = False\n",
    "    print((\"✅ \" if ok else \"❌ \") + f\"{label}: got {got!r}, expected {expected!r}\")\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyogrio\n",
    "from pyproj import Geod\n",
    "\n",
    "GPKG = \"../data/geo/nc_geo.gpkg\"\n",
    "pyogrio.list_layers(GPKG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-002",
   "metadata": {},
   "outputs": [],
   "source": [
    "counties = gpd.read_file(GPKG, layer=\"counties\")\n",
    "places = gpd.read_file(GPKG, layer=\"places\")\n",
    "coast = gpd.read_file(GPKG, layer=\"coastline\")\n",
    "print(counties.crs, len(counties), len(places), len(coast))\n",
    "counties.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-003",
   "metadata": {},
   "source": [
    "## 1. Coordinate reference systems: measure in a projected CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-004",
   "metadata": {},
   "outputs": [],
   "source": [
    "cities = gpd.GeoDataFrame({\"city\": [\"Raleigh\", \"Wilmington\"]},\n",
    "                          geometry=gpd.points_from_xy([-78.6382, -77.9447], [35.7796, 34.2257]), crs=\"EPSG:4326\")\n",
    "print(\"degrees:\", round(cities.distance(cities.iloc[[1]].geometry.iloc[0]).iloc[0], 4))\n",
    "for epsg in (3857, 32119):\n",
    "    c = cities.to_crs(epsg=epsg)\n",
    "    print(f\"EPSG:{epsg}: {c.geometry.iloc[0].distance(c.geometry.iloc[1]) / 1000:,.2f} km\")\n",
    "print(\"geodesic:\", round(Geod(ellps=\"WGS84\").inv(-78.6382, 35.7796, -77.9447, 34.2257)[2] / 1000, 2), \"km\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-005",
   "metadata": {},
   "source": [
    "## 2. Points from latitude/longitude (longitude is x!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-006",
   "metadata": {},
   "outputs": [],
   "source": [
    "abt = pd.read_parquet(\"../data/processed/abt.parquet\")\n",
    "pts = gpd.GeoDataFrame(abt, geometry=gpd.points_from_xy(abt[\"longitude\"], abt[\"latitude\"]),\n",
    "                       crs=\"EPSG:4326\").to_crs(counties.crs)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(11, 4.5))\n",
    "counties.plot(ax=ax, color=\"#f3f0e8\", edgecolor=\"#c9c2b1\", linewidth=0.4)\n",
    "places.plot(ax=ax, color=\"#0e6b6b\", alpha=0.35)\n",
    "coast.plot(ax=ax, color=\"#1b2330\", linewidth=0.8)\n",
    "pts.sample(3000, random_state=1).plot(ax=ax, markersize=1, color=\"#c26a1b\")\n",
    "ax.set_axis_off()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-007",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — which county is each vehicle in?\n",
    "\n",
    "Spatial-join `pts` to `counties` with `predicate=\"within\"` (keep `county_fips` and\n",
    "`county_name`). Count the rows and the unmatched vehicles, then rescue the misses with\n",
    "`gpd.sjoin_nearest(..., max_distance=2_000)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-008",
   "metadata": {},
   "outputs": [],
   "source": [
    "j = None\n",
    "n_unmatched = n_rescued = None\n",
    "...  # your code here\n",
    "check(\"rows after sjoin\", None if j is None else len(j), 23993)\n",
    "check(\"unmatched\", n_unmatched, 12)\n",
    "check(\"rescued by sjoin_nearest\", n_rescued, 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-009",
   "metadata": {},
   "source": [
    "## Your turn 2 — urban flag and distance to the coast\n",
    "\n",
    "* `urban`: is the vehicle inside any incorporated place?\n",
    "* `coast_km`: distance to the coastline (union all pieces first), in km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-010",
   "metadata": {},
   "outputs": [],
   "source": [
    "if j is None:\n",
    "    j = pts.copy()\n",
    "urban_share = median_coast_km = None\n",
    "...  # your code here\n",
    "check(\"urban share\", None if urban_share is None else round(urban_share, 4), 0.5943)\n",
    "check(\"median distance to coast (km)\", None if median_coast_km is None else round(median_coast_km, 1), 164.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-011",
   "metadata": {},
   "source": [
    "## Your turn 3 — territories and relativities\n",
    "\n",
    "Territory = **Coastal** if within 25 km of the coast, else **Urban inland** if urban,\n",
    "else **Rural inland**. Relativity = the territory's pure premium ÷ the portfolio's."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-012",
   "metadata": {},
   "outputs": [],
   "source": [
    "relativity = None\n",
    "if \"coast_km\" in j:\n",
    "    ...  # your code here\n",
    "check(\"Coastal relativity\", None if relativity is None else round(relativity[\"Coastal\"], 3), 1.073)\n",
    "check(\"Rural inland relativity\", None if relativity is None else round(relativity[\"Rural inland\"], 3), 0.870)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}