<jupyter_start><jupyter_text>Notation:SAL- small areaPP- police precinct<jupyter_code>from random import shuffle, randint
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Polygon, Point, MultiPoint, MultiPolygon, LineString, mapping, shape
from descartes import PolygonPatch
import random
import fiona
import numpy as np
import csv
from fiona import collection

import geopandas as gpd
from geopandas.tools import sjoin # rtree index in-build, used with inner, intersection
import pandas as pd

from collections import defaultdict

def csv_dict_reader(file_obj):

    """

    Read a CSV file using csv.DictReader

    """

    reader = csv.DictReader(file_obj, delimiter=',')
    data_list=[]
    for line in reader:
        data_list.append(line)
        
    return data_list<jupyter_output><empty_output><jupyter_text>Adding crime numbers to the shapefiles<jupyter_code>with open("data/pp_murder_2015.csv") as f_obj:
        pp_data = csv_dict_reader(f_obj)
        
# rearrange the dict so the districts are directly the keys
pp_data_d = {}
pp_data_prov = {}

for i in pp_data:
    pp_data_d[i['PoliceStation'].lower()] = int(i['Incidents']) 
    pp_data_prov[i['PoliceStation'].lower()] = i['Province'] 

# police precincts
polSHP  = 'shapefiles/oryginal/pol_prec.shp'


# aggregate values from array x to a newly created key in d
def aggregate(data, key1, xdata, key2): 
    i = 1 # check of the num of crimes is present
    try:
            temp = xdata[data[key2]]
    except:
            temp = 0
            
    data[key1] = temp 
    return i

with collection(polSHP, 'r') as input:
    sink_schema = input.schema.copy()
    sink_schema['properties']['murd_cnt'] = 'int'
    sink_schema['properties']['province'] = 'str:25'

    # not specifying the .shp creates a dir with all files in it
    with collection("shapefiles/updated/polPrec_murd2015_prov", "w", "ESRI Shapefile", sink_schema) as sink:
        
        for pt in input:
                pt_tosave = pt
                # will be adding a property
                pt_prop = pt['properties']       
                pt_prop['compnt_nm'] = pt_prop['compnt_nm'].lower() # change to lower case to be able to compare
        
                check_exists = aggregate(pt_prop, 'murd_cnt', pp_data_d, 'compnt_nm') # modifies d on the spot
                check_exists = aggregate(pt_prop, 'province', pp_data_prov, 'compnt_nm') # modifies d on the spot

                pt_tosave['properties'].update(pt_prop)
                
                sink.write(pt_tosave)

geo_sal.head()
polSHP  = 'shapefiles/oryginal/pol_prec.shp'

polSHP.format_map<jupyter_output><empty_output><jupyter_text>Adding population data to the SAL shapefiles<jupyter_code>salSHP = 'shapefiles/oryginal/SAL_SA_2013.shp'

with open("data/sal_processed_head.csv") as f_obj:

      sal_data = csv_dict_reader(f_obj)

sal_data_d = {}

# put it in a dictionary
for i in sal_data:
    sal_data_d[int(i["district"])]=int(i["count"]) 



# aggregate values from array x to a newly created key in d- it enforces int on the key
def aggregate_int(d, key1, x, key2): 
    i = 1 # check of the num of crimes is present
    try:
            temp = x[int(d[key2])]
    except:
            temp = 0
            
    d[key1] = temp #x[d[key2]]
    return i

with collection(salSHP, 'r') as input:
    sink_schema_sal = input.schema.copy()
    sink_schema_sal['properties']['PPL_CNT'] = 'int'
    
    # not specifying the .shp creates a dir with all files in it
    with collection("shapefiles/updated/sal_population", "w", "ESRI Shapefile", sink_schema_sal) as sink:
        
        for pt in input:
            d = pt['properties']
            # this bit is for city of Cape Town
            #if (d['MN_NAME'] == 'City of Cape Town'): # we restrict to the City of CT
            f = pt       
        
            check_exists = aggregate_int(d, 'PPL_CNT', sal_data_d, 'SAL_CODE') # modifies d on the spot
                
                #sal_shp_all_data[i]['PPL_CNT'] = sal_data_d[int(sal_shp_all_data[i]['SAL_CODE'])]

            f['properties'].update(d)
                
            sink.write(f)<jupyter_output><empty_output><jupyter_text>checking which column entries to keep:<jupyter_code>#geo_sal = gpd.GeoDataFrame.from_file(salSHP)
geo_sal.filter(regex=("CODE")).head(n=2)
import geopandas as gpd
from geopandas.tools import sjoin # rtree index in-build, used with inner, intersection

def find_intersections(o):
    
    paired_ind = [o.polPrecincts_index, o.sal_index]

    d_over_ind = defaultdict(list)

    # creating a dictionary that has prescints as keys and associated small areas as values
    for i in range(len(paired_ind[0].values)):
        if not paired_ind[0].values[i]==paired_ind[1].values[i]: # it shows itself as intersection
            d_over_ind[paired_ind[0].values[i]].append(paired_ind[1].values[i])

    # get rid of the pol precincts with no small areas associated to them- not the most efficient way
    d_temp = {}
    for l in d_over_ind:
        if len(d_over_ind[l]):
            d_temp[l] = d_over_ind[l]

    return d_temp
    
    
