import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import rasterio as rio
import json
import raster_tools as rtpd.set_option('display.max_columns', 500)
pd.set_option('display.max_rows', 500)# for each geojson in data/missoula/missoula_poi, load the geojson into a geodataframe with the name of the file as the name of the geodataframe
# get filenames with glob
import glob
filenames = glob.glob('data/missoula/missoula_poi/*.geojson')
# create a dictionary to hold the geodataframes
gdfs = {}
# loop through the filenames
for filename in filenames:
# get the name of the file
name = filename.split('/')[-1].split('.')[0]
# load the geojson into a geodataframe
gdfs[name] = gpd.read_file(filename)
# add a column to the geodataframe with the name of the file
gdfs[name]['b_type'] = name# combine the poi geodataframes into a single geodataframe, with a new column 'type' that is the name of the original geodataframe
poi = gpd.GeoDataFrame(pd.concat(gdfs.values(), ignore_index=True))
poi['b_type'] = poi['b_type'].astype(str)
# poi['b_type'] = poi['b_type'].str.split('_').str[1:2].str[0
# split and keep everything after the first _
poi['b_type'] = poi['b_type'].str.split('_').str[1:].str.join('_')
poi['lat'] = poi.geometry.y
poi['lon'] = poi.geometry.x
poi.b_type# distribution of business types
poi.b_type.value_counts()# write pois to geojosn
poi.to_file('data/missoula/missoula_poi_all.geojson', driver='GeoJSON')poi.to_crs(epsg=5070, inplace=True)# load missoula_landuse_fixed.geojson
landuse = gpd.read_file('data/missoula/missoula_landuse_fixed.geojson')import pandas as pd
import folium
from folium.plugins import MarkerCluster
# Load the POI data
# poi = pd.read_csv('path_to_your_poi_data.csv') # Replace with actual path to your POI data
# Create a map centered on Missoula
m_pois = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a color map for business types
unique_types = poi['b_type'].unique()
color_map = {b_type: f'#{hash(b_type) % 0xFFFFFF:06x}' for b_type in unique_types}
# Create a feature group for all POIs
all_pois = folium.FeatureGroup(name="All POIs", show=True)
# Add all POIs to the all_pois group
for _, row in poi.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=color_map[row['b_type']],
fill=True,
fillColor=color_map[row['b_type']]
).add_to(all_pois)
all_pois.add_to(m_pois)
# Create a feature group for each business type
for b_type in unique_types:
fg = folium.FeatureGroup(name=b_type, show=False)
type_df = poi[poi['b_type'] == b_type]
for _, row in type_df.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=color_map[b_type],
fill=True,
fillColor=color_map[b_type]
).add_to(fg)
fg.add_to(m_pois)
# Add layer control to the map
folium.LayerControl().add_to(m_pois)
# Add a legend
legend_html = '''
<div style="position: fixed; bottom: 50px; left: 50px; width: 220px; height: 180px;
border:2px solid grey; z-index:9999; font-size:14px; background-color:white;
overflow-y: auto;">
<p><strong>Business Types:</strong></p>
'''
for b_type, color in color_map.items():
legend_html += f'<p><span style="color:{color};">●</span> {b_type}</p>'
legend_html += '</div>'
m_pois.get_root().html.add_child(folium.Element(legend_html))
# Save the map
m_pois.save("maps/missoula_poi_map.html")
print("POI map has been saved as 'missoula_poi_map.html'")import pandas as pd
import folium
from folium.plugins import MarkerCluster
# Load the POI data
# poi = pd.read_csv('path_to_your_poi_data.csv') # Replace with actual path to your POI data
# Create a map centered on Missoula
m_pois = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a color map for business types
unique_types = poi['b_type'].unique()
color_map = {b_type: f'#{hash(b_type) % 0xFFFFFF:06x}' for b_type in unique_types}
# Create a feature group for all POIs
all_pois = folium.FeatureGroup(name="All POIs", show=True)
# Add all POIs to the all_pois group
for _, row in poi.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=color_map[row['b_type']],
fill=True,
fillColor=color_map[row['b_type']]
).add_to(all_pois)
all_pois.add_to(m_pois)
# Create a feature group for each business type
for b_type in unique_types:
fg = folium.FeatureGroup(name=b_type, show=False)
type_df = poi[poi['b_type'] == b_type]
for _, row in type_df.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=color_map[b_type],
fill=True,
fillColor=color_map[b_type]
).add_to(fg)
fg.add_to(m_pois)
# Add layer control to the map
folium.LayerControl().add_to(m_pois)
# Add a legend
legend_html = '''
<div style="position: fixed; bottom: 50px; left: 50px; width: 220px; height: 180px;
border:2px solid grey; z-index:9999; font-size:14px; background-color:white;
overflow-y: auto;">
<p><strong>Business Types:</strong></p>
'''
for b_type, color in color_map.items():
legend_html += f'<p><span style="color:{color};">●</span> {b_type}</p>'
legend_html += '</div>'
m_pois.get_root().html.add_child(folium.Element(legend_html))
# Add custom CSS and JavaScript for the layer control label
custom_css = """
<style>
.leaflet-control-layers-list {
margin-top: 30px !important;
}
.layer-control-label {
position: absolute;
top: 10px;
left: 10px;
font-weight: bold;
}
</style>
"""
custom_js = """
<script>
document.addEventListener('DOMContentLoaded', function() {
var layerControl = document.querySelector('.leaflet-control-layers');
if (layerControl) {
var label = document.createElement('div');
label.className = 'layer-control-label';
label.textContent = 'Select POI Type';
layerControl.insertBefore(label, layerControl.firstChild);
}
});
</script>
"""
m_pois.get_root().header.add_child(folium.Element(custom_css + custom_js))
# Save the map
m_pois.save("website/maps/missoula_poi_map.html")
print("POI map has been saved as 'missoula_poi_map.html'")# load the 911 calls data into a geodataframe
calls_911_2023 = gpd.read_file('data/missoula/missoula_911Calls_2023.csv')
calls_911_2023.geometry = gpd.points_from_xy(calls_911_2023.Lon, calls_911_2023.Lat)
calls_911_2024 = gpd.read_file('data/missoula/missoula_911_calls_2024.geojson')
calls_911_2024.geometry = gpd.points_from_xy(calls_911_2024.Lon, calls_911_2024.Lat)
# merge the 911 calls data
calls_911_2023_2024 = pd.concat([calls_911_2023, calls_911_2024], ignore_index=True)
# calls_911_2023_2024.to_crs(epsg=5070, inplace=True)
# List of Violent call types
violent_calls = [
'Assault', 'Assault with a Weapon', 'Domestic Animal', 'Homicide',
'Intimidation', 'Kidnapping', 'Partner Family Member Assault',
'Robbery', 'Sexual Assault', 'Shots Fired', 'Shots Heard',
'Suicidal Person', 'Suicide', 'Weapons - Carry/Possess'
]
# Function to classify call types
def classify_call(call_type):
return 'Violent' if call_type in violent_calls else 'Nonviolent'
# Add new column 'crime_type' to the DataFrame
calls_911_2023_2024['crime_type'] = calls_911_2023_2024['Call Type Description'].apply(classify_call)# load missoula_roads_segments.geojson
missoula_roads = gpd.read_file('data/missoula/missoula_roads_clipped.geojson')
# drop roads with roadname DRIVEWAY, ALLEY
# missoula_roads = missoula_roads[~missoula_roads.roadname.isin(['DRIVEWAY', 'ALLEY', 'SERVICE'])]
missoula_roads = missoula_roads[~missoula_roads.roadname.isin(['SERVICE'])]
missoula_roads.roadname.value_counts()
missoula_roads.isnull().sum()
# drop all columns with any missing values
missoula_roads.dropna(axis=1, inplace=True)
missoula_roads.reset_index(drop=True, inplace=True)
missoula_roads['uniquename'] = missoula_roads['roadname'] + '_' + missoula_roads['fromleft'].astype(str) + '_' + missoula_roads['toleft'].astype(str) + '_' + missoula_roads.index.astype(str)
# create n meter buffer around roads
missoula_roads_buffer = missoula_roads.copy()
missoula_roads_buffer.to_crs(epsg=5070, inplace=True)
missoula_roads_buffer['geometry'] = missoula_roads_buffer.buffer(30)
missoula_roads_buffer.reset_index(drop=True, inplace=True)
missoula_roads_buffer.plot()calls_911_2023_2024_buffer = calls_911_2023_2024.copy()
calls_911_2023_2024_buffer.to_crs(epsg=5070, inplace=True)calls_911_2023_2024_buffer_violent = calls_911_2023_2024_buffer[calls_911_2023_2024_buffer['crime_type'] == 'Violent']
roads_calls_violent = missoula_roads_buffer.merge(missoula_roads_buffer.sjoin(calls_911_2023_2024_buffer_violent).groupby('uniquename').size().rename('n_points').reset_index())
roads_calls_fixed_violent = roads_calls_violent[['uniquename', 'geometry', 'n_points', 'shape_Leng']]
# remove buffer for plotting
roads_calls_fixed_violent['geometry'] = roads_calls_fixed_violent.buffer(-20)
# remove roads starting with '90' in uniquename
roads_calls_fixed_violent = roads_calls_fixed_violent[~roads_calls_fixed_violent.uniquename.str.startswith('90')]
roads_calls = missoula_roads_buffer.merge(missoula_roads_buffer.sjoin(calls_911_2023_2024_buffer).groupby('uniquename').size().rename('n_points').reset_index())
roads_calls_fixed = roads_calls[['uniquename', 'geometry', 'n_points', 'shape_Leng']]
# remove buffer for plotting
roads_calls_fixed['geometry'] = roads_calls_fixed.buffer(-20)
# remove roads starting with '90' in uniquename
roads_calls_fixed = roads_calls_fixed[~roads_calls_fixed.uniquename.str.startswith('90')]# show the roads from roads_calls_fixed with the most 911 calls
roads_calls_fixed.sort_values('n_points', ascending=False).head(10)import pandas as pd
import geopandas as gpd
import folium
from folium.features import GeoJsonPopup, GeoJsonTooltip
from branca.colormap import LinearColormap
import numpy as np
# Ensure the GeoDataFrame is in WGS84 (EPSG:4326) for Folium
road_segments = roads_calls_fixed.to_crs(epsg=4326)
# Create a map centered on Missoula with Esri.WorldStreetMap as default
m = folium.Map(location=[46.87, -113.99], zoom_start=13, min_zoom=13, max_bounds=True, tiles='Esri.WorldStreetMap', attr='Esri')
# Define the bounding box for the initial view
min_lat, max_lat, min_lon, max_lon = 46.77854245661467, 46.97166694574133, -114.17662423519394, -113.85756122668695
# Set the max bounds to restrict panning
m.fit_bounds([[min_lat, min_lon], [max_lat, max_lon]])
# Create a color map with percentile-based scaling
min_points = road_segments['n_points'].min()
max_points = road_segments['n_points'].max()
p10 = np.percentile(road_segments['n_points'], 10)
p99 = np.percentile(road_segments['n_points'], 99.5)
colormap = LinearColormap(
colors = [
'#00FF00', # Green
'#90EE90', # Light Green
'#FFA500', # Orange
'#FF8C00', # Dark Orange
'#FF4500', # Red-Orange
'#FF0000', # Red
'#DC143C', # Crimson
'#B22222', # Fire Brick
'#8B0000', # Dark Red
'#000000' # Black
],
vmin=p10,
vmax=p99
)
# Function to get color based on number of points
def get_color(n_points):
if n_points <= p10:
return colormap(p10)
elif n_points >= p99:
return colormap(p99)
else:
return colormap(n_points)
# Add road segments to the map
folium.GeoJson(
road_segments,
style_function=lambda feature: {
'fillColor': get_color(feature['properties']['n_points']),
'color': 'black',
'weight': 0,
'fillOpacity': 0.7,
},
tooltip=GeoJsonTooltip(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True),
popup=GeoJsonPopup(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True)
).add_to(m)
# Add a color legend
colormap.add_to(m)
colormap.caption = f'Number of Calls near roads (10th-99.5th percentile: {p10:.0f}-{p99:.0f})'
# Add tile layers
folium.TileLayer(
tiles='Esri.WorldImagery',
attr='Esri',
name='Esri.WorldImagery',
overlay=False,
control=True
).add_to(m)
folium.TileLayer(
tiles='Esri.WorldGrayCanvas',
attr='Esri',
name='Esri.WorldGrayCanvas',
overlay=False,
control=True
).add_to(m)
# Add layer control
folium.LayerControl().add_to(m)
# Save the map
m.save("website/maps/missoula_road_segments_map.html")
print("Map of road segments has been saved as 'missoula_road_segments_map.html'")import pandas as pd
import geopandas as gpd
import folium
from folium.features import GeoJsonPopup, GeoJsonTooltip
from branca.colormap import LinearColormap
import numpy as np
# Ensure the GeoDataFrame is in WGS84 (EPSG:4326) for Folium
road_segments = roads_calls_fixed.to_crs(epsg=4326)
# Create a map centered on Missoula with Esri.WorldStreetMap as default
m = folium.Map(location=[46.87, -113.99], zoom_start=12, min_zoom=12, max_bounds=True, tiles='Esri.WorldStreetMap', attr='Esri')
# Define the bounding box for the initial view
min_lat, max_lat, min_lon, max_lon = 46.77854245661467, 46.97166694574133, -114.17662423519394, -113.85756122668695
# Set the max bounds to restrict panning
m.fit_bounds([[min_lat, min_lon], [max_lat, max_lon]])
# Calculate mean and standard deviation
mean_points = road_segments['n_points'].mean()
std_points = road_segments['n_points'].std()
# Define color scale range (mean +/- 2 standard deviations)
vmin = (mean_points - 1 * std_points) if (mean_points - 0.5 * std_points) > 0 else 0
vmax = mean_points + 2.0 * std_points
# Create a diverging color map
colormap = LinearColormap(
colors = [
'#440154', '#482878', '#3E4989', '#31688E', '#26828E',
'#1F9E89', '#35B779', '#6DCD59', '#B4DE2C', '#FDE725'
],
vmin=vmin,
vmax=vmax,
caption=f'Number of Points (Mean: {mean_points:.0f}, Std: {std_points:.0f})'
)
# Function to get color based on number of points
def get_color(n_points):
if n_points < vmin:
return colormap(vmin)
elif n_points > vmax:
return colormap(vmax)
else:
return colormap(n_points)
# Add road segments to the map
folium.GeoJson(
road_segments,
style_function=lambda feature: {
'fillColor': get_color(feature['properties']['n_points']),
'color': 'black',
'weight': 0,
'fillOpacity': 0.7,
},
tooltip=GeoJsonTooltip(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True),
popup=GeoJsonPopup(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True)
).add_to(m)
# Add the color scale to the map
colormap.add_to(m)
# Add tile layers
folium.TileLayer(
tiles='Esri.WorldImagery',
attr='Esri',
name='Esri.WorldImagery',
overlay=False,
control=True
).add_to(m)
folium.TileLayer(
tiles='Esri.WorldGrayCanvas',
attr='Esri',
name='Esri.WorldGrayCanvas',
overlay=False,
control=True
).add_to(m)
# Add layer control
folium.LayerControl().add_to(m)
# Save the map
m.save("maps/missoula_road_segments_map_std_scaled.html")
print("Map of road segments has been saved as 'missoula_road_segments_map_std_scaled.html'")import pandas as pd
import geopandas as gpd
import folium
from folium.features import GeoJsonPopup, GeoJsonTooltip
from branca.colormap import LinearColormap
import numpy as np
# Ensure the GeoDataFrame is in WGS84 (EPSG:4326) for Folium
road_segments = roads_calls_fixed_violent.to_crs(epsg=4326)
# Create a map centered on Missoula with Esri.WorldStreetMap as default
m = folium.Map(location=[46.87, -113.99], zoom_start=13, min_zoom=13, max_bounds=True, tiles='Esri.WorldStreetMap', attr='Esri')
# Define the bounding box for the initial view
min_lat, max_lat, min_lon, max_lon = 46.77854245661467, 46.97166694574133, -114.17662423519394, -113.85756122668695
# Set the max bounds to restrict panning
m.fit_bounds([[min_lat, min_lon], [max_lat, max_lon]])
# Create a color map with percentile-based scaling
min_points = road_segments['n_points'].min()
max_points = road_segments['n_points'].max()
p10 = np.percentile(road_segments['n_points'], 10)
p99 = np.percentile(road_segments['n_points'], 99.5)
colormap = LinearColormap(
colors = [
'#00FF00', # Green
'#90EE90', # Light Green
'#FFA500', # Orange
'#FF8C00', # Dark Orange
'#FF4500', # Red-Orange
'#FF0000', # Red
'#DC143C', # Crimson
'#B22222', # Fire Brick
'#8B0000', # Dark Red
'#000000' # Black
],
vmin=p10,
vmax=p99
)
# Function to get color based on number of points
def get_color(n_points):
if n_points <= p10:
return colormap(p10)
elif n_points >= p99:
return colormap(p99)
else:
return colormap(n_points)
# Add road segments to the map
folium.GeoJson(
road_segments,
style_function=lambda feature: {
'fillColor': get_color(feature['properties']['n_points']),
'color': 'black',
'weight': 0,
'fillOpacity': 0.7,
},
tooltip=GeoJsonTooltip(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True),
popup=GeoJsonPopup(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True)
).add_to(m)
# Add a color legend
colormap.add_to(m)
colormap.caption = f'Number of violent calls near roads (10th-99.5th percentile: {p10:.0f}-{p99:.0f})'
# Add tile layers
folium.TileLayer(
tiles='Esri.WorldImagery',
attr='Esri',
name='Esri.WorldImagery',
overlay=False,
control=True
).add_to(m)
folium.TileLayer(
tiles='Esri.WorldGrayCanvas',
attr='Esri',
name='Esri.WorldGrayCanvas',
overlay=False,
control=True
).add_to(m)
# Add layer control
folium.LayerControl().add_to(m)
# Save the map
m.save("website/maps/missoula_road_segments_violent_map.html")
print("Map of road segments has been saved as 'missoula_road_segments_violent_map.html'")import pandas as pd
import geopandas as gpd
import folium
from folium.features import GeoJsonPopup, GeoJsonTooltip
from branca.colormap import LinearColormap
import numpy as np
from shapely.ops import unary_union
# Ensure the data is in EPSG:5070 projection
missoula_roads = missoula_roads.to_crs(epsg=5070)
# Create a smart buffer that fills all space between polygons
def create_smart_buffer(gdf):
buffer_distance = 1
buffered = gdf.buffer(buffer_distance)
while not unary_union(buffered).is_valid:
buffer_distance += 1
buffered = gdf.buffer(buffer_distance)
return gpd.GeoDataFrame(geometry=buffered, crs=gdf.crs)
# Create the smart buffer
missoula_roads_buffer = create_smart_buffer(missoula_roads)
missoula_roads_buffer['uniquename'] = missoula_roads['uniquename']
# Process 911 calls data
calls_911_2023_2024_buffer = calls_911_2023_2024.copy()
calls_911_2023_2024_buffer.to_crs(epsg=5070, inplace=True)
# Perform spatial join and count points
roads_calls = missoula_roads_buffer.sjoin(calls_911_2023_2024_buffer).groupby('uniquename').size().rename('n_points').reset_index()
roads_calls = missoula_roads_buffer.merge(roads_calls, on='uniquename', how='left')
# Fill NaN values with 0 for roads with no calls
roads_calls['n_points'] = roads_calls['n_points'].fillna(0)
# Keep only necessary columns
roads_calls_fixed = roads_calls[['uniquename', 'geometry', 'n_points']]
# Remove roads starting with '90' in uniquename
roads_calls_fixed = roads_calls_fixed[~roads_calls_fixed.uniquename.str.startswith('90')]
# Create a 20m buffer for visualization
roads_calls_fixed['geometry'] = roads_calls_fixed.buffer(20)
# Ensure the GeoDataFrame is in WGS84 (EPSG:4326) for Folium
road_segments = roads_calls_fixed.to_crs(epsg=4326)
# Create a map centered on Missoula
m = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a color map with percentile-based scaling
min_points = road_segments['n_points'].min()
max_points = road_segments['n_points'].max()
p10 = np.percentile(road_segments['n_points'], 10)
p99 = np.percentile(road_segments['n_points'], 99)
colormap = LinearColormap(
colors = [
'#00FF00', # Green
'#90EE90', # Light Green
'#FFA500', # Orange
'#FF8C00', # Dark Orange
'#FF4500', # Red-Orange
'#FF0000', # Red
'#DC143C', # Crimson
'#B22222', # Fire Brick
'#8B0000', # Dark Red
'#000000' # Black
],
vmin=p10,
vmax=p99
)
# Function to get color based on number of points
def get_color(n_points):
if n_points <= p10:
return colormap(p10)
elif n_points >= p99:
return colormap(p99)
else:
return colormap(n_points)
# Add road segments to the map
folium.GeoJson(
road_segments,
style_function=lambda feature: {
'fillColor': get_color(feature['properties']['n_points']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7,
},
tooltip=GeoJsonTooltip(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True),
popup=GeoJsonPopup(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True)
).add_to(m)
# Add a color legend
colormap.add_to(m)
colormap.caption = f'Number of Points on Road Segment (10th-90th percentile: {p10:.0f}-{p99:.0f})'
# Save the map
m.save("maps/missoula_road_segments_map_smartbuffer.html")
print("Map of road segments has been saved as 'missoula_road_segments_map_smartbuffer.html'")import pandas as pd
import geopandas as gpd
import folium
from folium.features import GeoJsonPopup, GeoJsonTooltip
from branca.colormap import LinearColormap
import numpy as np
from shapely.ops import unary_union
# Ensure the data is in EPSG:5070 projection
missoula_roads = missoula_roads.to_crs(epsg=5070)
# Create a smart buffer that fills all space between polygons
def create_smart_buffer(gdf):
buffer_distance = 1
buffered = gdf.buffer(buffer_distance)
while not unary_union(buffered).is_valid:
buffer_distance += 1
buffered = gdf.buffer(buffer_distance)
return gpd.GeoDataFrame(geometry=buffered, crs=gdf.crs)
# Create the smart buffer
missoula_roads_buffer = create_smart_buffer(missoula_roads)
missoula_roads_buffer['uniquename'] = missoula_roads['uniquename']
# Process 911 calls data
calls_911_2023_2024_buffer = calls_911_2023_2024.copy()
calls_911_2023_2024_buffer.to_crs(epsg=5070, inplace=True)
# Perform spatial join and count points
roads_calls = missoula_roads_buffer.sjoin(calls_911_2023_2024_buffer).groupby('uniquename').size().rename('n_points').reset_index()
roads_calls = missoula_roads_buffer.merge(roads_calls, on='uniquename', how='left')
# Fill NaN values with 0 for roads with no calls
roads_calls['n_points'] = roads_calls['n_points'].fillna(0)
# Keep only necessary columns
roads_calls_fixed = roads_calls[['uniquename', 'geometry', 'n_points']]
# Remove roads starting with '90' in uniquename
roads_calls_fixed = roads_calls_fixed[~roads_calls_fixed.uniquename.str.startswith('90')]
# Create a copy for the 20m buffer
roads_calls_20m = roads_calls_fixed.copy()
roads_calls_20m['geometry'] = roads_calls_20m.buffer(20)
# Ensure the GeoDataFrames are in WGS84 (EPSG:4326) for Folium
road_segments_original = roads_calls_fixed.to_crs(epsg=4326)
road_segments_20m = roads_calls_20m.to_crs(epsg=4326)
# Create a map centered on Missoula
m = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a color map with percentile-based scaling
min_points = road_segments_original['n_points'].min()
max_points = road_segments_original['n_points'].max()
p10 = np.percentile(road_segments_original['n_points'], 10)
p99 = np.percentile(road_segments_original['n_points'], 99)
colormap = LinearColormap(
colors = [
'#00FF00', # Green
'#90EE90', # Light Green
'#FFA500', # Orange
'#FF8C00', # Dark Orange
'#FF4500', # Red-Orange
'#FF0000', # Red
'#DC143C', # Crimson
'#B22222', # Fire Brick
'#8B0000', # Dark Red
'#000000' # Black
],
vmin=p10,
vmax=p99
)
# Function to get color based on number of points
def get_color(n_points):
if n_points <= p10:
return colormap(p10)
elif n_points >= p99:
return colormap(p99)
else:
return colormap(n_points)
# Add original road segments to the map
folium.GeoJson(
road_segments_original,
style_function=lambda feature: {
'fillColor': get_color(feature['properties']['n_points']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7,
},
tooltip=GeoJsonTooltip(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True),
popup=GeoJsonPopup(fields=['uniquename', 'n_points'],
aliases=['Road Name', 'Number of Points'],
localize=True)
).add_to(m)
# Add 20m buffered road segments to the map
folium.GeoJson(
road_segments_20m,
style_function=lambda feature: {
'fillColor': 'none',
'color': 'white',
'weight': 1,
'fillOpacity': 0,
}
).add_to(m)
# Add a color legend
colormap.add_to(m)
colormap.caption = f'Number of Points on Road Segment (10th-90th percentile: {p10:.0f}-{p99:.0f})'
# Save the map
m.save("maps/missoula_road_segments_dual_buffer_map.html")
print("Map of road segments with dual buffers has been saved as 'missoula_road_segments_dual_buffer_map.html'")import pandas as pd
import geopandas as gpd
import folium
from folium.plugins import Draw
import random
import json
from branca.element import MacroElement, Template
# Load the data
gdf = calls_911_2023_2024 # Assuming this is already a GeoDataFrame
# Function to generate random color
def random_color():
return f'#{random.randint(0, 0xFFFFFF):06x}'
# Create a color map with random colors
unique_call_types = gdf['Call Type Description'].unique()
color_dict = {call_type: random_color() for call_type in unique_call_types}
# Convert GeoDataFrame to a list of dictionaries for JSON serialization
calls_data = gdf.copy()
calls_data['lat'] = calls_data.geometry.y
calls_data['lon'] = calls_data.geometry.x
calls_data = calls_data.drop(columns=['geometry']).to_dict('records')
# Create the map
m_911_radius = folium.Map(location=[46.87, -114.00], zoom_start=12)
# Add Draw plugin
draw = Draw()
draw.add_to(m_911_radius)
# Custom MacroElement to add JavaScript
class ClickForMarker(MacroElement):
def __init__(self, calls_data, color_dict):
super(ClickForMarker, self).__init__()
self._template = Template("""
{% macro script(this, kwargs) %}
var calls_data = {{ this.calls_data }};
var color_dict = {{ this.color_dict }};
var currentMarker = null;
function updateMap(lat, lon) {
// Clear previous update results
{{ this._parent.get_name() }}.eachLayer(function(layer) {
if (layer instanceof L.CircleMarker || layer instanceof L.Circle) {
{{ this._parent.get_name() }}.removeLayer(layer);
}
});
// Remove the current marker if it exists
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
// Create a 100m buffer around the clicked point (approximate)
var buffer_size = 0.00090; // Approx. 100m in degrees
// Filter calls within the buffer
var calls_in_buffer = calls_data.filter(function(call) {
var dx = call.lon - lon;
var dy = call.lat - lat;
return Math.sqrt(dx*dx + dy*dy) <= buffer_size;
});
// Add circle to show buffer
L.circle([lat, lon], {
radius: 100,
color: 'red',
fillColor: 'red',
fillOpacity: 0.2
}).addTo({{ this._parent.get_name() }});
// Add points for calls
calls_in_buffer.forEach(function(call) {
L.circleMarker([call.lat, call.lon], {
radius: 5,
color: color_dict[call['Call Type Description']],
fillColor: color_dict[call['Call Type Description']],
fillOpacity: 1
}).bindPopup(call['Call Type Description']).addTo({{ this._parent.get_name() }});
});
// Create popup with top 10 call types
var call_type_counts = {};
calls_in_buffer.forEach(function(call) {
var type = call['Call Type Description'];
call_type_counts[type] = (call_type_counts[type] || 0) + 1;
});
var top_calls = Object.entries(call_type_counts)
.sort((a, b) => b[1] - a[1])
.slice(0, 10);
var popup_html = "<h4>Top Call Types:</h4>";
top_calls.forEach(function([call_type, count]) {
popup_html += call_type + ": " + count + "<br>";
});
L.popup()
.setLatLng([lat, lon])
.setContent(popup_html)
.openOn({{ this._parent.get_name() }});
}
{{ this._parent.get_name() }}.on('click', function(e) {
var lat = e.latlng.lat.toFixed(6);
var lon = e.latlng.lng.toFixed(6);
// Remove the current marker if it exists
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
// Create a new marker and store it in currentMarker
currentMarker = L.marker([lat, lon]).addTo({{ this._parent.get_name() }});
currentMarker.bindPopup('<button onclick="updateMap(' + lat + ',' + lon + ')">Analyze This Area</button>').openPopup();
});
{% endmacro %}
""")
self.calls_data = calls_data
self.color_dict = color_dict
# Add the custom MacroElement to the map
m_911_radius.add_child(ClickForMarker(json.dumps(calls_data), json.dumps(color_dict)))
# Save the map
m_911_radius.save("maps/911_calls_map.html")
print("Map saved as 911_calls_map.html")import pandas as pd
import geopandas as gpd
import folium
from folium.plugins import Draw
import random
import json
from branca.element import MacroElement, Template
from geopy.geocoders import Nominatim
from geopy.exc import GeocoderTimedOut
# Load the data
gdf = calls_911_2023_2024 # Assuming this is already a GeoDataFrame
# Function to generate random color
def random_color():
return f'#{random.randint(0, 0xFFFFFF):06x}'
# Create a color map with random colors
unique_call_types = gdf['Call Type Description'].unique()
color_dict = {call_type: random_color() for call_type in unique_call_types}
# Convert GeoDataFrame to a list of dictionaries for JSON serialization
calls_data = gdf.copy()
calls_data['lat'] = calls_data.geometry.y
calls_data['lon'] = calls_data.geometry.x
calls_data = calls_data.drop(columns=['geometry']).to_dict('records')
# Create the map
m_911_radius = folium.Map(location=[46.87, -114.00], zoom_start=12)
# # Add Draw plugin
# draw = Draw()
# draw.add_to(m_911_radius)
# Custom MacroElement to add JavaScript
class ClickForMarker(MacroElement):
def __init__(self, calls_data, color_dict):
super(ClickForMarker, self).__init__()
self._template = Template("""
{% macro script(this, kwargs) %}
var calls_data = {{ this.calls_data }};
var color_dict = {{ this.color_dict }};
var currentMarker = null;
var currentCircle = null;
function metersToDegreesLat(meters) {
return meters / 111320.0; // Roughly 111,320 meters per degree of latitude
}
function metersToDegreesLon(meters, lat) {
return meters / (111320.0 * Math.cos(lat * Math.PI / 180));
}
function updateMap(lat, lon) {
// Get the radius in meters from the input
var radiusMeters = parseFloat(document.getElementById('radius-input').value) || 100;
var radiusLat = metersToDegreesLat(radiusMeters);
var radiusLon = metersToDegreesLon(radiusMeters, lat);
// Clear previous update results
{{ this._parent.get_name() }}.eachLayer(function(layer) {
if (layer instanceof L.CircleMarker || layer instanceof L.Circle) {
{{ this._parent.get_name() }}.removeLayer(layer);
}
});
// Remove the current marker if it exists
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
// Remove the current circle if it exists
if (currentCircle) {
{{ this._parent.get_name() }}.removeLayer(currentCircle);
}
// Filter calls within the buffer
var calls_in_buffer = calls_data.filter(function(call) {
var dx = call.lon - lon;
var dy = call.lat - lat;
return Math.sqrt(dx*dx + dy*dy) <= Math.max(radiusLat, radiusLon);
});
// Add circle to show buffer
currentCircle = L.circle([lat, lon], {
radius: radiusMeters,
color: 'red',
fillColor: 'red',
fillOpacity: 0.2
}).addTo({{ this._parent.get_name() }});
// Add points for calls
calls_in_buffer.forEach(function(call) {
L.circleMarker([call.lat, call.lon], {
radius: 5,
color: color_dict[call['Call Type Description']],
fillColor: color_dict[call['Call Type Description']],
fillOpacity: 1
}).bindPopup(call['Call Type Description']).addTo({{ this._parent.get_name() }});
});
// Create popup with top 10 call types
var call_type_counts = {};
calls_in_buffer.forEach(function(call) {
var type = call['Call Type Description'];
call_type_counts[type] = (call_type_counts[type] || 0) + 1;
});
var top_calls = Object.entries(call_type_counts)
.sort((a, b) => b[1] - a[1])
.slice(0, 10);
var popup_html = "<h4>Top Call Types:</h4>";
top_calls.forEach(function([call_type, count]) {
popup_html += call_type + ": " + count + "<br>";
});
L.popup()
.setLatLng([lat, lon])
.setContent(popup_html)
.openOn({{ this._parent.get_name() }});
}
{{ this._parent.get_name() }}.on('click', function(e) {
var lat = e.latlng.lat.toFixed(6);
var lon = e.latlng.lng.toFixed(6);
// Remove the current marker if it exists
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
// Create a new marker and store it in currentMarker
currentMarker = L.marker([lat, lon]).addTo({{ this._parent.get_name() }});
currentMarker.bindPopup('<button onclick="updateMap(' + lat + ',' + lon + ')">Analyze This Area</button>').openPopup();
});
// Add search and radius control
var searchControl = L.Control.extend({
options: {
position: 'topright'
},
onAdd: function (map) {
var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control leaflet-control-custom');
container.style.backgroundColor = 'white';
container.style.padding = '5px';
container.innerHTML = '<input type="text" id="address-input" placeholder="Enter address in Missoula, MT">' +
'<input type="number" id="radius-input" placeholder="Radius (meters)" value="100">' +
'<button id="search-button">Search</button>';
container.style.fontSize = '14px';
return container;
}
});
{{ this._parent.get_name() }}.addControl(new searchControl());
document.getElementById('search-button').addEventListener('click', function() {
var address = document.getElementById('address-input').value + ', Missoula, MT';
fetch('https://nominatim.openstreetmap.org/search?format=json&q=' + encodeURIComponent(address))
.then(response => response.json())
.then(data => {
if (data.length > 0) {
var lat = parseFloat(data[0].lat);
var lon = parseFloat(data[0].lon);
{{ this._parent.get_name() }}.setView([lat, lon], 16);
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
currentMarker = L.marker([lat, lon]).addTo({{ this._parent.get_name() }});
currentMarker.bindPopup('<button onclick="updateMap(' + lat + ',' + lon + ')">Analyze This Area</button>').openPopup();
} else {
alert('Address not found');
}
})
.catch(error => {
console.error('Error:', error);
alert('An error occurred while searching for the address');
});
});
{% endmacro %}
""")
self.calls_data = calls_data
self.color_dict = color_dict
# Add the custom MacroElement to the map
m_911_radius.add_child(ClickForMarker(json.dumps(calls_data), json.dumps(color_dict)))
# Save the map
m_911_radius.save("9maps/11_calls_map_2.html")
print("Map saved as 911_calls_map_2.html")import pandas as pd
import geopandas as gpd
import folium
from folium.plugins import Draw
import random
import json
from branca.element import MacroElement, Template
# Load the data
gdf = calls_911_2023_2024 # Assuming this is already a GeoDataFrame
poi = poi.copy()
# Function to generate random color
def random_color():
return f'#{random.randint(0, 0xFFFFFF):06x}'
# Create color maps
unique_call_types = gdf['Call Type Description'].unique()
call_color_dict = {call_type: random_color() for call_type in unique_call_types}
unique_poi_types = poi['b_type'].unique()
poi_color_map = {b_type: f'#{hash(b_type) % 0xFFFFFF:06x}' for b_type in unique_poi_types}
# Convert GeoDataFrame to a list of dictionaries for JSON serialization
calls_data = gdf.copy()
calls_data['lat'] = calls_data.geometry.y
calls_data['lon'] = calls_data.geometry.x
calls_data = calls_data.drop(columns=['geometry']).to_dict('records')
# Create the map
m_combined = folium.Map(location=[46.87, -114.00], zoom_start=12)
# Add Draw plugin
draw = Draw()
draw.add_to(m_combined)
# Create feature groups for POIs
all_pois = folium.FeatureGroup(name="All POIs", show=False)
for _, row in poi.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=poi_color_map[row['b_type']],
fill=True,
fillColor=poi_color_map[row['b_type']]
).add_to(all_pois)
all_pois.add_to(m_combined)
for b_type in unique_poi_types:
fg = folium.FeatureGroup(name=f"POI: {b_type}", show=False)
type_df = poi[poi['b_type'] == b_type]
for _, row in type_df.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=5,
popup=f"{row['b_type']}",
color=poi_color_map[b_type],
fill=True,
fillColor=poi_color_map[b_type]
).add_to(fg)
fg.add_to(m_combined)
# Custom MacroElement to add JavaScript
class ClickForMarker(MacroElement):
def __init__(self, calls_data, call_color_dict):
super(ClickForMarker, self).__init__()
self._template = Template("""
{% macro script(this, kwargs) %}
var calls_data = {{ this.calls_data }};
var call_color_dict = {{ this.call_color_dict }};
var currentMarker = null;
var currentCircle = null;
var currentCallMarkers = L.layerGroup().addTo({{ this._parent.get_name() }});
function metersToDegreesLat(meters) {
return meters / 111320.0;
}
function metersToDegreesLon(meters, lat) {
return meters / (111320.0 * Math.cos(lat * Math.PI / 180));
}
function updateMap(lat, lon) {
var radiusMeters = parseInt(document.getElementById('radius-select').value);
var radiusLat = metersToDegreesLat(radiusMeters);
var radiusLon = metersToDegreesLon(radiusMeters, lat);
// change radiusLat and radiusLon based on radiusMeters (i.e. var buffer_size = 0.00090; Approx. 100m in degrees)
if (radiusMeters == 10) {
radiusLat = 0.00010;
radiusLon = 0.00010;
} else if (radiusMeters == 25) {
radiusLat = 0.00025;
radiusLon = 0.00025;
} else if (radiusMeters == 50) {
radiusLat = 0.00050;
radiusLon = 0.00050;
} else if (radiusMeters == 100) {
radiusLat = 0.00100;
radiusLon = 0.00100;
} else if (radiusMeters == 200) {
radiusLat = 0.00200;
radiusLon = 0.00200;;
}
currentCallMarkers.clearLayers();
if (currentCircle) {
{{ this._parent.get_name() }}.removeLayer(currentCircle);
}
currentCircle = L.circle([lat, lon], {
radius: radiusMeters,
color: 'red',
fillColor: 'red',
fillOpacity: 0.2
}).addTo({{ this._parent.get_name() }});
var calls_in_buffer = calls_data.filter(function(call) {
var dx = call.lon - lon;
var dy = call.lat - lat;
return Math.sqrt(dx*dx + dy*dy) <= Math.max(radiusLat, radiusLon);
});
calls_in_buffer.forEach(function(call) {
L.circleMarker([call.lat, call.lon], {
radius: 5,
color: call_color_dict[call['Call Type Description']],
fillColor: call_color_dict[call['Call Type Description']],
fillOpacity: 1
}).bindPopup(call['Call Type Description']).addTo(currentCallMarkers);
});
var call_type_counts = {};
calls_in_buffer.forEach(function(call) {
var type = call['Call Type Description'];
call_type_counts[type] = (call_type_counts[type] || 0) + 1;
});
var top_calls = Object.entries(call_type_counts)
.sort((a, b) => b[1] - a[1])
.slice(0, 10);
var popup_html = "<h4>Top Call Types:</h4>";
top_calls.forEach(function([call_type, count]) {
popup_html += call_type + ": " + count + "<br>";
});
L.popup()
.setLatLng([lat, lon])
.setContent(popup_html)
.openOn({{ this._parent.get_name() }});
}
{{ this._parent.get_name() }}.on('click', function(e) {
var lat = e.latlng.lat.toFixed(6);
var lon = e.latlng.lng.toFixed(6);
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
currentMarker = L.marker([lat, lon]).addTo({{ this._parent.get_name() }});
currentMarker.bindPopup('<button onclick="updateMap(' + lat + ',' + lon + ')">Analyze This Area</button>').openPopup();
});
var searchControl = L.Control.extend({
options: {
position: 'topright'
},
onAdd: function (map) {
var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control leaflet-control-custom');
container.style.backgroundColor = 'white';
container.style.padding = '5px';
container.innerHTML = '<input type="text" id="address-input" placeholder="Enter address in Missoula, MT">' +
'<select id="radius-select">' +
'<option value="10">10m</option>' +
'<option value="25">25m</option>' +
'<option value="50">50m</option>' +
'<option value="100" selected>100m</option>' +
'<option value="200">200m</option>' +
'</select>' +
'<button id="search-button">Search</button>';
container.style.fontSize = '14px';
return container;
}
});
{{ this._parent.get_name() }}.addControl(new searchControl());
document.getElementById('search-button').addEventListener('click', function() {
var address = document.getElementById('address-input').value + ', Missoula, MT';
fetch('https://nominatim.openstreetmap.org/search?format=json&q=' + encodeURIComponent(address))
.then(response => response.json())
.then(data => {
if (data.length > 0) {
var lat = parseFloat(data[0].lat);
var lon = parseFloat(data[0].lon);
{{ this._parent.get_name() }}.setView([lat, lon], 16);
if (currentMarker) {
{{ this._parent.get_name() }}.removeLayer(currentMarker);
}
currentMarker = L.marker([lat, lon]).addTo({{ this._parent.get_name() }});
currentMarker.bindPopup('<button onclick="updateMap(' + lat + ',' + lon + ')">Analyze This Area</button>').openPopup();
} else {
alert('Address not found');
}
})
.catch(error => {
console.error('Error:', error);
alert('An error occurred while searching for the address');
});
});
{% endmacro %}
""")
self.calls_data = calls_data
self.call_color_dict = call_color_dict
# Add the custom MacroElement to the map
m_combined.add_child(ClickForMarker(json.dumps(calls_data), json.dumps(call_color_dict)))
# Add layer control to the map
folium.LayerControl().add_to(m_combined)
# Add a legend for POIs
poi_legend_html = '''
<div style="position: fixed; bottom: 50px; left: 50px; width: 220px; height: 180px;
border:2px solid grey; z-index:9999; font-size:14px; background-color:white;
overflow-y: auto;">
<p><strong>POI Types:</strong></p>
'''
for b_type, color in poi_color_map.items():
poi_legend_html += f'<p><span style="color:{color};">●</span> {b_type}</p>'
poi_legend_html += '</div>'
m_combined.get_root().html.add_child(folium.Element(poi_legend_html))
# Add a legend for 911 call types
call_legend_html = '''
<div style="position: fixed; bottom: 50px; right: 50px; width: 220px; height: 180px;
border:2px solid grey; z-index:9999; font-size:14px; background-color:white;
overflow-y: auto;">
<p><strong>911 Call Types:</strong></p>
'''
for call_type, color in call_color_dict.items():
call_legend_html += f'<p><span style="color:{color};">●</span> {call_type}</p>'
call_legend_html += '</div>'
m_combined.get_root().html.add_child(folium.Element(call_legend_html))
# Save the map
m_combined.save("maps/911_poi_combined_map.html")
print("Combined map saved as 911_poi_combined_map.html")buffer_radius = 50import pandas as pd
import geopandas as gpd
# Assuming calls_911_2023_2024 and poi are already GeoDataFrames in CRS 5070
calls_911_gdf = calls_911_2023_2024.to_crs(epsg=5070)
poi_gdf = poi
# Ensure the geometry columns are properly set
calls_911_gdf = calls_911_gdf.set_geometry('geometry')
poi_gdf = poi_gdf.set_geometry('geometry')
# Count the number of each b_type
b_type_counts = poi_gdf['b_type'].value_counts().reset_index()
b_type_counts.columns = ['b_type', 'POI_Count']
# Create a buffer around each POI
poi_gdf['buffer'] = poi_gdf.geometry.buffer(buffer_radius)
poi_gdf = poi_gdf.set_geometry('buffer')
# Spatial join to find calls within 50 meters of each POI
calls_near_poi = gpd.sjoin(calls_911_gdf, poi_gdf[['buffer', 'b_type']], predicate='within', how='inner')
# Group by b_type and Call Type Description
grouped = calls_near_poi.groupby(['b_type', 'Call Type Description']).size().unstack(fill_value=0)
# Calculate total for each b_type and add as a column
grouped['Total'] = grouped.sum(axis=1)
# Sort rows by total calls descending
poi_911_calls = grouped.sort_values('Total', ascending=False)
# Sort columns by total calls descending, keeping 'Total' as the last column
column_order = poi_911_calls.sum().sort_values(ascending=False).index.tolist()
column_order.remove('Total')
column_order.append('Total')
poi_911_calls = poi_911_calls[column_order]
# Add the POI count column
poi_911_calls = poi_911_calls.reset_index().merge(b_type_counts, on='b_type', how='left').set_index('b_type')
# Reorder columns to put POI_Count right after b_type
columns = poi_911_calls.columns.tolist()
columns.remove('POI_Count')
columns.insert(0, 'POI_Count')
poi_911_calls_by_b_type = poi_911_calls[columns]
# Rename the index to 'b_type'
poi_911_calls_by_b_type.index.name = 'b_type'
print(poi_911_calls_by_b_type)
# Save to CSV
poi_911_calls_by_b_type.to_csv('data/missoula/poi_911_calls_summary.csv')import pandas as pd
import geopandas as gpd
# Assuming calls_911_2023_2024 and poi are already GeoDataFrames in CRS 5070
calls_911_gdf = calls_911_2023_2024.to_crs(epsg=5070)
poi_gdf = poi
# Ensure the geometry columns are properly set
calls_911_gdf = calls_911_gdf.set_geometry('geometry')
poi_gdf = poi_gdf.set_geometry('geometry')
# Create a buffer around each POI
poi_gdf['buffer'] = poi_gdf.geometry.buffer(buffer_radius)
poi_gdf = poi_gdf.set_geometry('buffer')
# Combine b_types for POIs with the same address
poi_grouped = poi_gdf.groupby('formatted').agg({
'b_type': lambda x: ', '.join(sorted(set(x))),
'buffer': 'first' # Take the first buffer for each address
}).reset_index()
# set it as a gpd
poi_grouped = gpd.GeoDataFrame(poi_grouped, geometry='buffer')
poi_grouped.crs = poi_gdf.crs
# Spatial join to find calls within 50 meters of each POI
calls_near_poi = gpd.sjoin(calls_911_gdf, poi_grouped[['buffer', 'formatted', 'b_type']], predicate='within', how='inner')
# Group by formatted address and Call Type Description
grouped = calls_near_poi.groupby(['formatted', 'Call Type Description']).size().unstack(fill_value=0)
# Calculate total for each address and add as a column
grouped['Total'] = grouped.sum(axis=1)
# Sort rows by total calls descending
poi_911_calls = grouped.sort_values('Total', ascending=False)
# Sort columns by total calls descending, keeping 'Total' as the last column
column_order = poi_911_calls.sum().sort_values(ascending=False).index.tolist()
column_order.remove('Total')
column_order.append('Total')
poi_911_calls = poi_911_calls[column_order]
# Add the b_type column
b_type_map = poi_grouped.set_index('formatted')['b_type']
poi_911_calls['b_type'] = poi_911_calls.index.map(b_type_map)
# Reorder columns to put b_type right after the index
columns = poi_911_calls.columns.tolist()
columns.remove('b_type')
columns.insert(0, 'b_type')
poi_911_calls_by_poi_address = poi_911_calls[columns]
# Rename the index to 'Address'
poi_911_calls_by_poi_address.index.name = 'Address'
print(poi_911_calls_by_poi_address)
# Save to CSV
poi_911_calls_by_poi_address.to_csv('data/missoula/poi_911_calls_summary_by_address.csv')poi_911_calls_by_poi_address# plot poi_911_calls for each b_type, sum across all call types
poi_911_calls_by_b_type.loc[:, 'Total'].sort_values(ascending=False).plot(kind='bar', figsize=(10, 8), color='skyblue')
plt.title('Total 911 Calls by Business Type within {} Meters'.format(buffer_radius))
plt.ylabel('Total Calls')
plt.xlabel('Business Type')
plt.show()# plot poi_911_calls for each b_type, sum across all call types, scaled by the number of POIs of that type
poi_911_calls_by_b_type['Total_per_POI'] = poi_911_calls_by_b_type['Total'] / poi_911_calls_by_b_type['POI_Count']
poi_911_calls_by_b_type.loc[:, 'Total_per_POI'].sort_values(ascending=False).plot(kind='bar', figsize=(14, 10), color='skyblue')
plt.title('Total 911 Calls by Business Type per POI, scaled by number of POIs within {} Meters'.format(buffer_radius))
plt.ylabel('Total Calls per POI')
plt.xlabel('Business Type')
plt.show()import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# Load the data (assuming you've already created this CSV in the previous step)
poi_911_calls = pd.read_csv('data/missoula/poi_911_calls_summary.csv', index_col='b_type')
# Remove 'POI_Count' and 'Total' columns for this analysis
call_type_data = poi_911_calls.drop(['POI_Count', 'Total'], axis=1)
# Create dropdown menu options
dropdown_options = [{'label': b_type, 'value': b_type} for b_type in call_type_data.index]
# Function to get top 10 call types for a given b_type
def get_top_10_call_types(b_type):
return call_type_data.loc[b_type].nlargest(10)
# Create the initial plot (we'll use the first b_type as default)
initial_b_type = call_type_data.index[0]
initial_data = get_top_10_call_types(initial_b_type)
fig = make_subplots(rows=1, cols=1)
bar = go.Bar(
x=initial_data.index,
y=initial_data.values,
name=initial_b_type
)
fig.add_trace(bar)
# Update layout
fig.update_layout(
title='Top 10 Call Types by POI Type within {} meters of POIs'.format(buffer_radius),
xaxis_title='Call Type',
yaxis_title='Number of Calls',
updatemenus=[{
'buttons': [
{
'method': 'update',
'label': b_type,
'args': [{'x': [get_top_10_call_types(b_type).index],
'y': [get_top_10_call_types(b_type).values]},
{'title': f'Top 10 Call Types for {b_type} within {buffer_radius} meters of POIs',}]
} for b_type in call_type_data.index
],
'direction': 'down',
'showactive': True,
}]
)
# Show the plot
fig.show()
# If you want to save the plot as an HTML file for later viewing:
fig.write_html("charts/interactive_911_calls_plot.html")import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# Load the data
poi_911_calls = pd.read_csv('data/missoula/poi_911_calls_summary.csv', index_col='b_type')
# Remove 'POI_Count' column for this analysis
call_type_data = poi_911_calls.drop('POI_Count', axis=1)
# Function to get top 10 call types proportions for a given b_type
def get_top_10_call_types_prop(b_type):
b_type_data = call_type_data.loc[b_type]
total_calls = b_type_data['Total']
proportions = (b_type_data.drop('Total') / total_calls * 100).nlargest(10)
return proportions
# Create dropdown menu options
dropdown_options = [{'label': b_type, 'value': b_type} for b_type in call_type_data.index]
# Create the initial plot (we'll use the first b_type as default)
initial_b_type = call_type_data.index[0]
initial_data = get_top_10_call_types_prop(initial_b_type)
fig = make_subplots(rows=1, cols=1)
bar = go.Bar(
x=initial_data.index,
y=initial_data.values,
name=initial_b_type,
text=[f'{val:.2f}%' for val in initial_data.values],
textposition='auto'
)
fig.add_trace(bar)
# Update layout
fig.update_layout(
title=f'Top 10 Call Types for {initial_b_type} (% of Total Calls)',
xaxis_title='Call Type',
yaxis_title='Percentage of Calls',
yaxis=dict(range=[0, 100]), # Set y-axis range from 0 to 100%
updatemenus=[{
'buttons': [
{
'method': 'update',
'label': b_type,
'args': [{'x': [get_top_10_call_types_prop(b_type).index],
'y': [get_top_10_call_types_prop(b_type).values],
'text': [[f'{val:.2f}%' for val in get_top_10_call_types_prop(b_type).values]]},
{'title': f'Top 10 Call Types for {b_type} (% of Total Calls)'}]
} for b_type in call_type_data.index
],
'direction': 'down',
'showactive': True,
}]
)
# Show the plot
fig.show()
# Save the plot as an HTML file for later viewing
fig.write_html("charts/interactive_911_calls_plot_proportional.html")import pandas as pd
import plotly.graph_objects as go
# Load the data
poi_911_calls = pd.read_csv('data/missoula/poi_911_calls_summary.csv', index_col='b_type')
# Remove 'POI_Count' column for this analysis
call_type_data = poi_911_calls.drop('POI_Count', axis=1)
# Function to get top 10 call types for a given b_type
def get_top_10_call_types(b_type):
return call_type_data.loc[b_type].drop('Total').nlargest(10)
# Function to get top 10 call types proportions for a given b_type
def get_top_10_call_types_prop(b_type):
b_type_data = call_type_data.loc[b_type]
total_calls = b_type_data['Total']
proportions = (b_type_data.drop('Total') / total_calls * 100).nlargest(10)
return proportions
# Create the initial plot (we'll use the first b_type as default)
initial_b_type = call_type_data.index[0]
initial_data_count = get_top_10_call_types(initial_b_type)
initial_data_prop = get_top_10_call_types_prop(initial_b_type)
fig = go.Figure()
# Add traces for both count and percentage views
fig.add_trace(go.Bar(
x=initial_data_count.index,
y=initial_data_count.values,
name='Count',
visible=True
))
fig.add_trace(go.Bar(
x=initial_data_prop.index,
y=initial_data_prop.values,
name='Percentage',
text=[f'{val:.2f}%' for val in initial_data_prop.values],
textposition='auto',
visible=False
))
# Update layout
fig.update_layout(
title=f'Top 10 Call Types for {initial_b_type}',
xaxis_title='Call Type',
yaxis_title='Number of Calls',
updatemenus=[
# Dropdown for POI type selection
{
'buttons': [
{
'method': 'update',
'label': b_type,
'args': [
{
'x': [get_top_10_call_types(b_type).index, get_top_10_call_types_prop(b_type).index],
'y': [get_top_10_call_types(b_type).values, get_top_10_call_types_prop(b_type).values],
'text': [None, [f'{val:.2f}%' for val in get_top_10_call_types_prop(b_type).values]]
},
{'title': f'Top 10 Call Types for {b_type}'}
]
} for b_type in call_type_data.index
],
'direction': 'down',
'showactive': True,
'x': 0.1,
'xanchor': 'left',
'y': 1.15,
'yanchor': 'top'
},
# Buttons for toggling between count and percentage views
{
'type': 'buttons',
'direction': 'right',
'x': 0.5,
'xanchor': 'center',
'y': 1.15,
'yanchor': 'top',
'buttons': [
{
'label': "Count",
'method': 'update',
'args': [
{'visible': [True, False]},
{'yaxis': {'title': 'Number of Calls', 'autorange': True}}
]
},
{
'label': "Percentage",
'method': 'update',
'args': [
{'visible': [False, True]},
{'yaxis': {'title': 'Percentage of Calls', 'range': [0, 100]}}
]
}
]
}
]
)
# Show the plot
fig.show()
# Save the plot as an HTML file for later viewing
fig.write_html("charts/interactive_911_calls_plot_combined.html")import pandas as pd
import plotly.graph_objects as go
# Load the data
poi_911_calls = pd.read_csv('data/missoula/poi_911_calls_summary.csv', index_col='b_type')
# Remove 'POI_Count' column for this analysis
call_type_data = poi_911_calls.drop('POI_Count', axis=1)
# Calculate overall top call types
overall_top_calls = call_type_data.drop('Total', axis=1).sum().sort_values(ascending=False)
# Function to get top 10 call types for a given b_type, ignoring top N overall
def get_top_10_call_types(b_type, ignore_top_n):
ignored_calls = set(overall_top_calls.nlargest(ignore_top_n).index)
filtered_data = call_type_data.loc[b_type].drop('Total').drop(ignored_calls, errors='ignore')
return filtered_data.nlargest(10)
# Function to get top 10 call types proportions for a given b_type, ignoring top N overall
def get_top_10_call_types_prop(b_type, ignore_top_n):
ignored_calls = set(overall_top_calls.nlargest(ignore_top_n).index)
b_type_data = call_type_data.loc[b_type]
total_calls = b_type_data['Total']
filtered_data = b_type_data.drop('Total').drop(ignored_calls, errors='ignore')
proportions = (filtered_data / total_calls * 100).nlargest(10)
return proportions
# Create the initial plot
initial_b_type = call_type_data.index[0]
initial_ignore_top_n = 0
initial_data_count = get_top_10_call_types(initial_b_type, initial_ignore_top_n)
initial_data_prop = get_top_10_call_types_prop(initial_b_type, initial_ignore_top_n)
fig = go.Figure()
# Add traces for both count and percentage views
fig.add_trace(go.Bar(
x=initial_data_count.index,
y=initial_data_count.values,
name='Count',
visible=True
))
fig.add_trace(go.Bar(
x=initial_data_prop.index,
y=initial_data_prop.values,
name='Percentage',
text=[f'{val:.2f}%' for val in initial_data_prop.values],
textposition='auto',
visible=False
))
# Function to create slider steps
def create_slider_steps(b_type):
return [
{
'method': 'update',
'label': str(n),
'args': [
{
'x': [get_top_10_call_types(b_type, n).index,
get_top_10_call_types_prop(b_type, n).index],
'y': [get_top_10_call_types(b_type, n).values,
get_top_10_call_types_prop(b_type, n).values],
'text': [None, [f'{val:.2f}%' for val in get_top_10_call_types_prop(b_type, n).values]]
},
{'title': f'Top 10 Call Types for {b_type} (Ignoring top {n} overall)'}
]
} for n in range(0, 11) # Allows ignoring 0 to 10 top overall call types
]
# Update layout
fig.update_layout(
title=f'Top 10 Call Types for {initial_b_type} (Ignoring top {initial_ignore_top_n} overall)',
xaxis_title='Call Type',
yaxis_title='Number of Calls',
updatemenus=[
# Dropdown for POI type selection
{
'buttons': [
{
'method': 'update',
'label': b_type,
'args': [
{
'x': [get_top_10_call_types(b_type, initial_ignore_top_n).index,
get_top_10_call_types_prop(b_type, initial_ignore_top_n).index],
'y': [get_top_10_call_types(b_type, initial_ignore_top_n).values,
get_top_10_call_types_prop(b_type, initial_ignore_top_n).values],
'text': [None, [f'{val:.2f}%' for val in get_top_10_call_types_prop(b_type, initial_ignore_top_n).values]]
},
{
'title': f'Top 10 Call Types for {b_type} (Ignoring top {initial_ignore_top_n} overall)',
'sliders[0].steps': create_slider_steps(b_type)
}
]
} for b_type in call_type_data.index
],
'direction': 'down',
'showactive': True,
'x': 0.5,
'xanchor': 'left',
'y': 1.15,
'yanchor': 'top'
},
# Buttons for toggling between count and percentage views
{
'type': 'buttons',
'direction': 'right',
'x': 0.8,
'xanchor': 'center',
'y': 1.15,
'yanchor': 'top',
'buttons': [
{
'label': "Count",
'method': 'update',
'args': [
{'visible': [True, False]},
{'yaxis': {'title': 'Number of Calls', 'autorange': True}}
]
},
{
'label': "Percentage",
'method': 'update',
'args': [
{'visible': [False, True]},
{'yaxis': {'title': 'Percentage of Calls for this POI', 'range': [0, 100]}}
]
}
]
}
],
# Add a slider to control the number of top overall call types to ignore
sliders=[{
'active': 0,
'currentvalue': {"prefix": "Ignore top N overall: "},
'pad': {"t": 50},
'steps': create_slider_steps(initial_b_type)
}]
)
# Save the plot as an HTML file for later viewing
fig.write_html("charts/interactive_911_calls_plot_with_ignore.html")# List of Violent call types
violent_calls = [
'Assault', 'Assault with a Weapon', 'Domestic Animal', 'Homicide',
'Intimidation', 'Kidnapping', 'Partner Family Member Assault',
'Robbery', 'Sexual Assault', 'Shots Fired', 'Shots Heard',
'Suicidal Person', 'Suicide', 'Weapons - Carry/Possess'
]
# Function to classify call types
def classify_call(call_type):
return 'Violent' if call_type in violent_calls else 'Nonviolent'
# Add new column 'crime_type' to the DataFrame
calls_911_2023['crime_type'] = calls_911_2023['Call Type Description'].apply(classify_call)
# Print the first few rows to verify the new column
print(calls_911_2023[['Call Type Description', 'crime_type']].head())
# Print value counts to see the distribution
print(calls_911_2023['crime_type'].value_counts())
# Optional: If you want to see the distribution as percentages
print(calls_911_2023['crime_type'].value_counts(normalize=True))import pandas as pd
import plotly.graph_objects as go
# Load the data
df = calls_911_2023_2024
df['Call Creation Date and Time'] = df['Call Creation Date and Time'].astype('datetime64[ns]')
# Extract date
df['Date'] = df['Call Creation Date and Time'].dt.date
# Create a dataframe with each day as a row and columns for each call type
daily_counts = df.groupby(['Date', 'Call Type Description']).size().unstack(fill_value=0)
# Add a total column
daily_counts['Total'] = daily_counts.sum(axis=1)
# Sort the index to ensure chronological order
daily_counts = daily_counts.sort_index()
# Function to get top 10 call types for a given date
def get_top_10_call_types(date):
day_data = daily_counts.loc[date].sort_values(ascending=False)
top_10 = day_data.head(10)
return '<br>'.join([f"{call_type}: {count}" for call_type, count in top_10.items() if call_type != 'Total'])
# Create hover text for each day
hover_text = daily_counts.index.map(get_top_10_call_types)
# Create the plot
fig = go.Figure()
# Add trace for total calls (visible by default)
fig.add_trace(
go.Scatter(
x=daily_counts.index,
y=daily_counts['Total'],
name='Total Calls',
line=dict(color='black', width=2),
hovertemplate='<b>Date</b>: %{x}<br>' +
'<b>Total Calls</b>: %{y}<br><br>' +
'<b>Top 10 Call Types:</b><br>%{text}',
text=hover_text
)
)
# Add traces for each call type (hidden by default)
for column in daily_counts.columns:
if column != 'Total':
fig.add_trace(
go.Scatter(
x=daily_counts.index,
y=daily_counts[column],
name=column,
visible=False,
hovertemplate='<b>Date</b>: %{x}<br>' +
f'<b>{column}</b>: %{{y}}<br><br>' +
'<b>Top 10 Call Types:</b><br>%{text}',
text=hover_text
)
)
# Create buttons for filtering
buttons = [dict(
label='Total Calls',
method='update',
args=[{'visible': [True] + [False] * (len(daily_counts.columns) - 1)},
{'title': 'Total 911 Calls by Day'}]
)]
for i, column in enumerate(daily_counts.columns):
if column != 'Total':
visibility = [False] * len(daily_counts.columns)
visibility[i + 1] = True # +1 because Total is the first trace
buttons.append(dict(
label=column,
method='update',
args=[{'visible': visibility},
{'title': f'911 Calls by Day: {column}'}]
))
# Update layout
fig.update_layout(
title='Total 911 Calls by Day',
xaxis_title='Date',
yaxis_title='Number of Calls',
updatemenus=[dict(
active=0,
buttons=buttons,
direction="down",
pad={"r": 10, "t": 10},
showactive=True,
x=0.5,
xanchor="left",
y=1.15,
yanchor="top"
)],
hovermode='closest'
)
# Show the plot
fig.show()
# Save the plot as an HTML file
fig.write_html("charts/911_calls_time_series_plot_by_day_with_popup.html")import pandas as pd
import plotly.graph_objects as go
# Load the data
df = calls_911_2023_2024
df['Call Creation Date and Time'] = df['Call Creation Date and Time'].astype('datetime64[ns]')
df = df[df['Call Creation Date and Time'].dt.strftime('%H-%M-%S') != '00-00-00'] # Remove rows with time 00:00:00 (to avoid issues with plotting a lot of calls at midnight)
# Extract month and day, ignoring year
df['Hour'] = df['Call Creation Date and Time'].dt.hour
# Create a dataframe with each day as a row and columns for each call type
daily_counts = df.groupby(['Hour', 'Call Type Description']).size().unstack(fill_value=0)
# Add a total column
daily_counts['Total'] = daily_counts.sum(axis=1)
# Sort the index to ensure chronological order
daily_counts = daily_counts.sort_index()
# Create the plot
fig = go.Figure()
# Add trace for total calls (visible by default)
fig.add_trace(
go.Scatter(
x=daily_counts.index,
y=daily_counts['Total'],
name='Total Calls',
line=dict(color='black', width=2)
)
)
# Add traces for each call type (hidden by default)
for column in daily_counts.columns:
if column != 'Total':
fig.add_trace(
go.Scatter(
x=daily_counts.index,
y=daily_counts[column],
name=column,
visible=False
)
)
# Create buttons for filtering
buttons = [dict(
label='Total Calls',
method='update',
args=[{'visible': [True] + [False] * (len(daily_counts.columns) - 1)},
{'title': 'Total 911 Calls by Hour'}]
)]
for i, column in enumerate(daily_counts.columns):
if column != 'Total':
visibility = [False] * len(daily_counts.columns)
visibility[i + 1] = True # +1 because Total is the first trace
buttons.append(dict(
label=column,
method='update',
args=[{'visible': visibility},
{'title': f'911 Calls by Hour: {column}'}]
))
# Update layout
fig.update_layout(
title='Total 911 Calls by Hour',
xaxis_title='Hour',
yaxis_title='Number of Calls',
updatemenus=[dict(
active=0,
buttons=buttons,
direction="down",
pad={"r": 10, "t": 10},
showactive=True,
x=0.5,
xanchor="left",
y=1.15,
yanchor="top"
)],
)
# Save the plot as an HTML file
fig.write_html("charts/911_calls_time_series_plot_by_hour.html")import pandas as pd
import folium
from folium.plugins import HeatMap
# Create DataFrame from the provided data
data = calls_911_2023_2024
df = pd.DataFrame(data)
# Create a map centered on Missoula
m = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a feature group for all calls
all_calls_group = folium.FeatureGroup(name="All Calls", show=False)
all_heat_data = df[['Lat', 'Lon']].values.tolist()
HeatMap(all_heat_data).add_to(all_calls_group)
all_calls_group.add_to(m)
# Get unique call types
call_types = df['Call Type Description'].unique()
# Create a feature group for each call type
for call_type in call_types:
filtered_data = df[df['Call Type Description'] == call_type]
heat_data = filtered_data[['Lat', 'Lon']].values.tolist()
feature_group = folium.FeatureGroup(name=call_type, show=False) # Set show=False to make layers off by default
HeatMap(heat_data, radius = 15, min_opacity = 0.4, blur = 15, gradient = {0.4: 'blue', 0.65: 'lime', 1: 'red'}, max_zoom=18).add_to(feature_group)
feature_group.add_to(m)
# Add layer control to the map
folium.LayerControl().add_to(m)
# Save the map
m.save("maps/missoula_911_heatmap_with_all_calls.html")
print("Heatmap with 'All Calls' option has been saved as 'missoula_911_heatmap_with_all_calls.html'")import pandas as pd
import folium
from folium.plugins import HeatMap
# Create DataFrame from the provided data
data = calls_911_2023_2024
df = pd.DataFrame(data)
# Create a map centered on Missoula
m = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a feature group for all calls
all_calls_group = folium.FeatureGroup(name="All Calls", show=False)
all_heat_data = df[['Lat', 'Lon']].values.tolist()
HeatMap(all_heat_data).add_to(all_calls_group)
all_calls_group.add_to(m)
# Get unique crime types
call_types = df['crime_type'].unique()
# Create a feature group for each crime type
for call_type in call_types:
filtered_data = df[df['crime_type'] == call_type]
heat_data = filtered_data[['Lat', 'Lon']].values.tolist()
feature_group = folium.FeatureGroup(name=call_type, show=False) # Set show=False to make layers off by default
HeatMap(heat_data, radius = 15, min_opacity = 0.4, blur = 15, gradient = {0.4: 'blue', 0.65: 'lime', 1: 'red'}, max_zoom=18).add_to(feature_group)
feature_group.add_to(m)
# Add layer control to the map
folium.LayerControl().add_to(m)
# Save the map
m.save("maps/missoula_911_heatmap_with_crime_type.html")
print("Heatmap with 'All Calls' option has been saved as 'missoula_911_heatmap_with_crime_type.html'")import pandas as pd
import folium
from folium.plugins import HeatMap
# Create DataFrame from the provided data
data = calls_911_2023_2024
df = pd.DataFrame(data)
# drop dates with no hms (cant use dt for this)
df = df[df['Call Creation Date and Time'].str.contains(':')]
# Convert the 'Date' column to datetime and extract hour
df['Date'] = pd.to_datetime(df['Call Creation Date and Time'])
df['Hour'] = df['Date'].dt.hour
# Create a map centered on Missoula
m = folium.Map(location=[46.87, -113.99], zoom_start=12)
# Create a feature group for each hour
for hour in range(24):
hour_df = df[df['Hour'] == hour]
heat_data = hour_df[['Lat', 'Lon']].values.tolist()
# Create a name for the layer (e.g., "12 AM" for hour 0, "1 PM" for hour 13)
layer_name = f"{hour % 12 or 12} {'AM' if hour < 12 else 'PM'}"
feature_group = folium.FeatureGroup(name=layer_name, show=False)
HeatMap(heat_data, radius = 15, min_opacity = 0.4, blur = 15, gradient = {0.4: 'blue', 0.65: 'lime', 1: 'red'}, max_zoom=18).add_to(feature_group)
feature_group.add_to(m)
# Add a feature group for all hours
all_heat_data = df[['Lat', 'Lon']].values.tolist()
all_hours_group = folium.FeatureGroup(name="All Hours", show=True)
HeatMap(all_heat_data).add_to(all_hours_group)
all_hours_group.add_to(m)
# Add layer control to the map
folium.LayerControl().add_to(m)
# Save the map
m.save("maps/missoula_911_heatmap_hourly_layers.html")
print("Heatmap with hourly layers has been saved as 'missoula_911_heatmap_hourly_layers.html'")import rasterio
import numpy as np
import pandas as pd
# Open the TIFF file
tif_file = "data/missoula/missoula_pop_density.tif"
with rasterio.open(tif_file) as src:
# Read the raster data
data = src.read(1) # Read the first band (assuming single-band raster)
# Replace NaN (no-data) values with 0
data = np.nan_to_num(data, nan=0) # Replace NaN with 0
# Print out the metadata to verify file details
print("Metadata:", src.meta)
# Check the data shape and some basic statistics
print("Data Shape:", data.shape)
print("Min Value:", np.nanmin(data)) # Should not be NaN anymore
print("Max Value:", np.nanmax(data))
print("Mean Value:", np.nanmean(data))
# Flatten the raster data to create a table (row for each pixel)
data_flat = data.flatten()
# Create a DataFrame from the flattened raster data
df = pd.DataFrame({'Pixel Value': data_flat})
# Show the first few rows of the DataFrame
print(df.head())
#------
import rasterio
import geopandas as gpd
from rasterstats import zonal_stats
import numpy as np
# Load the population density TIFF file
pop_density_file = 'data/missoula/missoula_pop_density.tif'
with rasterio.open(pop_density_file) as src:
population_density = src.read(1) # Read the first band (assuming single-band raster)
transform = src.transform # Affine transformation for pixel coordinates
pixel_size = transform[0] * transform[4] # Get pixel size (e.g., in square kilometers)
# Load the neighborhood boundaries (GeoJSON, Shapefile, etc.)
n_file = gpd.read_file('data/missoula/missoula_neighborhoods_fixed_3.geojson')
# Ensure both datasets are in the same coordinate reference system (CRS)
if n_file.crs != src.crs:
n_file = n_file.to_crs(src.crs)
# Perform zonal statistics to sum population density for each neighborhood
stats = zonal_stats(n_file, pop_density_file, stats="sum", affine=transform)
# Add the population density sums to the neighborhoods GeoDataFrame
n_file['population_density_sum'] = [stat['sum'] for stat in stats]
# Calculate the population for each neighborhood
# Assuming the population density is per square kilometer
# n_file['population'] = n_file['population_density_sum'] * pixel_size
# Display the results
# print(n_file) # Assuming 'name' is the neighborhood name column
print(n_file[['FeatureName', 'population_density_sum']])
n_file.to_file("pop_dens_sum.geojson", driver="GeoJSON")
#sums pops
i_sum = 0
for i in n_file["population_density_sum"]:
i_sum += i
print("totaled:", i_sum)
#save into file
#csv:
# n_file[['FeatureName', 'population_density_sum']].to_csv('neighbor_pop_dens.csv', index=False)
#------
n_file
#------
#1
#----------------------
import folium
import branca.colormap as cm
#merge at geometry
merged = n_file.merge(n_file.sjoin(calls_911_2023_2024).groupby('FeatureName').size().rename('CallAmount').reset_index())
merged
#----------------------
#2
#----------------------
merged.to_file("n_p_sum.geojson", driver="GeoJSON")
merged[["FeatureName", "CallAmount"]]
#----------------------
#3
#----------------------
m = folium.Map(location=[46.8721, -113.9940], zoom_start=12)
#colormap
colormap = cm.linear.YlOrRd_09.scale(merged["CallAmount"].min(), merged["CallAmount"].max())
colormap.caption = "Call Amount by Neighborhood"
#color neighs
def style_function(feature):
call_amount = feature['properties']['CallAmount']
return {
'fillOpacity': 0.7,
'weight': 0.5,
'fillColor': colormap(call_amount) if call_amount is not None else '#ffffff', # White if no data
'color': 'black'
}
folium.GeoJson(
merged,
style_function=style_function,
tooltip=folium.GeoJsonTooltip(fields=['FeatureName', 'CallAmount'], aliases=['Neighborhood:', 'Calls:']),
).add_to(m)
colormap.add_to(m)
#save map
m.save('maps/neighborhood_calls_map.html')
#----------------------# load zillow data
missoula_zillow = pd.read_csv('data/missoula/missoula_zillow_10_2024.csv')
# drop values where price is less than 3 characters
missoula_zillow = missoula_zillow[missoula_zillow['price'].str.len() > 3]
missoula_zillow.reset_index(drop=True, inplace=True)missoula_zillow# drop values where price is less than 3 characters
missoula_zillow = missoula_zillow[missoula_zillow['price'].str.len() > 3]
missoula_zillow.reset_index(drop=True, inplace=True)missoula_zillow.iloc[3]# show the first entry image using the image_url
from IPython.display import Image
Image(url=missoula_zillow.iloc[3].image_url)missoula_zillow.iloc[3].price# missoula_zillow.address
# write address to a file and open it in a text editor
missoula_zillow.address.to_csv('missoula_zillow_addresses.txt', index=False)# open missoula_zillow.csv and merge it with the missoula_zillow dataframe (original_address from missoula_zillow.csv and address from missoula_zillow df)
missoula_zillow_original = pd.read_csv('data/missoula/missoula_zillow.csv')
missoula_zillow_original
# merge the two dataframes
missoula_zillow_merged = pd.merge(missoula_zillow, missoula_zillow_original, left_on='address', right_on='original_address')
missoula_zillow_merged.to_csv('missoula_zillow_merged.csv', index=False)# format the price column to be just the number
missoula_zillow_merged['price'] = missoula_zillow_merged['price'].str.replace('$', '').str.replace(',', '').str.replace('+', '').astype(int)
missoula_zillow_merged['price']missoula_zillow_merged.to_csv('missoula_zillow_merged.csv', index=False)# read merged
missoula_zillow_merged = pd.read_csv('missoula_zillow_merged.csv')missoula_zillow_merged.price.plot()