Skip to content

Commit db26fcc

Browse files
stefancoestefancoe
committedFeb 8, 2017
adding code for isochrone project.
1 parent cf61519 commit db26fcc

File tree

4 files changed

+549
-0
lines changed

4 files changed

+549
-0
lines changed
 

‎isochrones/leaflet/Isochrone.htm

+94
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
<!DOCTYPE html>
2+
<html>
3+
<head>
4+
5+
<title>Choropleth Tutorial - Leaflet</title>
6+
7+
<meta charset="utf-8" />
8+
<meta name="viewport" content="width=device-width, initial-scale=1.0">
9+
10+
<link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />
11+
12+
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.3/dist/leaflet.css" />
13+
<script src="https://unpkg.com/leaflet@1.0.3/dist/leaflet.js"></script>
14+
<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.1.1/jquery.min.js"></script>
15+
16+
17+
<style>
18+
#map {
19+
width: 600px;
20+
height: 400px;
21+
}
22+
</style>
23+
24+
<style>#map { width: 800px; height: 500px; }
25+
.info { padding: 6px 8px; font: 14px/16px Arial, Helvetica, sans-serif; background: white; background: rgba(255,255,255,0.8); box-shadow: 0 0 15px rgba(0,0,0,0.2); border-radius: 5px; } .info h4 { margin: 0 0 5px; color: #777; }
26+
.legend { text-align: left; line-height: 18px; color: #555; } .legend i { width: 18px; height: 18px; float: left; margin-right: 8px; opacity: 0.7; }</style>
27+
</head>
28+
<body>
29+
30+
<div id='map'></div>
31+
32+
33+
34+
<script type="text/javascript">
35+
36+
var map = L.map('map').setView([48.25014782463433, -122.3304104651998], 12);
37+
//var map = L.map('map').setView([-122.394, 47.833], 12);
38+
39+
40+
L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpandmbXliNDBjZWd2M2x6bDk3c2ZtOTkifQ._QA7i5Mpkd_m30IGElHziw', {
41+
maxZoom: 18,
42+
attribution: 'Map data &copy; <a href="http://openstreetmap.org">OpenStreetMap</a> contributors, ' +
43+
'<a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>, ' +
44+
'Imagery © <a href="http://mapbox.com">Mapbox</a>',
45+
id: 'mapbox.light'
46+
}).addTo(map);
47+
48+
//$.getJSON(".parcels2.geojson", function(data) {
49+
// var geojson = L.geoJson(data, {
50+
// onEachFeature: function (feature, layer) {
51+
// layer.bindPopup(feature.properties.psrc_id);
52+
// }
53+
// });
54+
// geojson.addTo(map);
55+
//});
56+
57+
58+
var n2k_dh_geojson = L.geoJson(null);
59+
60+
map.on('click', function (e){
61+
//n2k_dh_geojson.clearLayers();
62+
map.removeLayer(n2k_dh_geojson);
63+
var url = "http://localhost:5000/get_isochrone/" + String(e.latlng.lng) + "/" + String(e.latlng.lat);
64+
console.log(url);
65+
66+
$.ajax({
67+
type: "GET",
68+
//url: "http://localhost:5000/todo/get_parcels/-122.394/47.833",
69+
url: url,
70+
dataType: 'json',
71+
cache: false,
72+
success: function (response) {
73+
n2k_dh_geojson = L.geoJson(response, {
74+
onEachFeature: function (feature, layer) {
75+
//console.log(feature.properties.psrc_id);
76+
layer.bindPopup (String(feature.properties.psrc_id));
77+
}
78+
});
79+
n2k_dh_geojson.addTo(map);
80+
}
81+
});
82+
});
83+
84+
85+
86+
87+
88+
89+
</script>
90+
91+
92+
93+
</body>
94+
</html>

‎isochrones/leaflet/Parcels.htm

