In [2]:
# --- Imports ---
import sys,os
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt 
import seaborn as sns
%matplotlib inline

SF drug geography

Data pulling and cleaning -

Let's pull all SF Crime data provided by SF data:

Let's pull it in and peek at the schema.

In [3]:
d_crime=pd.read_csv('/Users/lmartin/Desktop/Misc/Blog_Archive/GIS/Map__Crime_Incidents_-_from_1_Jan_2003.csv')
print d_crime.shape
d_crime.head(1)
(1723961, 12)
Out[3]:
IncidntNum Category Descript DayOfWeek Date Time PdDistrict Resolution Address X Y Location
0 50436712 ASSAULT BATTERY Wednesday 04/20/2005 12:00:00 AM 04:00 MISSION NONE 18TH ST / CASTRO ST -122.435003 37.760888 (37.7608878061245, -122.435002864271)

It provides ~ 1M records with:

Here's a nice map of the districts: http://sf-police.org/index.aspx?page=796

Let's create an easy handle (days) for timeseries analysis.

In [4]:
date=pd.to_datetime(d_crime['Date'])
print date.min()
print date.max()
t_delta=(date-date.min()).astype('timedelta64[D]')
d_crime['days']=t_delta
d_crime.head(1)
2003-01-01 00:00:00
2015-02-13 00:00:00
Out[4]:
IncidntNum Category Descript DayOfWeek Date Time PdDistrict Resolution Address X Y Location days
0 50436712 ASSAULT BATTERY Wednesday 04/20/2005 12:00:00 AM 04:00 MISSION NONE 18TH ST / CASTRO ST -122.435003 37.760888 (37.7608878061245, -122.435002864271) 840

The first recorded request is 2003-01-01 and most recent is 2015-02-13. Nice.

Distributions and exploration -

In [411]:
def plotdat(data,cat):
    l=data.groupby(cat).size()
    l.sort()
    fig=plt.figure(figsize=(10,5))
    plt.yticks(fontsize=8)
    l.plot(kind='bar',fontsize=12,color='k')
    plt.xlabel('')
    plt.ylabel('Number of reports',fontsize=10)

plotdat(d_crime,'Category')

Let's get a more detailed view by examining Descript, which is the particular crime type.

In [405]:
l=d_crime.groupby('Descript').size()
l.sort()
print l.shape
(912,)

Since there's 912 different crime types, let's slice by percentile and peek at the top types of crime for each PdDistrict.

In [409]:
def types_districts(d_crime,per):
    
    # Group by crime type and district 
    hoods_per_type=d_crime.groupby('Descript').PdDistrict.value_counts(sort=True)
    t=hoods_per_type.unstack().fillna(0)
    
    # Sort by hood sum
    hood_sum=t.sum(axis=0)
    hood_sum.sort(ascending=False)
    t=t[hood_sum.index]
    
    # Filter by crime per district
    crime_sum=t.sum(axis=1)
    crime_sum.sort()
    
    # Large number, so let's slice the data.
    p=np.percentile(crime_sum,per)
    ix=crime_sum[crime_sum>p]
    t=t.loc[ix.index]
    return t
    
t=types_districts(d_crime,96)

Cluster the non-normalized data across the top percentile reports and each PdDistrict.

In [9]:
sns.clustermap(t) 
Out[9]:
<seaborn.matrix.ClusterGrid at 0x10299dd10>

Normalize verically across PdDistrict.

In [18]:
sns.clustermap(t,standard_scale=1) 
Out[18]:
<seaborn.matrix.ClusterGrid at 0x10a782c10>

Normalize horizontally across crime types.

In [10]:
sns.clustermap(t,standard_scale=0) 
Out[10]:
<seaborn.matrix.ClusterGrid at 0x1096d1110>

(1) GTA is the most common crime in most PdDistricts.

  • Tenderloin is an outlier, enriched in base/rock crack and narcotics.

(2) For the distribution of crime across areas:

  • Southern: Theft, including theft from auto.
  • Tenderloin: Base / rock crack and narcotics.
  • Bayview: Violence and threats.

Now, let's drill down on a specific question -

Lets, re-examine the crime types.

In [242]:
plotdat(d_crime,'Category')

I'm interested in DRUG/NARCOTIC:

  • I think it will show some interesting dynamics.
  • I think different areas of the city will have different distributions.
