script.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computation of isochrones with a local instance of OpenRouteService using online OSM data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation libraries pyhton\n",
    "\n",
    "```bash\n",
    "conda env create -f environment.yml \n",
    "conda activate ORS_stations\n",
    "pip install -r requirements.txt\n",
    "wget https://gist.githubusercontent.com/perrygeo/5667173/raw/763e1e50208e8c853f46e57cd07bb07b424fed10/zonal_stats.py\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "from IPython.display import display\n",
    "from ipywidgets import IntProgress\n",
    "\n",
    "import folium\n",
    "from folium.plugins import MarkerCluster\n",
    "\n",
    "import overpass\n",
    "\n",
    "import time\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "# import fiona as fn\n",
    "from shapely.geometry import shape, mapping, LineString\n",
    "from shapely.ops import cascaded_union, unary_union, polygonize\n",
    "import shapely\n",
    "import geojson\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import requests\n",
    "import json\n",
    "\n",
    "from zonal_stats import * # import zonal stats function from python file, get it here: https://gist.github.com/perrygeo/5667173\n",
    "\n",
    "api_ov = overpass.API(timeout=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is needed to donwload the OpenRouteService engine with:\n",
    "\n",
    "```bash\n",
    "wget https://github.com/GIScience/openrouteservice/archive/refs/tags/v7.1.0.zip\n",
    "unzip v7.1.0.zip\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, it is neeed to have installed [docker](https://www.docker.com/) and download the OpenStreetMap data of the area of study. This could be done using the [Geofabrick](https://download.geofabrik.de/) website.\n",
    "\n",
    "\n",
    "The file should be in the format similar to `italy-latest.osm.pbf`, the one used in this demo. Directly available at the [link](https://download.geofabrik.de/europe/italy.html), or at the [link](https://download.geofabrik.de/europe/italy-latest.osm.pbf). \n",
    "\n",
    "The file size is around 2 GB, the download speed depends on your connection. \n",
    "\n",
    "Copy the file inside the folder [openrouteservice-7.1.0\\docker](openrouteservice-7.1.0\\docker)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then it is need to modify the file [docker-compose.yml](openrouteservice-7.1.0/docker/docker-compose.yml) uncommenting the build option and modify the text to:\n",
    "\n",
    "```yml\n",
    "OSM_FILE: ./italy-latest.osm.pbf\n",
    "```\n",
    "```yml\n",
    "- ./italy-latest.osm.pbf:/home/ors/ors-core/data/osm_file.pbf\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then opening a terminal launch the following commands:\n",
    "\n",
    "```bash\n",
    "cd openrouteservice-7.1.0\n",
    "cd docker \n",
    "mkdir conf elevation_cache graphs logs\n",
    "cd logs\n",
    "mkdir ors tomcat \n",
    "cd ..\n",
    "docker compose up\n",
    "```\n",
    "\n",
    "At the first starts the streets graphs will be generated. The task required some minutes depending on the file and the machine used.\n",
    "\n",
    "The creation could be checked at the link [http://localhost:8080/ors/v2/status](http://localhost:8080/ors/v2/status)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name_region = \"Lombardia\"\n",
    "admin_level = \"4\" # \"2\" = State, \"4\" = Province\n",
    "charging_station_flag = 1 # 1 = charging stations, 0 = fuel stations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download of station from OSM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if charging_station_flag == 1:\n",
    "    query = \"[out:json][timeout:60];area[\\\"name\\\"=\\\"\" + name_region + \"\\\"][\\\"admin_level\\\"=\" + admin_level + \"]->.searchArea;node[\\\"amenity\\\"=\\\"charging_station\\\"](area.searchArea);out center;\"\n",
    "    text_stations = \"charging\"\n",
    "else:\n",
    "    query = \"[out:json][timeout:60];area[\\\"name\\\"=\\\"\" + name_region + \"\\\"][\\\"admin_level\\\"=\" + admin_level + \"]->.searchArea;node[\\\"amenity\\\"=\\\"fuel\\\"](area.searchArea);out center;\"\n",
    "    text_stations = \"petrol\"\n",
    "stations = api_ov.get(query, build = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "charging_dict = {}\n",
    "for feature in stations[\"elements\"]:\n",
    "    facility_id = int(feature['id'])\n",
    "    charging_dict[facility_id] = {\n",
    "        'geometry': {\"type\": \"Point\", \"coordinates\": [ feature[\"lon\"], feature[\"lat\"] ] }\n",
    "    }\n",
    "last = [feature[\"lat\"], feature[\"lon\"]]\n",
    "print(\"Created dictionary with \" + str(len(charging_dict)) + \" \" + text_stations + \" stations\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_outline = folium.Map(tiles='openstreetmap', location=(last), zoom_start=9)\n",
    "\n",
    "cluster = MarkerCluster().add_to(map_outline)  \n",
    "\n",
    "for facility_id in charging_dict:\n",
    "    folium.Marker(list(reversed(charging_dict[facility_id]['geometry']['coordinates']))).add_to(cluster)\n",
    "\n",
    "map_outline.save(os.path.join('results', \"1_\" + text_stations + \"_station_overview_\" + name_region + \".html\"))\n",
    "map_outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation of isochrones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "headers = {'Accept': 'application/json, application/geo+json, application/gpx+xml, img/png; charset=utf-8','Content-Type': 'application/json; charset=utf-8'}\n",
    "\n",
    "i = 0\n",
    "check = 0\n",
    "error_counter = 0\n",
    "\n",
    "iso_car_1 = []\n",
    "iso_car_5 = []\n",
    "iso_car_10 = []\n",
    "list_loc = []\n",
    "\n",
    "for facility_id in charging_dict.keys():\n",
    "    loc = charging_dict[facility_id]['geometry']['coordinates']\n",
    "    list_loc.append(loc)\n",
    "l = len(list_loc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = IntProgress(min=0, max=l) # instantiate the bar\n",
    "display(f) # display the bar\n",
    "\n",
    "while i < l:\n",
    "        iso_params = {'locations': [list_loc[i]],\n",
    "                'range_type': 'distance',\n",
    "                'range': [1000, 5000, 10000]  # meters\n",
    "                }\n",
    "                \n",
    "        call = requests.post('http://localhost:8080/ors/v2/isochrones/driving-car', json=iso_params, headers=headers)\n",
    "        request = json.loads(call.text)\n",
    "        check += 1\n",
    "        \n",
    "        try:\n",
    "                iso_car_1.append(shape(request['features'][0]['geometry']))\n",
    "                iso_car_5.append(shape(request['features'][1]['geometry']))\n",
    "                iso_car_10.append(shape(request['features'][2]['geometry']))\n",
    "                i += 1\n",
    "        except:\n",
    "                print(\"Error in \" + text_stations + \" stations at location: \" + str(list_loc[i]))\n",
    "                print(\"Check the error on the map at https://www.openstreetmap.org/#map=19/\" + str(list_loc[i][1]) + \"/\" + str(list_loc[i][0]))\n",
    "                error_counter += 1\n",
    "                i += 1\n",
    "\n",
    "        f.value += 1\n",
    "\n",
    "\n",
    "print(\"Have been requested \" + str(check) + \" \" + text_stations + \"stations\")\n",
    "if error_counter > 0:\n",
    "       print(\"Of which \" + str(error_counter) + \" ended with an error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iso_union_car_1 = cascaded_union(iso_car_1)\n",
    "iso_union_car_5 = cascaded_union(iso_car_5)\n",
    "iso_union_car_10 = cascaded_union(iso_car_10)\n",
    "geojson_iso_1 = shapely.to_geojson(iso_union_car_1)\n",
    "geojson_iso_5 = shapely.to_geojson(iso_union_car_5)\n",
    "geojson_iso_10 = shapely.to_geojson(iso_union_car_10)\n",
    "\n",
    "isochrones_car_filename_1 = \"results/iso_union_car_1_\" + text_stations + \"_\" + name_region + \".geojson\"\n",
    "isochrones_car_filename_5 = \"results/iso_union_car_5_\" + text_stations + \"_\" + name_region + \".geojson\"\n",
    "isochrones_car_filename_10 = \"results/iso_union_car_10_\" + text_stations + \"_\" + name_region + \".geojson\"\n",
    "\n",
    "# save isochrones to geojson\n",
    "with open(isochrones_car_filename_1, 'w') as f:\n",
    "    f.write(geojson_iso_1)\n",
    "with open(isochrones_car_filename_5, 'w') as f:\n",
    "    f.write(geojson_iso_5)\n",
    "with open(isochrones_car_filename_10, 'w') as f:\n",
    "    f.write(geojson_iso_10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_isochrones = folium.Map(tiles='openstreetmap', location=(last), zoom_start=9)\n",
    "\n",
    "def style_function(color):  # To style isochrones\n",
    "    return lambda feature: dict(color=color)\n",
    "\n",
    "folium.GeoJson(geojson_iso_10, name=\"10 km\", style_function = style_function('#fde725ff')).add_to(map_isochrones)\n",
    "folium.GeoJson(geojson_iso_5, name=\"5 km\",style_function = style_function('#1f968bff')).add_to(map_isochrones)\n",
    "folium.GeoJson(geojson_iso_1, name=\"1 km\",style_function = style_function('#440154ff')).add_to(map_isochrones)\n",
    "\n",
    "cluster = MarkerCluster().add_to(map_isochrones)  \n",
    "\n",
    "for facility_id in charging_dict:\n",
    "    folium.Marker(list(reversed(charging_dict[facility_id]['geometry']['coordinates']))).add_to(cluster)\n",
    "\n",
    "map_isochrones.save(os.path.join('results', \"2_\"  + text_stations + \"_isochrones_\" + name_region + \".html\"))\n",
    "map_isochrones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download of highways from OSM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_2 = \"[out:json][maxsize:1000000000][timeout:600];area[\\\"name\\\"=\\\"\" + name_region + \"\\\"][\\\"admin_level\\\"=\" + admin_level +\"]->.searchArea;(way[\\\"highway\\\"=\\\"motorway\\\"](area.searchArea);way[\\\"highway\\\"=\\\"trunk\\\"](area.searchArea);way[\\\"highway\\\"=\\\"primary\\\"](area.searchArea);way[\\\"highway\\\"=\\\"secondary\\\"](area.searchArea);way[\\\"highway\\\"=\\\"tertiary\\\"](area.searchArea);way[\\\"highway\\\"=\\\"motorway_link\\\"](area.searchArea);way[\\\"highway\\\"=\\\"trunk_link\\\"](area.searchArea);way[\\\"highway\\\"=\\\"primary_link\\\"](area.searchArea);way[\\\"highway\\\"=\\\"secondary_link\\\"](area.searchArea);way[\\\"highway\\\"=\\\"tertiary_link\\\"](area.searchArea);way[\\\"highway\\\"=\\\"unclassified\\\"](area.searchArea);way[\\\"highway\\\"=\\\"residential\\\"](area.searchArea);way[\\\"highway\\\"=\\\"living_street\\\"](area.searchArea);way[\\\"highway\\\"=\\\"service\\\"](area.searchArea););out geom;\"\n",
    "\n",
    "highways = api_ov.get(query_2, build = False)\n",
    "\n",
    "strade_dict = [{\n",
    "    'id': element['id'],\n",
    "    'class': element['tags']['highway'],\n",
    "    'geometry': LineString([(coords['lon'], coords['lat']) for coords in element['geometry']]),\n",
    "} for element in highways['elements']]\n",
    "\n",
    "highways_region = gpd.GeoDataFrame(strade_dict)\n",
    "highways_region = highways_region.set_crs('epsg:4326')\n",
    "highways_region.to_file(\"results/highways_\"+ name_region + \".geojson\", driver=\"GeoJSON\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intersection between isochrones and highways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iso_car_1_geo = gpd.read_file(\"results/iso_union_car_1_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")\n",
    "iso_car_5_geo = gpd.read_file(\"results/iso_union_car_5_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")\n",
    "iso_car_10_geo = gpd.read_file(\"results/iso_union_car_10_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")\n",
    "\n",
    "strade_iso_car_1 = highways_region[highways_region.intersects(iso_car_1_geo.geometry.iloc[0])]\n",
    "strade_iso_car_5 = highways_region[highways_region.intersects(iso_car_5_geo.geometry.iloc[0])]\n",
    "strade_iso_car_10 = highways_region[highways_region.intersects(iso_car_10_geo.geometry.iloc[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strade_iso_car_1.to_file(\"results/strade_iso_car_1_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")  \n",
    "strade_iso_car_5.to_file(\"results/strade_iso_car_5_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")  \n",
    "strade_iso_car_10.to_file(\"results/strade_iso_car_10_\" + text_stations + \"_\" + name_region + \".geojson\", driver=\"GeoJSON\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_streets = folium.Map(tiles='openstreetmap', location=(last), zoom_start=9)\n",
    "\n",
    "folium.GeoJson(highways_region.to_json(), name=\">10 km\", style_function = style_function('#F08080')).add_to(map_streets)\n",
    "folium.GeoJson(strade_iso_car_10.to_json(), name=\"10 km\", style_function = style_function('#fde725ff')).add_to(map_streets)\n",
    "folium.GeoJson(strade_iso_car_5.to_json(), name=\"5 km\",style_function = style_function('#1f968bff')).add_to(map_streets)\n",
    "folium.GeoJson(strade_iso_car_1.to_json(), name=\"1 km\",style_function = style_function('#440154ff')).add_to(map_streets)\n",
    "\n",
    "map_streets.save(os.path.join('results', \"3_streets_isocrones\" + text_stations + \"_\" + name_region + \".html\"))\n",
    "# The date in the map are significant and the load of the map could slow al lot the results\n",
    "# it is suggested to open the html in a broser\n",
    "# map_streets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l_1 = strade_iso_car_1.to_crs('epsg:3857').geometry.length.sum()\n",
    "l_5 = strade_iso_car_5.to_crs('epsg:3857').geometry.length.sum() - l_1\n",
    "l_10 = strade_iso_car_10.to_crs('epsg:3857').geometry.length.sum() - l_5 - l_1\n",
    "l_region = highways_region.to_crs('epsg:3857').geometry.length.sum()\n",
    "l_over_10 = l_region - l_10 - l_5 - l_1\n",
    "\n",
    "data = {\n",
    "    \"Distanza da colonnina\": [\"< 1 km\", \"1 - 5 km\", \"5 - 10 km\", \"> 10 km\"],\n",
    "    \"Lunghezza totale (km)\": [l_1/1000, l_5/1000, l_10/1000, l_over_10/1000],\n",
    "    \"Percentuale classe (%)\": [l_1/l_region*100, l_5/l_region*100, l_10/l_region*100, l_over_10/l_region*100]\n",
    "}\n",
    "\n",
    "pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_distribution = pd.DataFrame(data)\n",
    "final_distribution.to_csv(\"results/distribution_\" + text_stations + \"_\" + name_region + \".csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gis2",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}