+94
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
<!DOCTYPE html>
2+
<html>
3+
<head>
4+
5+
<title>Choropleth Tutorial - Leaflet</title>
6+
7+
<meta charset="utf-8" />
8+
<meta name="viewport" content="width=device-width, initial-scale=1.0">
9+
10+
<link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />
11+
12+
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.3/dist/leaflet.css" />
13+
<script src="https://unpkg.com/leaflet@1.0.3/dist/leaflet.js"></script>
14+
<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.1.1/jquery.min.js"></script>
15+
16+
17+
<style>
18+
#map {
19+
width: 600px;
20+
height: 400px;
21+
}
22+
</style>
23+
24+
<style>#map { width: 800px; height: 500px; }
25+
.info { padding: 6px 8px; font: 14px/16px Arial, Helvetica, sans-serif; background: white; background: rgba(255,255,255,0.8); box-shadow: 0 0 15px rgba(0,0,0,0.2); border-radius: 5px; } .info h4 { margin: 0 0 5px; color: #777; }
26+
.legend { text-align: left; line-height: 18px; color: #555; } .legend i { width: 18px; height: 18px; float: left; margin-right: 8px; opacity: 0.7; }</style>
27+
</head>
28+
<body>
29+
30+
<div id='map'></div>
31+
32+
33+
34+
<script type="text/javascript">
35+
36+
var map = L.map('map').setView([48.25014782463433, -122.3304104651998], 12);
37+
//var map = L.map('map').setView([-122.394, 47.833], 12);
38+
39+
40+
L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpandmbXliNDBjZWd2M2x6bDk3c2ZtOTkifQ._QA7i5Mpkd_m30IGElHziw', {
41+
maxZoom: 18,
42+
attribution: 'Map data &copy; <a href="http://openstreetmap.org">OpenStreetMap</a> contributors, ' +
43+
'<a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>, ' +
44+
'Imagery © <a href="http://mapbox.com">Mapbox</a>',
45+
id: 'mapbox.light'
46+
}).addTo(map);
47+
48+
//$.getJSON(".parcels2.geojson", function(data) {
49+
// var geojson = L.geoJson(data, {
50+
// onEachFeature: function (feature, layer) {
51+
// layer.bindPopup(feature.properties.psrc_id);
52+
// }
53+
// });
54+
// geojson.addTo(map);
55+
//});
56+
57+
58+
var n2k_dh_geojson = L.geoJson(null);
59+
60+
map.on('click', function (e){
61+
//n2k_dh_geojson.clearLayers();
62+
map.removeLayer(n2k_dh_geojson);
63+
var url = "http://localhost:5000/get_parcels/" + String(e.latlng.lng) + "/" + String(e.latlng.lat);
64+
console.log(url);
65+
66+
$.ajax({
67+
type: "GET",
68+
//url: "http://localhost:5000/todo/get_parcels/-122.394/47.833",
69+
url: url,
70+
dataType: 'json',
71+
cache: false,
72+
success: function (response) {
73+
n2k_dh_geojson = L.geoJson(response, {
74+
onEachFeature: function (feature, layer) {
75+
//console.log(feature.properties.psrc_id);
76+
layer.bindPopup (String(feature.properties.psrc_id));
77+
}
78+
});
79+
n2k_dh_geojson.addTo(map);
80+
}
81+
});
82+
});
83+
84+
85+
86+
87+
88+
89+
</script>
90+
91+
92+
93+
</body>
94+
</html>

‎isochrones/scripts/app.py