In [80]:
# Let's drill down onto one
cat=d_crime[d_crime['Category']=='DRUG/NARCOTIC']
c=cat['Descript'].value_counts()
c.sort(ascending=False)
c.head(10)
Out[80]:
POSSESSION OF NARCOTICS PARAPHERNALIA       19795
POSSESSION OF BASE/ROCK COCAINE             14003
POSSESSION OF MARIJUANA                     10815
SALE OF BASE/ROCK COCAINE                    8691
POSSESSION OF BASE/ROCK COCAINE FOR SALE     7204
POSSESSION OF METH-AMPHETAMINE               7112
POSSESSION OF MARIJUANA FOR SALES            5456
POSSESSION OF CONTROLLED SUBSTANCE           4123
POSSESSION OF HEROIN                         3954
POSSESSION OF COCAINE                        2897
dtype: int64

We can use what we had above, but we simply slice the input data on a category first (above).

In [81]:
t=types_districts(cat,70)
In [82]:
sns.clustermap(t)
Out[82]:
<seaborn.matrix.ClusterGrid at 0x141545450>
In [83]:
sns.clustermap(t,standard_scale=1)
Out[83]:
<seaborn.matrix.ClusterGrid at 0x151752f90>
In [84]:
sns.clustermap(t,standard_scale=0)
Out[84]:
<seaborn.matrix.ClusterGrid at 0x1127195d0>

Nice. We could study these for a while.

But, here's the point:

I think we can simplify this if we compress different types of drug into groups.

Then, we can examine both temporal and spatail profiles.

Drug dynamics -

We'll create a 30 day window.

In [376]:
# Let's drill down onto one
cat=d_crime[d_crime['Category']=='DRUG/NARCOTIC']

# Bin crime by 30 day window
cat['Month']=np.floor(cat['days']/30) # Approximate month (30 day window)

# Default
district='All'
In [377]:
def timeseries(dat,per):
    ''' Category grouped by month '''
    
    # Group by crime type and district 
    cat_per_time=dat.groupby('Month').Descript.value_counts(sort=True)
    t=cat_per_time.unstack().fillna(0)
        
    # Filter by crime per district
    crime_sum=t.sum(axis=0)
    crime_sum.sort()
    
    # Large number, so let's slice the data.
    p=np.percentile(crime_sum,per)
    ix=crime_sum[crime_sum>p]
    t=t[ix.index]
    return t
    
t_all=timeseries(cat,0)

Let's group the drug categories to make this easier to examine.

In [190]:
barituate_features=['SALE OF BARBITUATES',
                    'POSSESSION OF BARBITUATES FOR SALES',
                    'ENCOURAGE MINOR TO USE BARBITUATES',
                    'POSSESSION OF BARBITUATES'] 
In [212]:
coke_features=['ENCOURAGING MINOR TO USE COCAINE',
               'SALES COCAINE BASE/SCHOOLYARD TRAFFICKING ACT VIO',
               'TRANSPORTATION OF COCAINE',
               'SALE OF COCAINE',
               'POSSESSION OF COCAINE FOR SALES',
               'POSSESSION OF COCAINE']
In [210]:
weed_features=['ENCOURAGING MINOR TO USE MARIJUANA',
               'FURNISHING MARIJUANA',
               'PLANTING/CULTIVATING MARIJUANA',
               'TRANSPORTATION OF MARIJUANA',
               'SALE OF MARIJUANA',
               'POSSESSION OF MARIJUANA FOR SALES',
               'POSSESSION OF MARIJUANA']
In [171]:
metadone_features=['TRANSPORTATION OF METHADONE',
                   'SALE OF METHADONE',
                   'POSSESSION OF METHADONE FOR SALES',
                   'POSSESSION OF METHADONE']
In [172]:
hallu_features=['TRANSPORTATION OF OPIATES',
                'SALE OF HALLUCINOGENIC',
                'POSSESSION OF OPIUM',
                'POSSESSION OF OPIUM DERIVATIVE',
                'POSSESSION OF OPIUM',
                'SALE OF OPIUM',
                'SALE OF OPIUM DERIVATIVE',
                'TRANSPORTATION OF OPIATES',
                'POSSESSION OF OPIUM FOR SALES',
                'TRANSPORTATION OF HALLUCINOGENIC',
                'POSSESSION OF OPIUM DERIVATIVE FOR SALES',
                'SALE OF OPIATES',
                'SALE OF HALLUCINOGENIC',
                'POSSESSION OF OPIUM DERIVATIVE',
                'POSSESSION OF OPIUM',
                'POSSESSION OF OPIATES FOR SALES',
                'POSSESSION OF HALLUCINOGENIC FOR SALES',
                'POSSESSION OF OPIATES',
                'POSSESSION OF HALLUCINOGENIC']