def calculate_join_indices(g1, g2):

        # A: region of the police data with criminal record
        # C: small area with population data
        # we look for all small areas intersecting a given C_i, calculate the fraction of inclusion, scale the
        # population accordingly: area(A_j, where A_j crosses C_i)/area(A_j)* popul(A_j)

        g1_reind = g1.reset_index().rename(columns={'index':'polPrecincts_index'})
        g2_reind = g2.reset_index().rename(columns={'index':'sal_index'})

        # the actual indexing:
        out = sjoin(g1_reind, g2_reind, how ="inner", op = "intersects")
        out.drop('index_right', axis=1, inplace=True) # there is a double index fo smal areas, so we drop one
        #out_sorted = out.sort(columns='polPrecincts_index', ascending=True) # guess sorting is not necessary, cause we are
        # using doctionaries at later stages
        #dict_over_ind = find_intersections(out_sorted)

        # output retains only 1 area (left or right join), and gives no intersection area.
        # so we create an array with paired indices: police precincts with associated small areas
        # we use it in a loop in a function below
        dict_over_ind = find_intersections(out) 
        
        return dict_over_ind, g1_reind, g2_reind
    
def calculate_join(dict_over_ind, g1_reind, g2_reind, out_name):
        area_total = 0
        data_aggreg = []

        # note to self: make sure to import shapely Polygon
        for index1, crim in g1_reind.iterrows():
            try:
                index1 = crim.polPrecincts_index
                sals_found = dict_over_ind[index1]

                for sal in range(len(sals_found)):
                    pom = g2_reind[g2_reind.sal_index == sals_found[sal]]['geometry']        

                    #if pom.intersects(crim['geometry']).values[0]:
                    area_int = pom.intersection(crim['geometry']).area.values[0]
                    if area_int>0:
                        area_total += area_int 
                        area_crim = crim['geometry'].area

                        area_popu = pom.values[0].area

                        popu_count = g2_reind[g2_reind.sal_index == sals_found[sal]]['PPL_CNT'].values[0]
                        murd_count = crim['murd_cnt']
                        pol_province = crim['province']
                        popu_frac = (area_int / area_popu) * popu_count# fraction of the pop area contained inside the crim
                        #print(popu_frac)
                        extra_info_col_names = ['DC_NAME','MN_NAME','MP_NAME','PR_NAME','SP_NAME']
                        
                        extra_info_col_codes = ['MN_CODE','MP_CODE','PR_CODE','SAL_CODE','SP_CODE']

                        extra_names = g2_reind[g2_reind.sal_index == sals_found[sal]][extra_info_col_names]#.filter(regex=("NAME"))
                        extra_codes = g2_reind[g2_reind.sal_index == sals_found[sal]][extra_info_col_codes]#.filter(regex=("NAME"))

                        data_aggreg.append({'geometry': pom.intersection(crim['geometry']).values[0], 'id1': index1,\
                                      'id2': sals_found[sal] ,'area_crim': area_crim,'area_pop': area_popu,\
                                  'area_inter': area_int, 'popu_frac' : popu_frac, 'popu_count': popu_count,\
                                  'murd_cnt': murd_count,'province': pol_province,
                                  'DC_NAME': extra_names.DC_NAME.values[0],\
                                  'MN_NAME': extra_names.MN_NAME.values[0], 'MP_NAME': extra_names.MP_NAME.values[0],\
                                  'PR_NAME': extra_names.PR_NAME.values[0],'SP_NAME': extra_names.SP_NAME.values[0],\
                                  'MN_CODE': extra_codes.MN_CODE.values[0],'MP_CODE': extra_codes.MP_CODE.values[0],\
                                  'PR_CODE': extra_codes.PR_CODE.values[0],'SAL_CODE': extra_codes.SAL_CODE.values[0],\
                                  'SP_CODE': extra_codes.SP_CODE.values[0]} )
            except:
                pass
            
        df_t = gpd.GeoDataFrame(data_aggreg,columns=['geometry', 'id1','id2','area_crim',\
                                       'area_pop','area_inter', 'popu_frac',\
                                       'popu_count', 'murd_cnt','province','DC_NAME',\
                                       'MN_NAME','MP_NAME','PR_NAME','SP_NAME',\
                                      'MN_CODE','MP_CODE','PR_CODE','SAL_CODE','SP_CODE'])
        df_t.to_file(out_name)
        return df_t, area_total, data_aggreg<jupyter_output><empty_output><jupyter_text>Main functions to find intersection:<jupyter_code>salSHP_upd = 'shapefiles/updated/sal_population/sal_population.shp'
polSHP_upd = 'shapefiles/updated/polPrec_murd2015_prov/polPrec_murd2015_prov.shp'

geo_pol = gpd.GeoDataFrame.from_file(polSHP_upd)
geo_sal = gpd.GeoDataFrame.from_file(salSHP_upd)