+199
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,199 @@
1+
import pandana as pdna
2+
import pandas as pd
3+
import numpy as np
4+
import os
5+
import re
6+
import sys
7+
from pyproj import Proj, transform
8+
from flask import Flask, jsonify
9+
from flask import abort
10+
import json
11+
import geojson
12+
import sqlite3
13+
from shapely.geometry import Point, MultiPoint
14+
from geopandas import GeoDataFrame, GeoSeries
15+
16+
#os.chdir(r'D:\Stefan\Isochrone')
17+
app = Flask(__name__)
18+
19+
parcels_file_name = r'D:\Stefan\Isochrone\repository\data\parcels_urbansim.txt'
20+
nodes_file_name = r'D:\Stefan\Isochrone\repository\data\all_streets_nodes_2014.csv'
21+
links_file_name = r'D:\Stefan\Isochrone\repository\data\all_streets_links_2014.csv'
22+
23+
max_dist = 5280 # 3 miles
24+
parcels = pd.DataFrame.from_csv(parcels_file_name, sep = " ", index_col = None )
25+
#parcels = parcels.columns.map(lambda x: x.upper())
26+
27+
# nodes must be indexed by node_id column, which is the first column
28+
nodes = pd.DataFrame.from_csv(nodes_file_name)
29+
links = pd.DataFrame.from_csv(links_file_name, index_col = None )
30+
31+
# get rid of circular links
32+
links = links.loc[(links.from_node_id <> links.to_node_id)]
33+
34+
# assign impedance
35+
imp = pd.DataFrame(links.Shape_Length)
36+
imp = imp.rename(columns = {'Shape_Length':'distance'})
37+
38+
# create pandana network
39+
network = pdna.network.Network(nodes.x, nodes.y, links.from_node_id, links.to_node_id, imp)
40+
41+
def assign_nodes_to_dataset(dataset, network, column_name, x_name, y_name):
42+
"""Adds an attribute node_ids to the given dataset."""
43+
dataset[column_name] = network.get_node_ids(dataset[x_name].values, dataset[y_name].values)
44+
45+
assign_nodes_to_dataset(parcels, network, 'node_ids', 'xcoord_p', 'ycoord_p')
46+
network.init_pois(1, max_dist, 1)
47+
48+
#class MyEncoder(json.JSONEncoder):
49+
# def default(self, obj):
50+
# if isinstance(obj, numpy.integer):
51+
# return int(obj)
52+
# elif isinstance(obj, numpy.floating):
53+
# return float(obj)
54+
# elif isinstance(obj, numpy.ndarray):
55+
# return obj.tolist()
56+
# else:
57+
# return super(MyEncoder, self).default(obj)
58+
59+
60+
61+
def dict_factory(cursor, row):
62+
d = {}
63+
for idx,col in enumerate(cursor.description):
64+
d[col[0]] = row[idx]
65+
return d
66+
def reproject_to_state_plane(longitude, latitude, ESPG = "+init=EPSG:2926", conversion = 3.28084):
67+
68+
# Remember long is x and lat is y!
69+
prj_wgs = Proj(init='epsg:4326')
70+
prj_wgs.is_latlong()
71+
prj_sp = Proj(init='EPSG:2926', preserve_units = True)
72+
73+
x, y = transform(prj_wgs, prj_sp, longitude, latitude)
74+
75+
return x, y
76+
77+
78+
@app.route('/get_parcels/<string:x>/<string:y>', methods=['GET'])
79+
def get_parcels(x,y):
80+
print 'Getting data'
81+
#Get coords from url
82+
x = float(x)
83+
y = float(y)
84+
#Go from lat, long to state plane
85+
coords = reproject_to_state_plane(x, y)
86+
x = pd.Series(coords[0])
87+
y = pd.Series(coords[1])
88+
#Set as a point of interest on the pandana network
89+
network.set_pois('tstop', x, y)
90+
#Find distance to stop from all nodes, everything over a mile gets a value of 99999
91+
res = network.nearest_pois(max_dist, 'tstop', num_pois=1, max_distance=99999)
92+
res[res <> 999] = (res[res <> 99999]/5280.).astype(res.dtypes) # convert to miles
93+
res_name = "dist_tstop"
94+
df_parcels = parcels
95+
df_parcels[res_name] = res.loc[df_parcels.node_ids].values
96+
df_parcels = df_parcels.loc[(df_parcels.dist_tstop<99999)]
97+
parcel_id_list = df_parcels.parcelid.values.astype(int).tolist()
98+
parcel_id_string = ''
99+
parcel_id_string = ",".join(map(str,parcel_id_list))
100+
parcel_id_string = '(' + parcel_id_string + ')'
101+
102+
psqlite='D:/Stefan/Isochrone/sqlite'
103+
os.environ["PATH"] = psqlite + os.pathsep + os.environ["PATH"]
104+
#no import sqlite3
105+
106+
#connect to a sqlite database
107+
dbpath=os.path.abspath(r'D:\Stefan\Isochrone\repository\data\parcels.sqlite')
108+
conn=sqlite3.connect(dbpath)
109+
#load the spatialite extension - if everything is fine, you should not get any errors
110+
conn.enable_load_extension(True)
111+
conn.execute("SELECT load_extension('mod_spatialite.dll')")
112+
# apply the function to the sqlite3 enginemy
113+
conn.row_factory = dict_factory
114+
115+
drop_table_query = """drop table if exists parcels1"""
116+
conn.execute(drop_table_query)
117+
118+
create_table_query = """create table parcels1 as select TAZ, psrc_id, st_unaryunion(st_collect(geometry)) as geometry from parcels WHERE psrc_id in %s group by TAZ""" % parcel_id_string
119+
conn.execute(create_table_query)
120+
121+
getResultsQuery = """SELECT AsGeoJSON(geometry), psrc_id FROM parcels1"""
122+
#print getResultsQuery
123+
124+
125+
# fetch the results in form of a list of dictionaries
126+
results = conn.execute(getResultsQuery).fetchall()
127+
#print results
128+
# create a new list which will store the single GeoJSON features
129+
featureCollection = []
130+
131+
# iterate through the list of result dictionaries
132+
for row in results:
133+
134+
# create a single GeoJSON geometry from the geometry column which already contains a GeoJSON string
135+
geom = geojson.loads(row['AsGeoJSON(geometry)'])
136+
137+
# remove the geometry field from the current's row's dictionary
138+
row.pop('AsGeoJSON(geometry)')
139+
140+
# create a new GeoJSON feature and pass the geometry columns as well as all remaining attributes which are stored in the row dictionary
141+
feature = geojson.Feature(geometry=geom, properties=row)
142+
143+
# append the current feature to the list of all features
144+
featureCollection.append(feature)
145+
146+
# when single features for each row from the database table are created, pass the list to the FeatureCollection constructor which will merge them together into one object
147+
featureCollection = geojson.FeatureCollection(featureCollection)
148+
# print the FeatureCollection as string
149+
150+
GeoJSONFeatureCollectionAsString = geojson.dumps(featureCollection)
151+
print 'Done getting data'
152+
#print GeoJSONFeatureCollectionAsString
153+
return GeoJSONFeatureCollectionAsString
154+
# return featureCollection
155+
156+
@app.route('/get_isochrone/<string:x>/<string:y>', methods=['GET'])
157+
def get_isochrone(x,y):
158+
print 'Getting data'
159+
#Get coords from url
160+
x = float(x)
161+
y = float(y)
162+
#Go from lat, long to state plane
163+
coords = reproject_to_state_plane(x, y)
164+
x = pd.Series(coords[0])
165+
y = pd.Series(coords[1])
166+
#Set as a point of interest on the pandana network
167+
network.set_pois('tstop', x, y)
168+
#Find distance to stop from all nodes, everything over a mile gets a value of 99999
169+
res = network.nearest_pois(max_dist, 'tstop', num_pois=1, max_distance=99999)
170+
res[res <> 999] = (res[res <> 99999]/5280.).astype(res.dtypes) # convert to miles
171+
res_name = "dist_tstop"
172+
df_parcels = parcels
173+
df_parcels[res_name] = res.loc[df_parcels.node_ids].values
174+
df_parcels = df_parcels.loc[(df_parcels.dist_tstop<99999)]
175+
176+
geometry = [(xy) for xy in zip (df_parcels.xcoord_p, df_parcels.ycoord_p)]
177+
geo_series = GeoSeries(MultiPoint(geometry))
178+
geo_series.crs = {'init' :'epsg:2285'}
179+
poly = geo_series.geometry.convex_hull
180+
poly.crs = {'init' :'epsg:2285'}
181+
poly = poly.to_crs({'init' :'epsg:4326'})
182+
g2 = geojson.Feature(geometry=poly.geometry[0], properites={})
183+
#f = open('d:/poly4.geojson', 'w')
184+
#geojson.dump(g2, f)
185+
#f.close()
186+
featureCollection = []
187+
featureCollection.append(poly)
188+
featureCollection = geojson.FeatureCollection(featureCollection)
189+
GeoJSONFeatureCollectionAsString = geojson.dumps(featureCollection)
190+
print 'Done getting data'
191+
#print GeoJSONFeatureCollectionAsString
192+
return GeoJSONFeatureCollectionAsString
193+
194+
195+
196+
197+
198+
if __name__ == '__main__':
199+
app.run(debug=True)

