State County Boundaries.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73b6ecf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utz import *\n",
    "from IPython.display import HTML, SVG, display\n",
    "\n",
    "import fiona\n",
    "import geopandas as gpd\n",
    "from geopy import distance\n",
    "import shapefile\n",
    "from shapely.geometry import Point\n",
    "from shapely.ops import cascaded_union, unary_union\n",
    "\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "import sys\n",
    "sys.executable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0d580b",
   "metadata": {},
   "source": [
    "## Mapbox configs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1a395ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "token = open(\".mapbox-token\").read()\n",
    "mapbox = dict(mapbox=dict(\n",
    "    style=\"dark\",\n",
    "    accesstoken=token,\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bf674d3",
   "metadata": {},
   "source": [
    "## Load county shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "237b8fc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fiona.listlayers(\"tlgdb_2022_a_34_nj.gdb/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64657bd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fiona.listlayers(\"tlgdb_2022_a_34_nj_edges.gdb/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aab203b",
   "metadata": {},
   "outputs": [],
   "source": [
    "counties = gpd.read_file(\"tlgdb_2022_a_34_nj.gdb/\", layer='County')\n",
    "counties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e155ce16",
   "metadata": {},
   "outputs": [],
   "source": [
    "cgeos = counties.geometry.tolist()\n",
    "cgeos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d6956bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "[ len(c.boundary.geoms) for c in cgeos ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a5614d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "counties_union = unary_union(counties.geometry.tolist())\n",
    "counties_union"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6d06fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(counties_union)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ba9204b",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(counties_union.boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95c68e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "nj_boundary = counties_union.boundary\n",
    "nj_boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9c65e5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "nj_bnd_lines = list(nj_boundary.geoms)\n",
    "nj_bnd_lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7edc7edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "[ display(ln) for ln in nj_bnd_lines ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f7b682d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln = nj_bnd_lines[0]\n",
    "ln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f13061bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0835417e",
   "metadata": {},
   "outputs": [],
   "source": [
    "[ ln.length for ln in nj_bnd_lines ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4740248",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(nj_boundary.geoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "643d8e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(counties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a38a45b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_nj_ll(lat, lon):\n",
    "    point = Point(lon, lat)\n",
    "    return counties_union.contains(point)\n",
    "\n",
    "def is_nj(r):\n",
    "    return is_nj_ll(r.LAT, r.LON)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baaf5e1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_county(r):\n",
    "    hits = counties[counties.geometry.contains(Point(r.LON, r.LAT))].NAMELSAD\n",
    "    if len(hits) > 1:\n",
    "        err(f\"{r}: {len(hits)} counties: {hits}\")\n",
    "        return\n",
    "    elif hits.empty:\n",
    "        return\n",
    "    [county] = hits.tolist()\n",
    "    return county"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d226fd43",
   "metadata": {},
   "source": [
    "## County / State Boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b9549fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = counties.geometry.boundary.to_list()\n",
    "bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40a6d64f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "[len(b.geoms[0].coords) for b in bounds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf125f34",
   "metadata": {},
   "outputs": [],
   "source": [
    "counties.geometry.boundary.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b350ddd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def county_points(r):\n",
    "    [name] = r.index.unique()\n",
    "    [multilinestr] = r.geometry.boundary.tolist()\n",
    "    linestrs = multilinestr.geoms\n",
    "    if len(linestrs) > 1:\n",
    "        err(f'{name}: {len(linestrs)} linestrings')\n",
    "    return Series([ c for l in linestrs for c in l.coords ], name='point')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd79c6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "county_coords = (\n",
    "    counties\n",
    "    .rename(columns={'NAMELSAD': 'name'})\n",
    "    .groupby('name')\n",
    "    .apply(county_points)\n",
    "    .reset_index(level=1, drop=True)\n",
    "    .reset_index()\n",
    ")\n",
    "county_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "195e513b",
   "metadata": {},
   "outputs": [],
   "source": [
    "county_coords.point.value_counts().value_counts().sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "438f0621",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "bnd_lls = sxs(county_coords.name, county_coords.point.apply(lambda p: Series(p, index=['lon', 'lat'])))\n",
    "bnd_lls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3a1de23",
   "metadata": {},
   "outputs": [],
   "source": [
    "ll_hist = bnd_lls[['lat','lon']].value_counts()\n",
    "ll1s = ll_hist[ll_hist == 1]\n",
    "ll_hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b5b3d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ll1s.index.to_series()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5b97def",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "p1s = bnd_lls.merge(ll1s, left_on=['lat', 'lon'], right_index=True).drop(columns='count')\n",
    "p1s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4347e2ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter_mapbox(\n",
    "    p1s,\n",
    "    lat=\"lat\", lon=\"lon\",\n",
    "    color='name',\n",
    "    #hover_name=\"City\",\n",
    "    hover_data=[\"name\"],\n",
    "    #color_discrete_sequence=[\"yellow\", \"orange\", \"red\"],\n",
    "    center=dict(lat=40.15, lon=-74.715),\n",
    "    zoom=7.6,\n",
    "    height=1000,\n",
    ")\n",
    "legend_bgcolor = '50'\n",
    "fig.update_layout(\n",
    "    **mapbox,\n",
    "    title=dict(\n",
    "        text=\"NJ boundary points, by county\",\n",
    "        x=0.5, y=0.98,\n",
    "        xanchor='center', yanchor='top',\n",
    "        font=dict(size=32, color=\"white\")\n",
    "    ),\n",
    "    legend=dict(\n",
    "        title=dict(text=''),\n",
    "        x=0.98, y=0.98,\n",
    "        xanchor=\"right\", yanchor=\"top\",\n",
    "        font=dict(\n",
    "            size=14,\n",
    "            color=\"white\"\n",
    "        ),\n",
    "        bgcolor=f\"rgba({legend_bgcolor},{legend_bgcolor},{legend_bgcolor},0.8)\",\n",
    "        bordercolor=\"white\",\n",
    "        borderwidth=2,\n",
    "    ),\n",
    "    margin={\"r\":0,\"t\":0,\"l\":0,\"b\":0},\n",
    ")\n",
    "fig.write_image('nj-boundary-by-county.png', width=1000, height=1000)\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda-geopandas",
   "language": "python",
   "name": "conda-geopandas"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}