#dict_int, geo_pol_reind, geo_sal_reind = calculate_join_indices(geo_pol,geo_sal)
#out_name ='data/intersection_int_ext.shp'
#df, sum_area, data = calculate_join(dict_int, geo_pol_reind, geo_sal_reind, out_name)

geo_sal.geometry.head()<jupyter_output><empty_output><jupyter_text>find pol precincts within WC boundary<jupyter_code>za_province = gpd.read_file('za-provinces.topojson',driver='GeoJSON')#.set_index('id')
za_province.crs={'init': '27700'}


wc_boundary = za_province.ix[8].geometry # WC
#pp_WC = geo_pol[geo_pol.geometry.within(wc_boundary)]
pp_WC_in = geo_pol[geo_pol.geometry.intersects(wc_boundary)]
#.unary_union
pp_WC_overlaps = pp_WC_in[pp_WC_in.province!="Western Cape"]
pp_WC_pol_annot = pp_WC_in[pp_WC_in.province=="Western Cape"]
#pp_test = pp_WC_in[pp_WC_in['compnt_nm'].isin(['atlantis','philadelphia','kraaifontein','brackenfell','kuilsriver','kleinvleveerste river','macassar','somerset west','fish hoek'])]
#pp_test = pp_WC_in[pp_WC_in['compnt_nm'].isin(['beaufort west','doring bay','murraysburg', 'strandfontein','nuwerus','lutzville'])]
%matplotlib inline
pp_WC_overlaps.plot()
sal_wc_union_bound = sal_WC_in.unary_union
df_int =gpd.GeoDataFrame.from_file('data/intersection_int_ext.shp')

df_int.head(n=2)
# add population data per police percinct to the main table
temp = df_int.groupby(by=['id1'])['popu_frac'].sum().reset_index()

data_with_pop = pd.merge(df_int, temp, on='id1', how='outer')\
        .rename(columns={'popu_frac_y':'popu_frac_per_pp', 'popu_frac_x':'popu_inter'})

# finally, update the murder rate per SAL : id2 is sal's id    

data_with_pop['murd_est_per_int'] = data_with_pop['popu_inter']/data_with_pop['popu_frac_per_pp']\
           * data_with_pop['murd_cnt']
data_mur_per_int = data_with_pop.groupby(by=['id2'])['murd_est_per_int'].sum().reset_index()

data_mur_per_sal = data_mur_per_int.rename(columns={'murd_est_per_int':'murd_est_per_sal'})

data_with_pop['ratio_per_int'] = data_with_pop['popu_inter']/data_with_pop['popu_frac_per_pp']\
           
data_full = pd.merge(data_with_pop, data_mur_per_sal, on='id2', how='outer')\
        .rename(columns={'id1':'index_PP', 'id2':'index_SAL'})
    
# the geometry can be simplified using Line Simpl. algorithms
# data_with_pop['geo_simplified'] = data_with_pop.geometry.simplify(1000)
# for use one switches to a chosen simplified geometry:
# data_with_pop.set_geometry('geo_simplified', inplace=True)

    
# sth if off when saving data_with_pop directly- complaines about being a Series object. 
#test = gpd.GeoDataFrame(data_with_pop)
#test['geo_simplified'] = test.geometry.simplify(1000)
#test.to_file('crime/GIS-analyses/pooled_data/intersection.shp')

data_full.head(n=2)<jupyter_output><empty_output><jupyter_text>Small areas and police districts might have different annotation<jupyter_code>data_prov = data_full[['PR_NAME','province','murd_est_per_int']]<jupyter_output><empty_output><jupyter_text>Calculation of murders per province wrt police annotation of provinces:<jupyter_code>data_prov.groupby('province')['murd_est_per_int'].sum()<jupyter_output><empty_output><jupyter_text>same as above, but wrt to SAL annotation:<jupyter_code>data_prov.groupby('PR_NAME')['murd_est_per_int'].sum()
# check over small areas
pom = {}
for ind, row in data_full.iterrows():
    pom[row['index_SAL']] = row['murd_est_per_sal'] 
s=0
for key in pom:
    s = s + pom[key]
print(s)
#f, ax = plt.subplots(1, figsize=(12, 12))
#wc.plot(ax=ax, color='grey', linewidth=0)
#gpd.plotting.plot_multipolygon(ax, sal_wc_union_bound, facecolor='green')

#gpd.plotting.plot_multipolygon(ax, sal_wc_union, facecolor='blue')
f, ax = plt.subplots(1, figsize=(12, 12))
gpd.plotting.plot_multipolygon(ax, pp_WC_overlaps.unary_union, facecolor='red')
wc.plot(ax=ax, color='black', linewidth=0) 
data_full[data_full.MN_NAME=="City of Cape Town"]['murd_est_per_int'].sum()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11641791044776119, 'fraction_numerical': 0.008884150675195452, 'mean_word_length': 3.233152827918171, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '514645', 'license': 'MIT', 'n_tokens_mistral': 4986, 'n_tokens_neox': 4739, 'n_words': 961, 'path': 'crime_stats_compute.ipynb', 'repo_name': 'OpenUpSA/crime-stats-demystifed', 'size': '14070'}