notebooks/public_data_osm_deck.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "# **Toilet & Bench Coverage***\n",
    "###### *inside an area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "# Imports and credentials\n",
    "For this recipe :) we will need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd, numpy as np\n",
    "#geo libraries\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point, Polygon\n",
    "from shapely import speedups #optional speed-ups\n",
    "speedups.enable() \n",
    "import h3\n",
    "import overpass\n",
    "#graphics\n",
    "import pydeck as pdk\n",
    "#a mapbox api key for tiles (get your own at mapbox.com)\n",
    "with open('./.to_ignore/keys/v1_mapbox.key', 'r') as f:\n",
    "    mapbox_key= f.read()[:-1]\n",
    "#and a choice of mapbox tiles\n",
    "map_style='mapbox://styles/mapbox/light-v10'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "## OpenStreetmap's Data + Overpass API\n",
    "&nbsp;\n",
    "\n",
    "#### Below is a simple method that build geoDataFrames from:\n",
    "    * Querying OSM for all nodes/ways\n",
    "    * Tagged as \"amenity\"with \"bench\" or \"toilets\" as a value\n",
    "    * In a geographical area given by name and administrative level (see [here](https://www.openstreetmap.org/relation/52822) )\n",
    "\n",
    "&nbsp;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def fetch_amenity_gdf(amenity_types, area_name, adm_lvl, timeout_sec=90):\n",
    "    #exports a shapefile with OSM's response to a request of data for\n",
    "    #amenities of amenity_type in area_name / adm_lvl.\n",
    "    query=f'( area[admin_level={adm_lvl}][\"name\"=\"{area_name}\"]; ) ->.a; \\n('\n",
    "    node_req='\\nnode[\"amenity\"=\"{0}\"](area.a);\\nway[\"amenity\"=\"{0}\"](area.a);'\n",
    "    for x in amenity_types:\n",
    "        query=query+node_req.format(x)\n",
    "    query=query+'\\n);'\n",
    "    print(query)\n",
    "    api = overpass.API(debug=False, timeout=timeout_sec)\n",
    "    try:\n",
    "        result = api.Get(query, responseformat=\"geojson\")\n",
    "        result = gpd.GeoDataFrame.from_features(result['features'])\n",
    "        result.to_file('./.to_ignore/out/'\n",
    "                       +f'gdf_amenities_{area_name}_{adm_lvl}') \n",
    "        print('Downloaded to ./.to_ignore/out/'\n",
    "              +f'gdf_amenities_{area_name}_{adm_lvl}')\n",
    "    except Exception as HELL:\n",
    "        print(HELL)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "\n",
    "# Fetch & clean geoDataFrame from OSM\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "#### With features on toilets and benches\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Avoid querying if file already exists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "area='Sverige'#'Göteborgs Stad'\n",
    "lvl='2'#'7' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    gdf_amenity=gpd.read_file('./.to_ignore/out/'+f'gdf_amenities_{area}_{lvl}')\n",
    "except:\n",
    "    print('file_not_found, fetching amenities from OSM')\n",
    "    fetch_amenity_gdf([\"toilets\",\"bench\"], area,lvl)\n",
    "    gdf_amenity=gpd.read_file('./.to_ignore/out/'+f'gdf_amenities_{area}_{lvl}')\n",
    "gdf_amenity.dropna(subset=['geometry'], inplace=True) #remove nans on geometry\n",
    "gdf_amenity = gdf_amenity.set_crs(epsg=4326) # this is the CRS in OSM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "\n",
    "## What features do we have here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Features: amenity,created_by,wheelchair,backrest,source,material,fee,layer,access,table,tourism,seamark_sm,seamark_ty,fixme,name,descriptio,note,operator,emergency,indoor,contact_ph,opening_ho,drinking_w,viewpoint,colour,fee_price,toilets_di,unisex,descript_1,women,check_date,direction,seats,amount,changing_t,capacity,capacity_d,level,snowmobile,female,sanitary_d,public_tra,shelter,toilets_wh,bench_coun,leisure,toilets_po,supervised,network,toilets_ha,ele,route_ref,ref,building,survey_dat,artwork,diaper,male,toilets_ac,designatio,natural,note_sv,highway,fireplace,bench,toilets__1,toilets_pa,water,descript_2,addr_city,addr_count,addr_house,addr_stree,wheelcha_1,charge,history,covered,shop,entrance,heated,place,addr_postc,website,historic,type,toilets_ch,Bing,surface,water_sour,public,artist_nam,start_date,payment_ca,payment_cr,seasonal,changing_1,levels,source_geo,bin,ele_msl,ele_source,source_pos,tactile_pa,phone,service_bi,service__1,service__2,service__3,service__4,toilets,internet_a,lit,abandoned_,openfire,changing_2,mapillary,capacity_w,source_nam,shower,roof,height,wikimedia_,lamp_type,image,handwashin,handwash_1,toilets_ho,payment_sw,payment_co,payment_de,location,board_type,informatio,undefined,geocaching,name_en,landuse,advertisin,roof_shape,barrier,source_add,toilets_nu,marked,permanent,building_p,building_l,area,outdoor_se,addr_place,amenity_1,building_1,est_height,name_de,camp_site,building_c,roof_colou,roof_heigh,capacity_m,wikidata,local_ref,building_m,roof_mater,backcountr,shelter_ty,toilets__2,room,geometry\n"
     ]
    }
   ],
   "source": [
    "print('Features: '+','.join(gdf_amenity.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Woa. Cambrian explosion of Features.\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "#### Let's only keep a reduced set for today's analysis :). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>amenity</th>\n",
       "      <th>water</th>\n",
       "      <th>drinking_w</th>\n",
       "      <th>fee</th>\n",
       "      <th>wheelchair</th>\n",
       "      <th>wheelcha_1</th>\n",
       "      <th>backrest</th>\n",
       "      <th>name</th>\n",
       "      <th>geocaching</th>\n",
       "      <th>note</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POINT (14.21335 56.59857)</td>\n",
       "      <td>toilets</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>POINT (18.05524 59.33863)</td>\n",
       "      <td>toilets</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    geometry  amenity water drinking_w   fee wheelchair  \\\n",
       "0  POINT (14.21335 56.59857)  toilets  None       None  None       None   \n",
       "1  POINT (18.05524 59.33863)  toilets  None       None  None       None   \n",
       "\n",
       "  wheelcha_1 backrest  name geocaching  note  \n",
       "0       None     None  None       None  None  \n",
       "1       None     None  None       None  None  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf_amenity=gdf_amenity[['geometry', 'amenity', 'water' ,'drinking_w','fee',\n",
    "        'wheelchair', 'wheelcha_1', 'backrest', 'name', 'geocaching','note']]\n",
    "gdf_amenity[:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "# Now lets have a look at this features in a map\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "#### We will use PyDeck for Jupyter (Deck's Jupyter plugin)\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "**To enable pydeck for Jupyter, run this in command before starting the notebook:**\n",
    "\n",
    "jupyter nbextension install --sys-prefix --symlink --overwrite --py pydeck\n",
    "\n",
    "jupyter nbextension enable --sys-prefix --py pydeck\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "##### For more info and examples, [check here](https://github.com/visgl/deck.gl/blob/master/bindings/pydeck/examples/geopandas_integration.py)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "#create separate lng and lat facilitates plotting different things\n",
    "gdf_amenity['lng']= list(map(lambda p: p.x, gdf_amenity.geometry ))\n",
    "gdf_amenity['lat']= list(map(lambda p: p.y, gdf_amenity.geometry ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Now, let's prepare the layers\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Say we want to start with simple points in different colours (see docs [here](https://deck.gl/docs/api-reference/layers/scatterplot-layer))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "opacity=0.5\n",
    "radius=40\n",
    "kwargs={ 'pickable':True,  'get_position':['lng', 'lat'],  'opacity':opacity}\n",
    "toilet_point_ly = pdk.Layer(\"ScatterplotLayer\", getFillColor=[237,102,99],#ed6663\n",
    "                        data=gdf_amenity.loc[gdf_amenity.amenity=='toilets'],\n",
    "                        getRadius=radius*1.5, **kwargs) \n",
    "bench_point_ly = pdk.Layer(\"ScatterplotLayer\", getFillColor=[78,137,174],#4e89ae\n",
    "                        data=gdf_amenity.loc[gdf_amenity.amenity=='bench'],\n",
    "                        getRadius=radius, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### And now, we also want some way to aggregate points\n",
    "\n",
    "&nbsp;\n",
    "#### We can try using Hexagon heatmaps from H3\n",
    "There are a couple of ways to do that, today we will use the simplest Hexagon layer.\n",
    "\n",
    "&nbsp;\n",
    "##### ...Then, we need to do one more thing: \n",
    "&nbsp;\n",
    "That is assigning H3 ids to each point.\n",
    "\n",
    "&nbsp;\n",
    "##### Trick question: Which resolutions do we want?\n",
    "For observational purposes I chose the following cases:\n",
    "\n",
    "| Resolution | Hex Edge Km  Avg|\n",
    "| --- | --- |\n",
    "| 8 | 0.46 |\n",
    "| 6 | 3.23 |\n",
    "| 4 | 22.61 |\n",
    "\n",
    "For more on resolutions, see [this link](https://h3geo.org/docs/core-library/restable)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### To colour the hexes, we might want to use the count of points belonging to that hexagon. \n",
    "\n",
    "&nbsp;\n",
    "\n",
    "#### This is how we construct such dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "code_folding": [],
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for resolution in (4,6,8):\n",
    "    gdf_amenity[f'h3id_{resolution}']=[h3.geo_to_h3(gdf_amenity.lat[x],\n",
    "                    gdf_amenity.lng[x],resolution) for x in gdf_amenity.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>h3id_4</th>\n",
       "      <th>cnt</th>\n",
       "      <th>norm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8408867ffffffff</td>\n",
       "      <td>246</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>841f05bffffffff</td>\n",
       "      <td>77</td>\n",
       "      <td>0.310204</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>841f251ffffffff</td>\n",
       "      <td>70</td>\n",
       "      <td>0.281633</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            h3id_4  cnt      norm\n",
       "0  8408867ffffffff  246  1.000000\n",
       "1  841f05bffffffff   77  0.310204\n",
       "2  841f251ffffffff   70  0.281633"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cnt_dfs_bench={} #and we store a dataframe per resolution\n",
    "cnt_dfs_toilet={}\n",
    "nrm_func=lambda x: (x-mn)/dffx\n",
    "for res in (4,6,8):\n",
    "    cnt_dfs_bench[res]=gdf_amenity.loc[gdf_amenity.amenity=='bench']\\\n",
    "            [f'h3id_{res}']\\\n",
    "            .value_counts().to_frame().rename(columns={f'h3id_{res}':f'cnt'})\\\n",
    "            .reset_index().rename(columns={'index':f'h3id_{res}'})\n",
    "    cnt_dfs_toilet[res]=gdf_amenity.loc[gdf_amenity.amenity=='toilets']\\\n",
    "            [f'h3id_{res}']\\\n",
    "            .value_counts().to_frame().rename(columns={f'h3id_{res}':f'cnt'})\\\n",
    "            .reset_index().rename(columns={'index':f'h3id_{res}'})\n",
    "    #feature scale normalization\n",
    "    mn=cnt_dfs_bench[res].cnt.min()\n",
    "    dffx=cnt_dfs_bench[res].cnt.max()-mn\n",
    "    cnt_dfs_bench[res]['norm']=cnt_dfs_bench[res].cnt.apply(nrm_func)\n",
    "    mn=cnt_dfs_toilet[res].cnt.min()\n",
    "    dffx=cnt_dfs_toilet[res].cnt.max()-mn\n",
    "    cnt_dfs_toilet[res]['norm']=cnt_dfs_toilet[res].cnt.apply(nrm_func)\n",
    "cnt_dfs_toilet[4][:3] #for resolution 4, show first 3 lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### (Optional) Let's generate a viewport pointing somewhere that makes sense in the map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "view_state = pdk.ViewState(latitude=57.71, longitude=11.94, zoom=6, bearing=0, \n",
    "                           pitch=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Now we create our Hex layers :)\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "For the sake of simplicity, let's look at the biggest resolution first (4)\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Check out the docs for Deck's  Hexagon Layer [here](https://pydeck.gl/gallery/h3_hexagon_layer.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "resolution=6 #4,6,8\n",
    "kwargs={'get_line_color':'[50,130,184]', 'line_width_min_pixels':0.3, \n",
    "        'pickable':True, 'filled':True, 'stroked':True, 'extruded':False, \n",
    "        'opacity':0.28, 'get_hexagon':f'h3id_{resolution}' }\n",
    "toilet_hex_ly = pdk.Layer(\"H3HexagonLayer\",  cnt_dfs_toilet[resolution], \n",
    "            getFillColor=\"[78+10/(norm*count),160-100*(norm+0.2),4000*norm]\",\n",
    "            **kwargs) \n",
    "bench_hex_ly = pdk.Layer(\"H3HexagonLayer\", cnt_dfs_bench[resolution], \n",
    "            getFillColor=\"[78+10/(norm*count),160-100*(norm+0.2),4000*norm]\",\n",
    "            **kwargs) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## And now let's draw our map!\n",
    "\n",
    "Note: Be careful, playing with the fill colors **can** be a bit addictive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ze_map=pdk.Deck(layers=[bench_point_ly,bench_hex_ly], map_style=map_style,\n",
    "initial_view_state=view_state, mapbox_key=mapbox_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ze_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " #### Exercises: \n",
    " * How does the world look like when tracking toilets insted of benches?\n",
    " * Can you find the optimal tiles given max number of toilets and benches?\n",
    "     * We will treat this one in another blog post, but have fun on your own first if you like :)!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "&nbsp;\n",
    "\n",
    "# Thanks for reading!\n",
    "\n",
    "&nbsp;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "&nbsp;\n",
    "\n",
    "##### If you are curious about _what else_ is well documented in OpenStreetMap, check out [overpass-turbo](http://overpass-turbo.eu/) and the [overpass python wrapper](https://github.com/mvexel/overpass-api-python-wrapper/tree/master/examples).\n",
    "\n",
    "&nbsp;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}