<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Uncensored chicago data
After asking around the research community, it appears that the chicago dataset has been re-geocoded sometime in the last year (or so).  In the past, each event was geocoded to a fairly exact address (perhaps building entrance).
A little searching can find copies of extracts of this data.  As there must have been a good reason for changing the geocoding (presumably privacy based) we will not release a copy of this data, and we will be careful when drawing maps etc. not to show more detail than can already by found in the literature.
However, it is interesting to compare the data, and gain some insight into how the new geocoding is performed.
Step1:  Load the uncensored data
As usual, we'll drop the "HOMICIDE" crimes.  
The case number "HV612588" is reported twice, and one looks like an error.  We'll drop in
Step2:  Load matching "new" data
We load the same case numbers from the dataset (downloaded in May 2017), perform a little manual cleaning, and check that we then have a unique set of case numbers.
We are now mising about 1000 cases from the "old" data.
Step3:  Check the time range
We now look at the current file, and find all events in the same time range.  There are about 1000 more; so very good agreement overall.
Step4:  Check for any differences, except in geometry
Turns out there are quite a few differences.
A lot of differences in "crime", and a huge number in "type".  I have not looked closely, but from a brief look, it seems that the differences are mostly down to reasonable re-classification.
Step5:  Compute difference in distance between the points
In lon/lat space, which is areal, but not unreasonable on a small scale like this.
We then visualise, with a log scale.  So almost all points are between 10e-3 and 10e-4 distance apart.  I.e. most move a bit, but not too much.
Step6:  Assemble the points into a geoDataFrame
Step7:  Look at one block
We'll pick some blocks at random (without giving their names, to protect privacy) and plot the difference in the geocoding as lines joining the old and new points.
What we see is that
Step8:  Look at the TIGER data
The TIGER/Line dataset from the US Census bureau provides a great resource for getting block-level address data.  See the other notebook for more details.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
import numpy as np
import open_cp.sources.chicago as chicago
import os, lzma
filename = os.path.join("..", "..", "..", "..", "Data", "two.csv")
file_all = os.path.join("..", "..", "open_cp", "sources", "chicago_all.csv.xz")
Explanation: Uncensored chicago data
After asking around the research community, it appears that the chicago dataset has been re-geocoded sometime in the last year (or so).  In the past, each event was geocoded to a fairly exact address (perhaps building entrance).
A little searching can find copies of extracts of this data.  As there must have been a good reason for changing the geocoding (presumably privacy based) we will not release a copy of this data, and we will be careful when drawing maps etc. not to show more detail than can already by found in the literature.
However, it is interesting to compare the data, and gain some insight into how the new geocoding is performed.
End of explanation
# Load the uncensored data
data = chicago.load_to_geoDataFrame(filename, type="all_other")
data = data[data.case != "HV612588"]
data = data[data.crime != "HOMICIDE"]
len(data)
# So every row has a unique case
assert len(data.case.unique()) == len(data)
data = data.set_index("case")
Explanation: Load the uncensored data
As usual, we'll drop the "HOMICIDE" crimes.  
The case number "HV612588" is reported twice, and one looks like an error.  We'll drop in
End of explanation
cases = set(data.index)
with lzma.open(file_all, "rt") as file:
    data_all = [event for event in chicago.generate_GeoJSON_Features(file, type="all")
               if event["properties"]["case"] in cases]
data_all = gpd.GeoDataFrame.from_features(data_all)
len(data_all), len(data_all.case.unique())
data_all = data_all[(data_all.case != "") |
     ((data_all.case == "") & (data_all.crime == "BURGLARY"))]
data_all = data_all[(data_all.case != "HT155546") |
    ((data_all.case == "HT155546") & (data_all.location == "TAVERN/LIQUOR STORE"))]
len(data_all), len(data_all.case.unique())
data_all = data_all.set_index("case")
len(set(data.index) - set(data_all.index)), len(set(data_all.index) - set(data.index))
Explanation: Load matching "new" data
We load the same case numbers from the dataset (downloaded in May 2017), perform a little manual cleaning, and check that we then have a unique set of case numbers.
We are now mising about 1000 cases from the "old" data.
End of explanation
times = data.timestamp.map(chicago._date_from_iso)
start, end = min(times), max(times)
start, end
with lzma.open(file_all, "rt") as file:
    data_all_in_time = []
    for event in chicago.generate_GeoJSON_Features(file, type="all"):
        ts = chicago._date_from_iso(event["properties"]["timestamp"])
        if ts >= start and ts <= end:
            data_all_in_time.append(event)
data_all_in_time = gpd.GeoDataFrame.from_features(data_all_in_time)
len(set(data_all_in_time.case))
data_all_in_time = None
Explanation: Check the time range
We now look at the current file, and find all events in the same time range.  There are about 1000 more; so very good agreement overall.
End of explanation
both = pd.merge(data, data_all, left_index=True, right_index=True)
print("'address' differences", len(both[~(both.address_x == both.address_y)]))
print("'crime' differences", len(both[~(both.crime_x == both.crime_y)]))
print("'type' differences", len(both[~(both.type_x == both.type_y)]))
print("'location' differences", len(both[~(both.location_x == both.location_y)]))
print("'timestamp' differences", len(both[~(both.timestamp_x == both.timestamp_y)]))
Explanation: Check for any differences, except in geometry
Turns out there are quite a few differences.
A lot of differences in "crime", and a huge number in "type".  I have not looked closely, but from a brief look, it seems that the differences are mostly down to reasonable re-classification.
End of explanation
dists = []
for _, row in both.iterrows():
    if row.geometry_x is not None and row.geometry_y is not None:
        dists.append(row.geometry_x.distance(row.geometry_y))
    else:
        dists.append(-1)
