Display VAW

# Import required modules
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx
import plotly.express as px
import folium
from folium import plugins

# Read the data
gdf = gpd.read_file('data/ageb/agbumge19gw.shp')
gdf.sample(2)

gdf = gdf[gdf['CVE_ENT'] == '09']

gdf = gdf[['CVEGEO', 'geometry']]

fig, ax = plt.subplots(figsize=(12,12))

gdf.plot(ax=ax)

plt.axis('off') 
plt.show()

gdf.explore()

# check the CRS
gdf.crs

# reproject to web mercator
gdf = gdf.to_crs(epsg=3857)

fig, ax = plt.subplots(figsize=(12,12))

gdf.plot(ax=ax,
         color='black', 
         edgecolor='black',
         lw=0.5,
         alpha=0.2,
         )

ax.axis('off')

# add basemap
ctx.add_basemap(ax,source=ctx.providers.OpenStreetMap.Mapnik)

df = pd.read_csv('data/census_data2020.csv')
df.head()

gdf = pd.merge(gdf, df[['CVEGEO', 'POBTOT', 'PRESOE15',
                        'P5_HLI', 'POB_AFRO', 'GRAPROES',
                        'PDESOCUP', 'PSINDER', 'VPH_SINLTC' ,
                        'VPH_SINCINT']], on="CVEGEO")

gdf = gdf[gdf['POBTOT']>100]

# Read data
vaw = pd.read_csv('data/vaw_filtered.csv')

# create a copy to compute required values for the plot
data = vaw.copy() 

# convert time to a datetime64 type 
data['Time'] = pd.to_datetime(data['Time'])
data['Month'] = data['Time'].dt.month # extract the numeric value of months

# group crime by type and count occurrence
occur_type = pd.DataFrame(data.groupby(['Crime', 'Month']).size())
occur_type = occur_type.rename(columns = {0:'Reported'}) # rename the column

# order from Jan to Dec
occur_type.sort_values(["Month"],axis=0, ascending=True, inplace=True)
occur_type.reset_index(inplace = True) # reset index

occur_type.tail(5)

# and plot it
fig = px.bar(occur_type, x="Month", y='Reported',
             color="Crime",
             hover_data={'Crime' : True ,
                         'Month' : False
                               },
             # costumize colours
             color_discrete_map={"DOMESTIC VIOLENCE": "darkblue",
                                 "FEMINICIDE": "royalblue",
                                 "SEXUAL HARASSMENT" : "cornflowerblue",
                                 "SEXUAL ABUSE": "deeppink",
                                },
             labels={"Crime" : "Type of violence",
                     "Reported": "Reported violence"},
             title='Reported Violence Against Women in CDMX by Type, 2021')


fig.update_layout(
        legend=dict(
        orientation="h",
        y = - 0.2
        ),
    legend_title_text=''

)


config = dict({'displayModeBar': False})

#fig.write_html("output/VAW_in_CDMX_byType.html",config=config)
#fig.show(config=config)

# get the data types
vaw.info()

# convert df to geodataframe
vaw_gdf = gpd.GeoDataFrame(vaw,
                         crs='EPSG:4326', # lat and long are in degree
                         geometry= gpd.points_from_xy(vaw.longitude, vaw.latitude))
vaw_gdf.head(3)

# plot hexbins with Matplotlib
fig, ax = plt.subplots(figsize=(12, 12))

hb= ax.hexbin(data=vaw_gdf,
    x = 'longitude',
    y= 'latitude',
    gridsize=50, # size of hexagonal cells
    alpha=0.5, 
    cmap='viridis_r'
)

# no axis
ax.axis('off')
# add bar side
fig.colorbar(hb, aspect=50, pad=0 );

# hexbin takes lat and lon as arguments --> specify crs
ctx.add_basemap(ax,
                source=ctx.providers.OpenStreetMap.Mapnik,
                crs=4326)

# reproject vaw_gdf to web mercator
vaw_gdf = vaw_gdf.to_crs(epsg=3857)

join = gpd.sjoin(vaw_gdf, gdf, how='left')
join.head(3)

reported_by_gdf = join.CVEGEO.value_counts().rename_axis('CVEGEO').reset_index(name='VAW_COUNT')
reported_by_gdf.head()

gdf= gdf.merge(reported_by_gdf,on='CVEGEO')
gdf.head()

gdf['VAW_per_1000'] = gdf['VAW_COUNT']/gdf['POBTOT']*1000

# Write to GeoJSON to avoid recalculation 

gdf.to_file("data/VAW_AGEB_clean.gjson", driver="GeoJSON")

# plot the top 10 geographies
fig,ax = plt.subplots(figsize=(15,15))

gdf.sort_values(by='VAW_per_1000',ascending=False)[:10].plot(ax=ax,
                                                                 color='red',
                                                                 edgecolor='white',
                                                                 alpha=0.5,legend=True)


# title
ax.set_title('Top 10 Most Dangerous Blocks for Women in CDMX, 2021', fontsize=22)

# no axis
ax.axis('off')

# add a basemap
ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)

fig,ax = plt.subplots(figsize=(15,15))

gdf.plot(ax=ax,
        column='VAW_per_1000',
        legend=True,
        alpha=0.6,
        cmap='RdYlGn_r',
        scheme='quantiles')