In [173]:
meth_features=['TRANSPORTATION OF AMPHETAMINE',
               'SALE OF AMPHETAMINE',
               'POSSESSION OF AMPHETAMINE',
               'SALE OF METH-AMPHETAMINE',
               'TRANSPORTATION OF METH-AMPHETAMINE',
               'POSSESSION OF AMPHETAMINE FOR SALES',
               'POSSESSION OF METH-AMPHETAMINE FOR SALE',
               'POSSESSION OF METH-AMPHETAMINE']
In [175]:
heroin_features=['SALE OF HEROIN',
                 'POSSESSION OF HEROIN',
                 'POSSESSION OF HEROIN FOR SALES',
                 'TRANSPORTATION OF HEROIN',
                 'SALE OF HEROIN',
                 'POSSESSION OF HEROIN FOR SALES',
                 'POSSESSION OF HEROIN']
In [174]:
crack_features=['POSSESSION OF BASE/ROCK COCAINE FOR SALE',
                'SALE OF BASE/ROCK COCAINE',
                'POSSESSION OF BASE/ROCK COCAINE']
In [283]:
# Lets use real dates for plotting
days_from_start=pd.Series(t_all.index*30).astype('timedelta64[D]')
dates_for_plot=date.min()+days_from_start
time_labels=dates_for_plot.map(lambda x: str(x.year)+'-'+str(x.month))
In [397]:
def drug_analysis(t,district,plot):
    t['BARBITUATES']=t[map(lambda s: s.strip(), barituate_features)].sum(axis=1)
    t['HEROIN']=t[map(lambda s: s.strip(), heroin_features)].sum(axis=1)
    t['HALLUCINOGENIC']=t[map(lambda s: s.strip(), hallu_features)].sum(axis=1)
    t['METH']=t[map(lambda s: s.strip(), meth_features)].sum(axis=1)
    t['WEED']=t[map(lambda s: s.strip(), weed_features)].sum(axis=1)
    t['COKE']=t[map(lambda s: s.strip(), coke_features)].sum(axis=1)
    t['METHADONE']=t[map(lambda s: s.strip(), metadone_features)].sum(axis=1)
    t['CRACK']=t[map(lambda s: s.strip(), crack_features)].sum(axis=1)
    drugs=t[['HEROIN','HALLUCINOGENIC','METH','WEED','COKE','CRACK']]
    if plot:
        drugs.index=[int(i) for i in drugs.index]
        colors = plt.cm.jet(np.linspace(0, 1, drugs.shape[1]))
        drugs.plot(kind='bar', stacked=True, figsize=(20,5), color=colors, width=1, title=district,fontsize=6)
    return drugs

drug_df_all=drug_analysis(t_all,district,True)
In [398]:
def drug_analysis_rescale(t,district,plot):
    t['BARBITUATES']=t[map(lambda s: s.strip(), barituate_features)].sum(axis=1)
    t['HEROIN']=t[map(lambda s: s.strip(), heroin_features)].sum(axis=1)
    t['HALLUCINOGENIC']=t[map(lambda s: s.strip(), hallu_features)].sum(axis=1)
    t['METH']=t[map(lambda s: s.strip(), meth_features)].sum(axis=1)
    t['WEED']=t[map(lambda s: s.strip(), weed_features)].sum(axis=1)
    t['COKE']=t[map(lambda s: s.strip(), coke_features)].sum(axis=1)
    t['METHADONE']=t[map(lambda s: s.strip(), metadone_features)].sum(axis=1)
    t['CRACK']=t[map(lambda s: s.strip(), crack_features)].sum(axis=1)
    drugs=t[['HEROIN','HALLUCINOGENIC','METH','WEED','COKE','CRACK']]
    if plot:
        drugs=drugs.div(drugs.sum(axis=1),axis=0)
        drugs.index=[int(i) for i in drugs.index]
        colors = plt.cm.GnBu(np.linspace(0, 1, drugs.shape[1]))
        colors = plt.cm.jet(np.linspace(0, 1, drugs.shape[1]))
        drugs.plot(kind='bar', stacked=True, figsize=(20,5), color=colors, width=1, title=district, legend=False)
        plt.ylim([0,1])
    return drugs

drug_df_all=drug_analysis_rescale(t_all,district,True)

Let's add the real dates.

In [375]:
dates_for_plot.index=dates_for_plot
sns.set_context(rc={"figure.figsize": (12.5,5.5)})
for d,c in zip(['METH','CRACK','HEROIN','WEED'],['b','r','c','g']):
    plt.plot(dates_for_plot.index,drug_df_all[d],'o-',color=c,ms=6,mew=1.5,mec='white',linewidth=0.5,label=d,alpha=0.75)
plt.legend(loc='upper left',scatterpoints=1,prop={'size':8})
Out[375]:
<matplotlib.legend.Legend at 0x156c6f290>