‎isochrones/scripts/test.py

+162
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
import pandana as pdna
2+
import pandas as pd
3+
import numpy as np
4+
import os
5+
import re
6+
import sys
7+
from pyproj import Proj, transform
8+
from flask import Flask, jsonify
9+
from flask import abort
10+
import json
11+
import geojson
12+
import sqlite3
13+
import shapely
14+
from shapely.geometry import Point, MultiPoint
15+
from geopandas import GeoDataFrame, GeoSeries
16+
17+
#os.chdir(r'D:\Stefan\Isochrone')
18+
#app = Flask(__name__)
19+
20+
parcels_file_name = r'D:\Stefan\Isochrone\parcels_urbansim.txt'
21+
nodes_file_name = r'D:\Stefan\Isochrone\all_streets_nodes_2014.csv'
22+
links_file_name = r'D:\Stefan\Isochrone\all_streets_links_2014.csv'
23+
24+
max_dist = 5280 # 3 miles
25+
parcels = pd.DataFrame.from_csv(parcels_file_name, sep = " ", index_col = None )
26+
#parcels = parcels.columns.map(lambda x: x.upper())
27+
28+
# nodes must be indexed by node_id column, which is the first column
29+
nodes = pd.DataFrame.from_csv(nodes_file_name)
30+
links = pd.DataFrame.from_csv(links_file_name, index_col = None )
31+
32+
# get rid of circular links
33+
links = links.loc[(links.from_node_id <> links.to_node_id)]
34+
35+
# assign impedance
36+
imp = pd.DataFrame(links.Shape_Length)
37+
imp = imp.rename(columns = {'Shape_Length':'distance'})
38+
39+
# create pandana network
40+
network = pdna.network.Network(nodes.x, nodes.y, links.from_node_id, links.to_node_id, imp)
41+
42+
def assign_nodes_to_dataset(dataset, network, column_name, x_name, y_name):
43+
"""Adds an attribute node_ids to the given dataset."""
44+
dataset[column_name] = network.get_node_ids(dataset[x_name].values, dataset[y_name].values)
45+
46+
assign_nodes_to_dataset(parcels, network, 'node_ids', 'xcoord_p', 'ycoord_p')
47+
network.init_pois(1, max_dist, 1)
48+
49+
#class MyEncoder(json.JSONEncoder):
50+
# def default(self, obj):
51+
# if isinstance(obj, numpy.integer):
52+
# return int(obj)
53+
# elif isinstance(obj, numpy.floating):
54+
# return float(obj)
55+
# elif isinstance(obj, numpy.ndarray):
56+
# return obj.tolist()
57+
# else:
58+
# return super(MyEncoder, self).default(obj)
59+
60+
61+
62+
def dict_factory(cursor, row):
63+
d = {}
64+
for idx,col in enumerate(cursor.description):
65+
d[col[0]] = row[idx]
66+
return d
67+
68+
69+
70+
#@app.route('/todo/get_parcels/<int:x>/<int:y>', methods=['GET'])
71+
#def get_parcels(x,y):
72+
x = pd.Series([1265627])
73+
y = pd.Series([255656])
74+
network.set_pois('tstop', x, y)
75+
res = network.nearest_pois(max_dist, 'tstop', num_pois=1, max_distance=99999)
76+
res[res <> 999] = (res[res <> 99999]/5280.).astype(res.dtypes) # convert to miles
77+
res_name = "dist_tstop"
78+
test = parcels
79+
test[res_name] = res.loc[test.node_ids].values
80+
test = test.loc[(test.dist_tstop<99999)]
81+
82+
83+
84+
geometry = [(xy) for xy in zip (test.xcoord_p, test.ycoord_p)]
85+
geo_series = GeoSeries(MultiPoint(geometry))
86+
geo_series.crs = {'init' :'epsg:2285'}
87+
poly = geo_series.geometry.convex_hull
88+
poly.crs = {'init' :'epsg:2285'}
89+
poly = poly.to_crs({'init' :'epsg:4326'})
90+
g2 = geojson.Feature(geometry=poly.geometry[0], properites={})
91+
f = open('d:/poly4.geojson', 'w')
92+
geojson.dump(g2, f)
93+
f.close()
94+
95+
96+
97+
98+
geometry = MultiPoint(geometry)
99+
poly = geometry.convex_hull
100+
g2 = geojson.Feature(geometry=poly, properites={})
101+
102+
geo_df = GeoDataFrame(df_parcels, crs=credits, geometry = geometry)
103+
geo_series = GeoSeries(MultiPoint(geometry))
104+
105+
106+
107+
108+
test = ''
109+
test = ",".join(map(str,x))
110+
test = '(' + test + ')'
111+
#x = test.parcelid.to_json
112+
#return json.dumps({'parcels': x})
113+
#return json.dumps(x)
114+
#return {'key': test.parcelid.to_dict(orient='record')}
115+
psqlite='D:/Stefan/Isochrone/sqlite'
116+
os.environ["PATH"] = psqlite + os.pathsep + os.environ["PATH"]
117+
#no import sqlite3
118+
119+
#connect to a sqlite database
120+
dbpath=os.path.abspath(r'D:\Stefan\Isochrone\parcels.sqlite')
121+
conn=sqlite3.connect(dbpath)
122+
#load the spatialite extension - if everything is fine, you should not get any errors
123+
conn.enable_load_extension(True)
124+
conn.execute("SELECT load_extension('mod_spatialite.dll')")
125+
# apply the function to the sqlite3 enginemy
126+
conn.row_factory = dict_factory
127+
128+
getResultsQuery = """SELECT AsGeoJSON(geometry), psrc_id FROM parcels WHERE psrc_id in %s""" % test
129+
#print getResultsQuery
130+
131+
132+
# fetch the results in form of a list of dictionaries
133+
results = conn.execute(getResultsQuery).fetchall()
134+
#print results
135+
# create a new list which will store the single GeoJSON features
136+
featureCollection = []
137+
138+
# iterate through the list of result dictionaries
139+
for row in results:
140+
141+
# create a single GeoJSON geometry from the geometry column which already contains a GeoJSON string
142+
geom = geojson.loads(row['AsGeoJSON(geometry)'])
143+
144+
# remove the geometry field from the current's row's dictionary
145+
row.pop('AsGeoJSON(geometry)')
146+
147+
# create a new GeoJSON feature and pass the geometry columns as well as all remaining attributes which are stored in the row dictionary
148+
feature = geojson.Feature(geometry=geom, properties=row)
149+
150+
# append the current feature to the list of all features
151+
print feature
152+
featureCollection.append(feature)
153+
154+
# when single features for each row from the database table are created, pass the list to the FeatureCollection constructor which will merge them together into one object
155+
featureCollection = geojson.FeatureCollection(featureCollection)
156+
# print the FeatureCollection as string
157+
158+
GeoJSONFeatureCollectionAsString = geojson.dumps(featureCollection)
159+
#return GeoJSONFeatureCollectionAsString
160+
f = open('d:/parcels2.geojson', 'w')
161+
geojson.dump(featureCollection, f)
162+
f.close()

0 commit comments

Comments
 (0)
Please sign in to comment.