ax.axis('off')
ax.set_title('Rates of violence against women per 1000 residents in CDMX, 2021',fontsize=22)

ctx.add_basemap(ax,source=ctx.providers.CartoDB.Positron)

# Read the data
gdf_mun = gpd.read_file('data/CDMX_municipalities.gjson')
# reproject to web mercator
gdf_mun = gdf_mun.to_crs(epsg=3857)

gdf_mun.head(3)

# Make a spatial join
join = gpd.sjoin(vaw_gdf, gdf_mun, how="left")
# count vaw per neighbourhood
reported_by_gdf = join.CVEGEO.value_counts().rename_axis('CVEGEO').reset_index(name='VAW_COUNT')

gdf_mun= gdf_mun.merge(reported_by_gdf,on='CVEGEO')

#Compute rate per 1000
gdf_mun['VAW_per_1000'] = gdf_mun['VAW_COUNT']/gdf['POBTOT']*1000
gdf_mun.head()

#Create map instance, set location
vaw_map =folium.Map(location=[19.32096,-99.15261], 
                   tiles = 'cartodbpositron', 
                   zoom_start=10, control_scale=True, overlay=False)
#set two layer
fg1 = folium.FeatureGroup(name='VAW rate per AGEB',overlay=False).add_to(vaw_map)
fg2 = folium.FeatureGroup(name='VAW rate per Neighbourhood',overlay=False).add_to(vaw_map)

#Add the first choropleth map layer to fg1
custom_scale1 = (gdf['VAW_per_1000'].quantile((0,0.2,0.4,0.6,0.8,1))).tolist()

ageb = folium.Choropleth(
        geo_data=gdf,
        name='vaw in 2021',
        data=gdf,
        columns=['CVEGEO', 'VAW_per_1000'],
        key_on='feature.properties.CVEGEO',
        fill_color='YlOrRd',
        fill_opacity=0.7,
        line_opacity=0.2,
        line_color='white', 
        line_weight=0, 
        smooth_factor=1,
        threshold_scale=custom_scale1,
        legend_name= 'Rates of violence against women per 1000 in CDMX, 2021',
    overlay=False).geojson.add_to(fg1)


folium.features.GeoJson(
                    data=gdf,
                    smooth_factor=2,
                    style_function=lambda x: {'color':'white','fillColor':'transparent','weight':0.5},
                    tooltip=folium.features.GeoJsonTooltip(
                        fields=['POBTOT',
                                'VAW_COUNT',
                                'VAW_per_1000',
                                'PDESOCUP',
                                'GRAPROES'
                               ],
                        aliases=["Total Population, <br>(Per AGEB):",
                                 "Reported cases:",
                                 "VAW rate per 1000:",
                                 "Pop Unemployed:",
                                 "Degree of Education:"
                                ], 
                        localize=True,
                        sticky=False,
                        labels=True,
                        style="""
                            background-color: #F0EFEF;
                            border: 2px solid black;
                            border-radius: 3px;
                            box-shadow: 3px;
                        """,
                        max_width=800,),
                            highlight_function=lambda x: {'weight':3,'fillColor':'grey'},
                        ).add_to(ageb)

#Add the 2nd choropleth map layer to fg2
custom_scale2 = (gdf_mun['VAW_per_1000'].quantile((0,0.2,0.4,0.6,0.8,1))).tolist()

neighb = folium.Choropleth(
    geo_data=gdf_mun,
    name='vaw per Neighbourhood',
    data=gdf_mun,
    columns=['CVEGEO', 'VAW_per_1000'],
    key_on='feature.properties.CVEGEO',
    fill_color='YlOrRd',
    fill_opacity=0.7,
    line_opacity=0.2,
    line_color='white', 
    line_weight=0, 
    smooth_factor=1,
    threshold_scale=custom_scale2,
    legend_name= 'Incidents of violence against women in CDMX, 2021',
    overlay=False).geojson.add_to(fg2)

folium.features.GeoJson(
                    data=gdf_mun,
                    smooth_factor=2,
                    style_function=lambda x: {'color':'white','fillColor':'transparent','weight':0.5},
                    tooltip=folium.features.GeoJsonTooltip(
                        fields=['NOMGEO',
                                'POBTOT',
                                'VAW_COUNT',
                                'VAW_per_1000',
                               ],
                        aliases=["Neighbourhood:",
                                 "Population:",
                                 "Reported Cases:",
                                 "VAW rate per 1000:"
                                 
                                ], 
                        localize=True,
                        sticky=False,
                        labels=True,
                        style="""
                            background-color: #F0EFEF;
                            border: 2px solid black;
                            border-radius: 3px;
                            box-shadow: 3px;
                        """,
                        max_width=800,),
                            highlight_function=lambda x: {'weight':3,'fillColor':'grey'},
                        ).add_to(neighb)   


#Add layer control to the map
folium.TileLayer('cartodbdark_matter',overlay=True,name="View in Dark Mode").add_to(vaw_map)
folium.TileLayer('cartodbpositron',overlay=True,name="Viw in Light Mode").add_to(vaw_map)
folium.LayerControl(collapsed=False).add_to(vaw_map)
#Add mini map
folium.plugins.MiniMap().add_to(vaw_map)

#vaw_map.save("output/CDMX_VAW_Folium.html") #save to a file
#vaw_map