{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use OSMnx to plot street network over place shape\n",
    "\n",
    "Author: [Geoff Boeing](https://geoffboeing.com/)\n",
    "\n",
    "  - [Documentation](https://osmnx.readthedocs.io/)\n",
    "  - [Journal article and citation info](https://doi.org/10.1111/gean.70009)\n",
    "  - [Code repository](https://github.com/gboeing/osmnx)\n",
    "  - [Examples gallery](https://github.com/gboeing/osmnx-examples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!uv pip install --system --quiet osmnx[all]\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import osmnx as ox\n",
    "\n",
    "ox.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example uses Portland, Maine - a city with several islands within its municipal boundaries. Thus, we set `retain_all=True` when getting the network so that we keep all the graph components, not just the largest connected component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the place boundaries\n",
    "place = \"Portland, Maine\"\n",
    "gdf = ox.geocoder.geocode_to_gdf(place)\n",
    "\n",
    "# get the street network, with retain_all=True to retain all the disconnected islands' networks\n",
    "G = ox.graph.graph_from_place(place, network_type=\"drive\", retain_all=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the network, but do not show it or close it yet\n",
    "fig, ax = ox.plot.plot_graph(\n",
    "    G,\n",
    "    show=False,\n",
    "    close=False,\n",
    "    bgcolor=\"#111111\",\n",
    "    edge_color=\"#ffcb00\",\n",
    "    edge_linewidth=0.3,\n",
    "    node_size=0,\n",
    ")\n",
    "\n",
    "# to this matplotlib axis, add the place shape(s)\n",
    "gdf.plot(ax=ax, fc=\"#444444\", ec=None, lw=1, alpha=1, zorder=-1)\n",
    "\n",
    "# optionally set up the axes extents\n",
    "margin = 0.02\n",
    "west, south, east, north = gdf.union_all().bounds\n",
    "margin_ns = (north - south) * margin\n",
    "margin_ew = (east - west) * margin\n",
    "ax.set_ylim((south - margin_ns, north + margin_ns))\n",
    "ax.set_xlim((west - margin_ew, east + margin_ew))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice this municipal boundary is an administrative boundary, not a physical boundary, so it represents jurisdictional bounds, not individual physical features like islands. To retrieve individual islands' geometries, use the `features` module to search for features matching certain OSM tags:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "islands = ox.features.features_from_place(place, tags={\"place\": [\"island\", \"islet\"]})\n",
    "islands.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do you remove the water and just get a coastal city's land?\n",
    "\n",
    "OpenStreetMap maps coastlines rather than water, so the task isn't as straightforward as \"download the political boundary, then subtract the water from it.\" Instead, the workflow is to get the political boundary and coastlines as linestrings, polygonize them, then drop the resulting water polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all municipalities and coastlines within Salvador\n",
    "place = \"Salvador, Bahía, Brasil\"\n",
    "tags = {\"natural\": \"coastline\", \"place\": \"municipality\"}\n",
    "gdf = ox.features.features_from_place(place, tags)\n",
    "\n",
    "# project and filter results to keep only Salvador municipality and coastlines\n",
    "mask1 = gdf[\"name\"] == \"Salvador\"\n",
    "mask2 = gdf[\"natural\"] == \"coastline\"\n",
    "gdf = ox.projection.project_gdf(gdf.loc[mask1 | mask2])\n",
    "\n",
    "# replace Salvador municipal boundary polygon with boundary linestring\n",
    "gdf.loc[mask1, \"geometry\"] = gdf.loc[mask1, \"geometry\"].boundary\n",
    "\n",
    "# polygonize the linestrings, drop the largest polygon (open water)\n",
    "gdf_land = gpd.GeoDataFrame(geometry=gdf.polygonize(), crs=gdf.crs)\n",
    "gdf_land = gdf_land.drop(gdf_land.area.idxmax())\n",
    "\n",
    "# plot it\n",
    "fig, ax = plt.subplots(facecolor=\"#111111\")\n",
    "ax = gdf_land.plot(ax=ax, color=\"#ffcb00\")\n",
    "_ = ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}