both["distance"] = dists
x = np.linspace(2, 5, 100)
xx = np.exp(-x * np.log(10)) 
y = [len(both[both["distance"] > t]) for t in xx]
fig, ax = plt.subplots()
ax.scatter(x,y)
Explanation: Compute difference in distance between the points
In lon/lat space, which is areal, but not unreasonable on a small scale like this.
We then visualise, with a log scale.  So almost all points are between 10e-3 and 10e-4 distance apart.  I.e. most move a bit, but not too much.
End of explanation
tmp = both.ix[both[both["distance"] < 0.001].index]
compare = gpd.GeoDataFrame()
compare["geo1"] = tmp["geometry_x"]
compare["geo2"] = tmp["geometry_y"]
compare["name1"] = tmp["address_x"]
compare["name2"] = tmp["address_y"]
compare.set_index(tmp.index)
tmp = None
compare = compare[compare.geo1.map(lambda pt : pt is not None)]
compare = compare[compare.geo2.map(lambda pt : pt is not None)]
import shapely.geometry
compare["geometry"] = [ shapely.geometry.LineString([ row.geo1, row.geo2 ])
    for _, row in compare.iterrows() ]
compare = compare.set_geometry("geometry")
compare.crs = {"init": "epsg:4326"}
import os
tmp = compare.drop(["geo1", "geo2"], axis=1)
tmp.to_file( os.path.join("..", "..", "..", "..", "Data", "Old_Chicago_SHP") )
tmp = None
bds = compare.total_bounds
lines = []
for geo in compare.geometry:
    xs, ys = geo.xy
    lines.append(list(zip(xs, ys)))
from matplotlib.collections import LineCollection
lc = LineCollection(lines, linewidths=0.1)
fig, ax = plt.subplots(figsize=(16,16))
ax.add_collection(lc)
ax.set(xlim=[bds[0], bds[2]], ylim=[bds[1], bds[3]])
None
Explanation: Assemble the points into a geoDataFrame
End of explanation
fig, axes = plt.subplots(ncols=5, figsize=(18,10))
blocks = np.random.randint(low=0, high=10000, size=5)
blocks = [3243, 7016, 6951, 5475,  804]
for i, ax in zip(blocks, axes):
    block = compare.ix[i].name1
    compare[compare.name1 == block].plot(ax=ax, color="black")
    ax.set_aspect(1)
    ax.axes.xaxis.set_visible(False)
    ax.axes.yaxis.set_visible(False)
Explanation: Look at one block
We'll pick some blocks at random (without giving their names, to protect privacy) and plot the difference in the geocoding as lines joining the old and new points.
What we see is that:
- Each point appears to be projected (orthogonally) to a line.  (Later we'll see that this line comes from the TIGER/Line dataset from the US Census bureau.)
- Then the points are squashed together.
- This has the effect of preserving the ordering and relative position along the block
- As you can see the last two plots, a single block might be split into more than one "edge", and this procedure appears to be applied to each "edge" independently.
We hence very clearly see how the pattern observed in the new/current dataset is generated: we will naturally end up with clumps of points in, roughly, the middle of the centre of the road of each block.
End of explanation
filename = os.path.join("/media", "disk", "TIGER Data", "tl_2013_17031_edges", "tl_2013_17031_edges.shp")
edges = gpd.read_file(filename)
want = {"geometry", "FULLNAME", "LFROMADD", "LTOADD", "RFROMADD", "RTOADD"}
edges = gpd.GeoDataFrame({key:edges[key] for key in want})
edges.crs={'init': 'epsg:4269'}
import rtree
gap = 0.001
def gen():
    for i, row in edges.iterrows():
        bds = list(row.geometry.bounds)
        bds = [bds[0]-gap, bds[1]-gap, bds[2]+gap, bds[3]+gap]
        yield i, bds, None
        
idx = rtree.index.Index(gen())
def find_match_via_rtree(point):
    possibles = edges.ix[list(idx.intersection(point.coords[0]))]
    if len(possibles) == 0:
        #raise ValueError("Found no candidates for {}".format(point))
        from collections import namedtuple
        Error = namedtuple("Error", ["name"])
        return Error(name=-1)
    i = possibles.geometry.distance(point).argmin()
    return edges.ix[i]
block = compare.ix[804].name1
block = compare[compare.name1 == block]
rows = set( find_match_via_rtree(row.geometry).name for _, row in block.iterrows() )
rows = edges.ix[rows]
rows = rows[~rows.FULLNAME.map(lambda x : x is None or x == "Alley")]
fig, ax = plt.subplots(figsize=(16,10))
rows.plot(ax=ax)
block.plot(color="black", ax=ax)
xcs, ycs = [], []
for geo in rows.geometry:
    x, y = geo.coords[0]
    xcs.append(x)
    ycs.append(y)
    x, y = geo.coords[-1]
    xcs.append(x)
    ycs.append(y)
ax.scatter(xcs, ycs, color="blue", marker="o")
ax.set_aspect(1)
ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
Explanation: Look at the TIGER data
The TIGER/Line dataset from the US Census bureau provides a great resource for getting block-level address data.  See the other notebook for more details.
End of explanation
<END_TASK>