notebooks/311-folium-with-density.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7e69053b-393d-49e6-91cb-2708bdaf6e1f",
   "metadata": {},
   "source": [
    "# Spatial Analysis and Map Visualizations\n",
    "\n",
    "This analysis is for the graffiti removal request type.  Five visualizations are used to tell the story:\n",
    "  1.  Dataframe Statistics\n",
    "  2.  Choropleth map showing annual 311 counts for LA\n",
    "  3.  Density choropleth\n",
    "  4.  Static heat map - high density neighborhood\n",
    "  5.  Animated heat map - high density neighborhood\n",
    "  6.  Annimated heat map showing daily 311 calls for LA\n",
    "  \n",
    "This nb is intended to show the ideas.  For starters I'm just using graffiti reports for 2021.  NC's are based on the certified NC dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43418e3e-4d3a-4765-b242-1063f3266881",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%run start.py\n",
    "\n",
    "import folium\n",
    "from folium import plugins\n",
    "from utils import read_new311_shape, dt_to_object, read_ncs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79b15a7f-6a29-4732-a196-2ed3837bb3d9",
   "metadata": {},
   "source": [
    "Two data sets provide the foundation data:\n",
    "\n",
    "  1. The graffiti subset of the cleaned 311 data\n",
    "  2. The cleaned version of the certified neighborhood councils\n",
    "\n",
    "I'm using the functions developed in other notebooks to read the data.  I've also added some date abstractions for upstream processing/querying."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "547b1fc6-e999-4768-8d95-833471889108",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#graffiti_gdf = read_new311_shape('../data/311/graffiti.geojson')\n",
    "graffiti_gdf = gpd.read_parquet('../data/311/graffiti.parq')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30bf3c4a-716f-4f3b-be62-e6ba429ff188",
   "metadata": {},
   "outputs": [],
   "source": [
    "graffiti_gdf['dayofyear'] = graffiti_gdf['created_dt'].apply(lambda dt: dt.dayofyear)\n",
    "graffiti_gdf['day_name'] = graffiti_gdf['created_dt'].apply(lambda dt: dt.day_name())\n",
    "graffiti_gdf['day'] = graffiti_gdf['created_dt'].apply(lambda dt: dt.date().strftime(\"%m-%d (%A)\"))\n",
    "graffiti_gdf['month'] = graffiti_gdf['created_dt'].apply(lambda dt: dt.month)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a22b620-334b-4b5c-a913-b1436e190a6e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "neighborhoods_gdf = read_ncs()\n",
    "\n",
    "graffiti_counts = graffiti_gdf['nc'].value_counts().to_frame().reset_index().rename(columns={'index': 'nc_id', 'nc': 'count'})\n",
    "\n",
    "graffiti_merged_gdf = pd.merge(neighborhoods_gdf, graffiti_counts, how=\"left\", on=[\"nc_id\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcc7c490-4d5c-4c65-a48d-cc893b2b7528",
   "metadata": {},
   "source": [
    "At this point we have the two baseline data frames that will be used in the following steps."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a828d39-2297-4feb-8f02-f46817404537",
   "metadata": {},
   "source": [
    "# 1 - Bar Chart\n",
    "\n",
    "We'll start with a very basic understanding of the 311 dataframe.  This simple graphic shows annual count / NC.  For simplicity I'm only showing the top 25."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b2832a8-df3d-4e47-a6d0-6076b199a7e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "graffiti_gdf['nc_name'].value_counts()[:25].plot.barh(figsize=(8, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a33e3a9-66cb-4875-92d4-b41f3081b662",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 2 - Choropleth for LA\n",
    " \n",
    "The Oxford Dictionary defines a choropleth map as \n",
    "\n",
    "`a map that uses differences in shading, coloring, or the placing of symbols within predefined areas to indicate the average values of a property or quantity in those areas.`\n",
    "\n",
    "Specifics of the model:\n",
    "\n",
    "  1. The coloring scheme ranges from Yellow -> Orange -> Red\n",
    "  2. Annual reported 311 reports (based on created date) are grouped by Neighborhood Council\n",
    "  3. These counts are then binned into 6 groups and mapped on the color line\n",
    "  4. Interpret as higher counts are \"brighter\" red, lower counts are yellow\n",
    "  5. This visualiztion supports high level analysis (City Council, empowerla, Service Regions, ...)\n",
    "\n",
    "The receipe for the choropleth - https://vverde.github.io/blob/interactivechoropleth.html\n",
    "\n",
    "Note:\n",
    "Data is a combination of counts from the 311 observations and the NC polygons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96e28fd3-4096-4c4e-baf6-6201f95d5295",
   "metadata": {},
   "outputs": [],
   "source": [
    "choropleth_map = folium.Map(location = [34.05, -118.25], zoom_start = 10)\n",
    "plugins.Fullscreen().add_to(choropleth_map)\n",
    "folium.Choropleth(\n",
    " geo_data=graffiti_merged_gdf,\n",
    " name='Choropleth',\n",
    " data=graffiti_merged_gdf,\n",
    " columns=['nc_id','count'],\n",
    " key_on=\"feature.properties.nc_id\",\n",
    " fill_color='YlOrRd',\n",
    " #threshold_scale=myscale,\n",
    " fill_opacity=1,\n",
    " line_opacity=0.2,\n",
    " legend_name='311 request counts',\n",
    " smooth_factor=0\n",
    ").add_to(choropleth_map)\n",
    "\n",
    "style_function = lambda x: {'fillColor': '#ffffff', \n",
    "                            'color':'#000000', \n",
    "                            'fillOpacity': 0.1, \n",
    "                            'weight': 0.1}\n",
    "highlight_function = lambda x: {'fillColor': '#000000', \n",
    "                                'color':'#000000', \n",
    "                                'fillOpacity': 0.50, \n",
    "                                'weight': 0.1}\n",
    "nc_info = folium.features.GeoJson(\n",
    "    graffiti_merged_gdf,\n",
    "    style_function=style_function, \n",
    "    control=False,\n",
    "    highlight_function=highlight_function, \n",
    "    tooltip=folium.features.GeoJsonTooltip(\n",
    "        fields=['name', 'nc_id', 'count'],\n",
    "        aliases=['Neighborhood: ','NC ID: ','311 count: '],\n",
    "        style=(\"background-color: white; color: #333333; font-family: arial; font-size: 12px; padding: 10px;\") \n",
    "    )\n",
    ")\n",
    "choropleth_map.add_child(nc_info)\n",
    "choropleth_map.keep_in_front(nc_info)\n",
    "folium.LayerControl().add_to(choropleth_map)\n",
    "\n",
    "\n",
    "choropleth_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21187566-8f52-444b-8969-bc1f29ccdbc9",
   "metadata": {},
   "source": [
    "Explore the Neighborhood Councils on the choropleth.  Hover over the darkest read NC.  That is the one we'll explore in a bit more detail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e183f7dc-5fe0-4e7d-b18d-1fd6fbc9a557",
   "metadata": {},
   "source": [
    "# 3 - Density Choropleth\n",
    "\n",
    "So I'm curious about normalizing the Graffiti Removal counts for the choropleth.\n",
    "\n",
    "Some thoughts:\n",
    "\n",
    "  1. The OCB might be interested in distribution of raw counts?\n",
    "  2. Residents may be interested in count per area?\n",
    "  \n",
    "Let's run with these thoughts and add a choropleth showing density (count/sq KM) in each NC.  Probably beating a dead horse but what the hey!\n",
    "\n",
    "**Note:** This section is a bit messy.  I'm using code snippets from the population density notebook and reusing the choropleth code from above.  This should be refactored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "211c0d03-c8bc-4929-8ed1-fd35008619f3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pyproj import Geod\n",
    "\n",
    "geod = Geod(ellps=\"WGS84\")\n",
    "\n",
    "def square_miles(geo):\n",
    "    square_meters = abs(geod.geometry_area_perimeter(geo)[0])\n",
    "    return (square_meters * 10.764) / 27878000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0692453f-9505-4200-acf8-3c6af8edeaf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "graffiti_merged_gdf['sq_miles'] = graffiti_merged_gdf.apply(lambda row: square_miles(row.geometry), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf0b11b7-b33b-4103-994e-27a8c2ac05ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "graffiti_merged_gdf['density'] = graffiti_merged_gdf.apply(lambda row: row['count'] / row['sq_miles'], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0fca443-476c-4774-a310-6eedcccb77fc",
   "metadata": {},
   "source": [
    "I just added a density column to the gdf.  We can use that for a new choropleth.  The results are interesting if not predictable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "026165c3-c60e-465d-a544-c7a395bfeefb",
   "metadata": {},
   "outputs": [],
   "source": [
    "density_choropleth_map = folium.Map(location = [34.05, -118.25], zoom_start = 10)\n",
    "plugins.Fullscreen().add_to(density_choropleth_map)\n",
    "folium.Choropleth(\n",
    " geo_data=graffiti_merged_gdf,\n",
    " name='Choropleth Density',\n",
    " data=graffiti_merged_gdf,\n",
    " columns=['nc_id','density'],\n",
    " key_on=\"feature.properties.nc_id\",\n",
    " fill_color='YlOrRd',\n",
    " #threshold_scale=myscale,\n",
    " fill_opacity=1,\n",
    " line_opacity=0.2,\n",
    " legend_name='311 request density',\n",
    " smooth_factor=0\n",
    ").add_to(density_choropleth_map)\n",
    "\n",
    "style_function = lambda x: {'fillColor': '#ffffff', \n",
    "                            'color':'#000000', \n",
    "                            'fillOpacity': 0.1, \n",
    "                            'weight': 0.1}\n",
    "highlight_function = lambda x: {'fillColor': '#000000', \n",
    "                                'color':'#000000', \n",
    "                                'fillOpacity': 0.50, \n",
    "                                'weight': 0.1}\n",
    "density_nc_info = folium.features.GeoJson(\n",
    "    graffiti_merged_gdf,\n",
    "    style_function=style_function, \n",
    "    control=False,\n",
    "    highlight_function=highlight_function, \n",
    "    tooltip=folium.features.GeoJsonTooltip(\n",
    "        fields=['name', 'nc_id', 'density'],\n",
    "        aliases=['Neighborhood: ','NC ID: ','311 density: '],\n",
    "        style=(\"background-color: white; color: #333333; font-family: arial; font-size: 12px; padding: 10px;\") \n",
    "    )\n",
    ")\n",
    "density_choropleth_map.add_child(density_nc_info)\n",
    "density_choropleth_map.keep_in_front(density_nc_info)\n",
    "folium.LayerControl().add_to(density_choropleth_map)\n",
    "\n",
    "\n",
    "density_choropleth_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b321ed3-34da-4f78-82a4-45d54e8db414",
   "metadata": {},
   "source": [
    "What I want to do is compare the two choropleths.  I could scroll up and down, break out the two output cells or ...\n",
    "\n",
    "Why don't we use two output widgets to look at them side by side!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24ab52c7-54c3-4ec3-8324-b29c9b4dd033",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "count_output = Output(layout={'border': '1px solid black',\n",
    "                            'width': '50%'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10a815b-4684-4af0-b404-14ca77cc5bf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "density_output = Output(layout={'border': '1px solid black',\n",
    "                            'width': '50%'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a78d573-fb2c-427e-a591-cb8bf1bff90b",
   "metadata": {},
   "outputs": [],
   "source": [
    "with count_output:\n",
    "    display(choropleth_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70e353e8-54f3-4449-a939-a1a398e23da6",
   "metadata": {},
   "outputs": [],
   "source": [
    "with density_output:\n",
    "    display(density_choropleth_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a9c9946-b3dd-4d04-97af-253c660b3091",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "HBox([count_output, density_output])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f81e858-7dda-4900-9527-9dc3505098aa",
   "metadata": {},
   "source": [
    "That's enough.  You tell me if this helps the story."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38b188d5-8408-449e-ac20-f308506af420",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 4 - Static Heat Map for South Central\n",
    "\n",
    "The definition I'll use for heatmaps comes from wikipedia - `A heat map is a data visualization technique that shows magnitude of a phenomenon as color.`  More specifically, in our case we're using the concept of geospatial heatmaps, `measure counts by color from green -> red so brighter red represents a greater density of 311 reports.`  First we'll look in 2-D.\n",
    "\n",
    "Processing steps for this visualization:\n",
    "\n",
    "1. To simplify I will just use the last quarter (2021) data\n",
    "2. Use the nc id (78) from the visual inspection above to select 311 observations\n",
    "3. Use the neighborhood polygon to center the map\n",
    "4. Build and display the heat map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "007ac03c-f6e1-4988-a6b5-bec1050a3aa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "south_central_gdf = graffiti_gdf.query(f\"nc == 78\", engine='python').reset_index().drop(columns=['index'])\n",
    "\n",
    "south_central_quarter_gdf = south_central_gdf.query(f\"month > 9\", engine='python').reset_index().drop(columns=['index'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8999b74-7936-4ac0-8897-1933394ea2d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "south_central_poly_geometry = neighborhoods_gdf.query(f\"nc_id == 78\", engine='python').reset_index().iloc[0]['geometry']\n",
    "\n",
    "map_center = [south_central_poly_geometry.centroid.y, south_central_poly_geometry.centroid.x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f8b8e5-875a-41d0-bcce-eaf6e2e86cb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "static_hmap = folium.Map(location = map_center, tiles='Stamen Toner', zoom_start = 15)\n",
    "plugins.Fullscreen().add_to(static_hmap)\n",
    "heat_data = [[point.xy[1][0], point.xy[0][0]] for point in south_central_quarter_gdf.geometry ]\n",
    "plugins.HeatMap(heat_data).add_to(static_hmap)\n",
    "static_hmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27d24b4f-795c-4770-8824-fe0abcc2fd33",
   "metadata": {},
   "source": [
    "Looking at these intensities, I see a pattern.  Do you?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96d31c04-7403-4156-8567-7d5cbfec5d17",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 5 - Animated Heat Map for South Central\n",
    "\n",
    "For this viusalization I want to \"see\" the 311 calls in the static heatmap over time.  Folium as a plugin to do just that.  I found a useful receipe here - https://stackoverflow.com/questions/64325958/heatmapwithtime-plugin-in-folium. \n",
    "\n",
    "\n",
    "\n",
    "  1. The data for the visualization uses the same data set for the static heatmap\n",
    "  2. The time component is when the 311 report was created\n",
    "  3. I do a simple mapping to day-of-the-year for the created_dt\n",
    "  4. The animation steps through the observations by day-of-the-year\n",
    "  5. The animation control is placed bottom left in the map display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9456b086-0d0b-4fed-932c-08f1480f6b86",
   "metadata": {},
   "outputs": [],
   "source": [
    "south_central_hmap = folium.Map(location = map_center, tiles='Stamen Toner', zoom_start = 14)\n",
    "plugins.Fullscreen().add_to(south_central_hmap);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f53ac2e6-b2a2-4666-be02-4474524d5c34",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict, OrderedDict\n",
    "\n",
    "south_central_data = defaultdict(list)\n",
    "for r in south_central_quarter_gdf.itertuples():\n",
    "    south_central_data[r.day].append([r.latitude, r.longitude])\n",
    "    \n",
    "south_central_data = OrderedDict(sorted(south_central_data.items(), key=lambda t: t[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c2f6484-84d8-4bf7-8209-3e0c969b560f",
   "metadata": {},
   "outputs": [],
   "source": [
    "south_central_hm = plugins.HeatMapWithTime(data=list(south_central_data.values()),\n",
    "                     index=list(south_central_data.keys()), \n",
    "                     radius=10,\n",
    "                     auto_play=True,\n",
    "                     max_opacity=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc64f48-a6d9-44dc-aaee-5fcb71a62487",
   "metadata": {},
   "outputs": [],
   "source": [
    "south_central_hm.add_to(south_central_hmap)\n",
    "south_central_hmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7e265b7-b7b1-4d8f-8087-f0985850015b",
   "metadata": {},
   "source": [
    "Look at the animation in full screen mode.  You can see that most reports seem to be on weekdays?  That may warrent further investigation at some point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b2dc278-76b7-4b98-ae9d-3c41171bc72a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 6 - Animated Heat Map for LA\n",
    "\n",
    "Now I return to the whole 311 data set for graffiti removal.  We started with the choropleth for all of LA, did an excursion to the high density neighborhood, and now we'll go back to the entire city for the entire year.\n",
    "\n",
    "  1. The data for the visualization uses 311 data for the graffiti\n",
    "  2. The time component is when the 311 report was created\n",
    "  3. I do a simple mapping to day-of-the-year for the created_dt\n",
    "  4. The animation steps through the observations by day-of-the-year\n",
    "  5. The animation control is placed bottom left in the map display\n",
    "\n",
    "\n",
    "\n",
    "This receipe helped.  https://stackoverflow.com/questions/64325958/heatmapwithtime-plugin-in-folium\n",
    "\n",
    "Note:  This is covering the whole city so it's a bit hard to get a good starting zoom level.  Once the animation starts you should go to full screen mode and the navigate around zoom in (with the mouse scoll).  That's a good way to see what's happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04e62153-efb2-45c2-8ba0-caae4f42c31d",
   "metadata": {},
   "outputs": [],
   "source": [
    "hmap = folium.Map(location = [34.05, -118.25], tiles='Stamen Toner', zoom_start = 10)\n",
    "plugins.Fullscreen().add_to(hmap);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4403f255-1b21-479b-904f-354c19bdb2e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict, OrderedDict\n",
    "\n",
    "data = defaultdict(list)\n",
    "for r in graffiti_gdf.itertuples():\n",
    "    data[r.day].append([r.latitude, r.longitude])\n",
    "    \n",
    "data = OrderedDict(sorted(data.items(), key=lambda t: t[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73e91fae-1297-44e4-b13d-fb974ceae7c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "hm = plugins.HeatMapWithTime(data=list(data.values()),\n",
    "                     index=list(data.keys()), \n",
    "                     radius=10,\n",
    "                     auto_play=True,\n",
    "                     max_opacity=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a50ebb8e-1df5-4660-82a6-814324a2f468",
   "metadata": {},
   "outputs": [],
   "source": [
    "hm.add_to(hmap)\n",
    "hmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "222a1f4e-67fa-48ee-b15f-83141dcd706b",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "I've created a basic story flow going from a large scale choropleth map to small scale, area specific maps, and finally back to the large scale animation.  This should give you some sense for what happened in 2021 for the graffiti removal request type.  It also highlights some other areas for exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bdc1685-3a57-4179-be5c-14cf51a5764d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}