-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcreate_area_map
executable file
·322 lines (270 loc) · 12.3 KB
/
create_area_map
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
#!/usr/bin/env python3
from optparse import OptionParser, OptionValueError
import os, sys
import geopandas as gpd
import pandas as pd
import shapely
import folium
import folium.plugins
import geopy.geocoders
import matplotlib
import matplotlib.cm
def get_schools(boundingbox = None):
"""Return the schools in optional boundingbox we are interested in.
boundingbox is
min lat, max lat, min long, max long
ymin, ymax, xmin, xmax
"""
# these came from
# https://catalogue.data.govt.nz/dataset/directory-of-educational-institutions/resource/bdfe0e4c-1554-4701-a8fe-ba1c8e0cc2ce
schools = pd.read_csv('schooldirectory-28-04-2019-233051.csv')
# only want to keep some info (knowing the Total number and number
# of Europeans might be useful in predicting Hoity-toity ness)
schools = schools[['School_Id', 'Org_Name', 'URL', 'Longitude', 'Latitude', 'Org_Type', 'CoEd_Status', 'Decile', 'Total', 'European']]
# some schools don't have location, but I have checked these by
# eye, none of them interest us, so just drop them
schools = schools[schools.Latitude.notnull()]
# since George is going into Intermediate, lets just consider
# schools that provide education for intermediate and older
#schools = schools.loc[(schools.Org_Type == 'Full Primary') | (schools.Org_Type == 'Intermediate') | (schools.Org_Type.str.startswith('Second'))]
schools = schools.loc[(schools.Org_Type == 'Contributing') | (schools.Org_Type == 'Full Primary') | (schools.Org_Type == 'Intermediate') | (schools.Org_Type.str.startswith('Second'))]
# we need a location as a proper Point type for cx slicing to work
schools['geometry'] = [shapely.geometry.Point(x, y) for x, y in zip(schools.Longitude, schools.Latitude)]
schools = gpd.GeoDataFrame(schools)
schools.crs = {'init' :'epsg:4326'} # this is basic latitude longitude
# boundingbox [ymin, ymax, xmin, xmax], but cx slicing takes # [xmin:xmax, ymin:ymax]
if boundingbox: schools = schools.cx[boundingbox[2]:boundingbox[3], boundingbox[0]:boundingbox[1]]
# add in a column to indicate if it is a high school (easier than
# checking Org_Type startswith 'Second')
schools['high'] = schools.Org_Type.str.startswith('Second')
# some schools have zones, lets put those in for those that do.
# Got the zone info from https://koordinates.com/layer/743-nz-school-zones-sept-2010/?ex=1
zones = gpd.read_file('school_zones/nz-school-zones-sept-2010.shp')
zones = zones.to_crs(epsg=4326) # basic lat/lon, everything else is in this
# merge with schools. Do a left merge so every school is retained,
# with optional zone info
schools = schools.merge(zones[['SchoolID', 'geometry']], how='left', left_on = 'School_Id', right_on = 'SchoolID')
schools.drop(columns='SchoolID', inplace=True)
schools.rename(columns={'geometry_y': 'zone'}, inplace=True)
# the zone is a shapely polygon, for mapping we need a folium polygon
for i, row in schools.iterrows():
if pd.notnull(row.zone):
latlons = [(i[1], i[0]) for i in shapely.geometry.mapping(row.zone)['coordinates'][0]]
schools.at[i, 'zonepoly'] = folium.vector_layers.Polygon(latlons)
# each school needs a marker for showing on the map
for i, row in schools.iterrows():
marker = folium.Marker(
location = [row.Latitude, row.Longitude],
# high schools get blue buildings, otherwise green I (for intermediate)
icon = folium.Icon(color='blue', prefix='fa', icon='building') if row.high else folium.Icon(color='green',icon='italic'),
popup="""
<a href="{url}" target="_blank">{name}</a>
<table>
<tr><td>Type</td><td>{stype}</td>
<tr><td>CoEd?</td><td>{coed}</td>
<tr><td>Decile</td><td>{decile}</td>
<tr><td># Students</td><td>{numstudents}</td>
<tr><td># European </td><td>{numeuropeans}</td>
<tr><td>Zone</td><td>{zone}</td>
</table>""".format(
url = row.URL if pd.notnull(row.URL) else '',
name = row.Org_Name.replace(' ', ' '),
stype = row.Org_Type.replace(' ', ' '),
coed = 'Y' if row.CoEd_Status == 'Co-Educational' else 'N',
decile = row.Decile,
numstudents = row.Total,
numeuropeans = row.European,
sid = row.School_Id,
zone = '<a id="zonelink%s" href="#">Show</a>' % row.School_Id if pd.notnull(row.zone) else 'None')
)
schools.at[i, 'marker'] = marker
return schools
def get_metro(city_name, boundingbox):
"""Return the bus routes and stops in optional boundingbox we are interested in."""
# the route data comes from
# chch
# http://opendata.canterburymaps.govt.nz/datasets/-bus-route-directions/data
#routes = gpd.read_file('metro/_Bus_Route_Directions.shp')
# welly
# http://data-gwrc.opendata.arcgis.com/datasets/metlink-public-transport-routes
if city_name == 'Christchurch':
routes = gpd.read_file('metro/_Bus_Route_Directions.shp')
else:
routes = gpd.read_file('metro/Metlink_Public_Transport_Routes.shp')
routes = routes.to_crs(epsg=4326) # basic lat/lon, everything else is in this
# only want routes in bounds of city
# boundingbox [ymin, ymax, xmin, xmax], but cx slicing takes # [xmin:xmax, ymin:ymax]
if boundingbox: routes = routes.cx[boundingbox[2]:boundingbox[3], boundingbox[0]:boundingbox[1]]
# I found this in a css file on the metro website. Bit of a pain
# it isn't included in the data file. This maps the RouteNo to
# colour for Christchurch
no2colour = {
'A': '#F18600',
'F': '#005899',
'B': '#3ebced',
'O': '#f37021',
'P': '#554588',
'Y': '#ffc20e',
'Or': '#79bc43',
'17': '#ec008c',
'28': '#f79328',
'29': '#00539f',
'44': '#0074ad',
'45': '#96643A',
'60': '#da6fab',
'80': '#717dbd',
'95': '#c41039',
'100': '#88807e',
'107': '#46a299',
'108': '#a285b3',
'120': '#faa61a',
'125': '#5f9937',
'130': '#9f3925',
'135': '#0db14b',
'140': '#00929e',
'145': '#956438',
'150': '#9853a1',
'155': '#9C6DA8',
'535': '#d2bfa5',
'820': '#46ba7c',
'951': '#009db2',
'952': '#80ceca',
'960': '#f05a4e' }
# maps service name to index
if city_name == 'Christchurch':
routes['colour'] = [no2colour[i] for i in routes.RouteNo]
else:
foo = routes.SERVICE.unique()
servicename_to_index = {foo[i]:i for i in range(len(foo))}
viridis = matplotlib.cm.get_cmap('prism', len(routes.SERVICE.unique()))
routes['colour'] = [matplotlib.colors.rgb2hex(viridis(servicename_to_index[i])[:3]) for i in routes.SERVICE]
# for mapping we need a folium polyline
for i, row in routes.iterrows():
latlons = [(i[1],i[0]) for i in shapely.geometry.mapping(row.geometry)['coordinates']]
routes.at[i, 'polyline'] = folium.vector_layers.PolyLine(latlons, color=row.colour, popup = row.RouteNo + ' ' + row.RouteName if city_name == 'Christchurch' else row.SERVICE)
# bus stops
# chch
# http://opendata.canterburymaps.govt.nz/datasets/b595b1365f8e47618658692539d946de_1/data
# wellington
# http://data-gwrc.opendata.arcgis.com/datasets/14d628fd0890409dbd2697262fb4c49b_1
if city_name == 'Christchurch':
stops = gpd.read_file('metro/Bus_Stops.shp')
else:
stops = gpd.read_file('metro/Metlink_Service_Stops.shp')
stops = stops.to_crs(epsg=4326) # basic lat/lon, everything else is in this
if boundingbox: stops = stops.cx[boundingbox[2]:boundingbox[3], boundingbox[0]:boundingbox[1]]
for i, row in stops.iterrows():
marker = folium.Marker(
location = [i[0] for i in reversed(row.geometry.coords.xy)],
icon = folium.Icon(color='blue', prefix='fa', icon='bus'),
popup="""
<b>{}</b> <br />Bus stop: {} <br />Routes: {}""".format(row.PlatformNa, int(row.PlatformNo) if pd.notnull(row.PlatformNo) else '-', row.RouteNos)
if city_name == 'Christchurch' else
"""<b>{}</b> <br />Bus stop: {}""".format(row.Service, row.Stop if pd.notnull(row.Stop) else '-')
)
stops.at[i, 'marker'] = marker
return (routes, stops)
# parse command line
p = OptionParser(usage="""usage: %prog [options] <city>
Create map with school zones for given city
Currently can only do Christchurch or Wellington
""")
p.add_option("-v", action="store_true", dest="verbose", help="Verbose")
(opts, args) = p.parse_args()
if len(args) != 1:
p.print_help()
sys.exit(1)
city_name = args.pop(0)
if city_name not in ('Christchurch', 'Wellington'):
p.print_help()
sys.exit(1)
# find location
loca = geopy.geocoders.Nominatim(user_agent='MyChChAp')
city = loca.geocode('%s, New Zealand' % city_name)
# make an interactive map of Christchurch
m = folium.Map(location=(city.latitude, city.longitude), zoom_start=13, control_scale=True)
######################################################################
# get the schools we are interested in
schools = get_schools([float(i) for i in city.raw['boundingbox']])
# add school markers to correct group
prim_group = folium.FeatureGroup("Full primary and intermediate schools")
high_group = folium.FeatureGroup("High schools")
for i, row in schools.iterrows():
(high_group if row.high else prim_group).add_child(row.marker)
# add the zonepoly's to the map
for i, row in schools.iterrows():
if pd.notnull(row.zonepoly): row.zonepoly.add_to(m)
######################################################################
######################################################################
# metro info
routes, stops = get_metro(city_name, [float(i) for i in city.raw['boundingbox']])
routes_group = folium.FeatureGroup('Bus routes')
for i, row in routes.iterrows():
routes_group.add_child(row.polyline)
stops_group = folium.FeatureGroup('Bus stops', show=False)
for i, row in stops.iterrows():
stops_group.add_child(row.marker)
######################################################################
# add the layers
m.add_child(prim_group)
m.add_child(high_group)
m.add_child(routes_group)
m.add_child(stops_group)
m.add_child(folium.plugins.MeasureControl(primary_length_unit='kilometers', primary_area_unit='hectares'))
folium.LayerControl().add_to(m)
######################################################################
# add the custom leaflet javascript code so when a marker is clicked
# the zonepoly gets shown
#
# the method followed is from FabeG on
# https://github.com/python-visualization/folium/issues/86
# (see near the end posted around the end of April 2018)
#
from branca.element import Element
e = []
# when zoomed in enough make the bus stops layer visible
e.append("""{0}.on('zoomend', function() {{
if ({0}.getZoom() < 17 ) {{
{0}.removeLayer({1});
}}
else {{
{0}.addLayer({1});
}}
}});""".format(m.get_name(), stops_group.get_name()))
# the zonepoly's start as visible, remove them all
for i, row in schools.iterrows():
if pd.notnull(row.zone): e.append("""{0}.remove();""".format(row.zonepoly.get_name()))
#
# using semiomant's suggestion on
# https://stackoverflow.com/questions/13698975/click-link-inside-leaflet-popup-and-do-javascript
# this makes a function that is called on every popupopen and it
# checks to make sure that their is a callback bound to every zonelink
#
# to stop the callback being bound repeatedly, use Kondra Garus's
# suggestion on
# https://stackoverflow.com/questions/6361465/how-to-check-if-click-event-is-already-bound-jquery
# which is the :not(.bound) thing
#
e.append("""{0}.on('popupopen', function() {{""".format(m.get_name()))
for i, row in schools.iterrows():
if pd.isnull(row.zone): continue
e.append("""
$("#zonelink{1}:not(.bound)").addClass('bound').on('click', function(event) {{
if({0}.hasLayer({2})) {{
$(this).text('Show');
{0}.removeLayer({2});
}} else {{
$(this).text('Hide');
{0}.addLayer({2});
}}
}}
)""".format(m.get_name(), row.School_Id, row.zonepoly.get_name()))
e.append("});")
# put the new JS at the end
e = Element("\n".join(e))
html = m.get_root()
html.script.get_root().render()
html.script._children[e.get_name()] = e
######################################################################
# write it out
m.save('docs/%s.html' % city_name)