import geopandas as gp
import folium
from shapely.geometry import Point
from geopy.geocoders import GoogleV3
from IPython.display import HTML
The first step is generating a base layer upon which we intend to overlay features.
# help(folium.Map)
geolocator=GoogleV3()
place,(latitude,longitude)=geolocator.geocode('Washington DC')
print place, latitude, longitude
Washington, DC, USA 38.9072309 -77.0364641
Word up, easy as pie. The next step is enabling the rendering of interactive maps within the Notebook. This is a really cool feature that ggmap cannot match. It does, however, involve initiating a local server. This is a trivial exercise, just type the following in the shell:
python -m SimpleHTTPServer 8015
The number is important here, because it is searched for directly (as seen in the iframe tag string below. Once this is going, we can define a couple convenience functions.
#Function that shows the map
def showmap(map_in,pth):
map_in.create_map(path=pth)
return HTML('<iframe src=http://localhost:8015/'+pth+' width=1000 height=500></iframe>')
#Function that resets the map
def resetmap():
return folium.Map(location=[40, -99], zoom_start=4)
And now we are ready to plot our interactive tile. (It should be noted that static HTML renderings of this Notebook are ... well, static. If interactivity is desired, it must be done within an IPython Notebook.)
dcmap=folium.Map(location=[latitude,longitude],tiles='Stamen Toner')
dcmap.create_map(path='index.html')
showmap(dcmap,'index.html')
Now that we have a handle on the base layer, let's see how we can include some features of interest. It is possible to build in individual markers, but who wants to input coordinates when many markers must be added? Given that there are many cameras, the tedium would be overwhelming (and make for incredibly boring code).
It turns out that folium supports GeoJSON overlays.
#Establish working directory
workdir='/home/choct155/ORA/misc/'
#Read in camera shp
camera=gp.GeoDataFrame.from_file(workdir+'TrafficCamera.shp')
print camera
camera.head()
<class 'geopandas.geodataframe.GeoDataFrame'> Int64Index: 339 entries, 0 to 338 Data columns (total 8 columns): CAMERAID 339 non-null values CAMERATYPE 339 non-null values FACILITYID 339 non-null values OWNER 339 non-null values POLEID 339 non-null values STREETJUNC 339 non-null values STREETSEGI 339 non-null values geometry 339 non-null values dtypes: int64(6), object(2)
CAMERAID | CAMERATYPE | FACILITYID | OWNER | POLEID | STREETJUNC | STREETSEGI | geometry | |
---|---|---|---|---|---|---|---|---|
0 | 329 | 1 | 329 | 1 | 25657 | 7507 | 2760 | POINT (400785.3803125023841858 129044.05624999... |
1 | 330 | 1 | 330 | 3 | 25673 | 2970 | 13183 | POINT (401033.8136250004172325 129436.28687499... |
2 | 332 | 1 | 332 | 3 | 26542 | 7467 | 13152 | POINT (402458.9065624997019768 132123.35637500... |
3 | 331 | 2 | 331 | 3 | 26425 | 4075 | 2895 | POINT (401478.2784375026822090 131751.35325000... |
4 | 333 | 2 | 333 | 3 | 27369 | 2813 | 2214 | POINT (402577.9644374996423721 133946.26999999... |
Sweet Jesus! There are 339 camera, which definitely means we need to leverage this overlay business. In any event, the coordinate reference system (CRS) for the camera data is currently in eastings/northings. If we actually want the cameras to appear anywhere in the vicinity of the District, it's probably in our interest to harmonize the CRSs.
camera_ll=camera.to_crs(crs={'proj':'longlat',
'datum':'WGS84'}).rename(columns={'geometry':'geo_temp'})
camera_ll.head()
CAMERAID | CAMERATYPE | FACILITYID | OWNER | POLEID | STREETJUNC | STREETSEGI | geo_temp | |
---|---|---|---|---|---|---|---|---|
0 | 329 | 1 | 329 | 1 | 25657 | 7507 | 2760 | POINT (-76.9909549637256276 38.8291823043240072) |
1 | 330 | 1 | 330 | 3 | 25673 | 2970 | 13183 | POINT (-76.9880932278285570 38.8327154426776175) |
2 | 332 | 1 | 332 | 3 | 26542 | 7467 | 13152 | POINT (-76.9716703595549774 38.8569189260698664) |
3 | 331 | 2 | 331 | 3 | 26425 | 4075 | 2895 | POINT (-76.9829692059042543 38.8535699575756723) |
4 | 333 | 2 | 333 | 3 | 27369 | 2813 | 2214 | POINT (-76.9702918311025286 38.8733401431466348) |
ESRI has created an extra hurdle for us. For whatever reason, they like to kick out coordinates in reverse order. It's challenging to plot features with reversed coordinates in the right place, so let's see what can be done about it. (This was the reason for renaming the 'geometry' variable above.)
Swapping these coordinates teaches us just a little about the geopandas machinery. The spatial information is actually captured with Shapely objects. For our camera point layer, the objects in the 'geometry' column are ... Point objects. A Point object has an interior set of one point, no points in the boundary set, and an exterior set of all other points. Just a little trivia...
In any case, it's important to know that Shapely's geometric objects are immutable, so we must make a copy that switches the (basically) tuple position of the coordinates from long-lat to lat-long. This is a tidy job for a list comprehension.
#Swap lat-long coords in Point objects
camera_ll['geometry']=[Point(x.y,x.x) for x in camera_ll['geo_temp'].values]
#Keep relevant columns
camera_ll2=camera_ll[['CAMERAID','CAMERATYPE','FACILITYID','OWNER','POLEID','STREETJUNC','STREETSEGI','geometry']]
camera_ll2.head()
CAMERAID | CAMERATYPE | FACILITYID | OWNER | POLEID | STREETJUNC | STREETSEGI | geometry | |
---|---|---|---|---|---|---|---|---|
0 | 329 | 1 | 329 | 1 | 25657 | 7507 | 2760 | POINT (38.8291823043240072 -76.9909549637256276) |
1 | 330 | 1 | 330 | 3 | 25673 | 2970 | 13183 | POINT (38.8327154426776175 -76.9880932278285570) |
2 | 332 | 1 | 332 | 3 | 26542 | 7467 | 13152 | POINT (38.8569189260698664 -76.9716703595549774) |
3 | 331 | 2 | 331 | 3 | 26425 | 4075 | 2895 | POINT (38.8535699575756723 -76.9829692059042543) |
4 | 333 | 2 | 333 | 3 | 27369 | 2813 | 2214 | POINT (38.8733401431466348 -76.9702918311025286) |
#Open file for writing
f=open(workdir+'camera_ll2.json','w')
#Write the GeoJSON info to disk
f.write(camera_ll2.to_json())
#Close file
f.close()
dcmap.geo_json(geo_path=workdir+'camera_ll2.json')
showmap(dcmap,'index.html')
Hmmm, definitely having some trouble here. My conversion from NAD83 to WGS84 seems to have been appropriate. I am thinking the real issue is that we are dealing with an object that has a topological dimension of zero. To test this theory, I am going to manually add a marker.
# dcmap.simple_marker(location=[38.8733401431466348,-76.9702918311025286],popup='Aesthetic Test - Camera 333')
# showmap(dcmap,'dcmap.html')
My suspicions are confirmed here. I need to somehow add markers for all of the points. In truth, I probably overstated how tedious this would be. We can cycle through the GeoDataFrame and extract the relevant information. First, however, we need to identify a less obtrusive marker than this Trulia-type nonsense.
dcmap.polygon_marker(location=[38.8733401431466348,-76.9702918311025286],popup='Test - Camera #333',
fill_color='#B80000',num_sides=4,radius=6)
showmap(dcmap,'index.html')
Now we are in good shape. We can roll though each of these points and populate the map with all cameras in the District. We need to extract three pieces of information:
This is clearly a loop job that pulls from two locations (columns) in the GDF. For pedagogical reasons, it's useful to highlight how I extract these values from the GDF. The first two (latitude and longitude) are extracted from the objects captured in the 'geometry' column. The last one (camera ID #) is captured from the column of the same name.
print camera_ll2['geometry'].values[0].x
print camera_ll2['geometry'].values[0].y
print camera_ll2['CAMERAID'].values[0]
38.8291823043 -76.9909549637 329
Now, let us evaluate all locations in the GDF.
#For every camera in the GDF...
for i in range(len(camera_ll2['geometry'])):
#...extract the coordinates...
lat=camera_ll2['geometry'].values[i].x
lon=camera_ll2['geometry'].values[i].y
#...and the pop up text...
text=str(camera_ll2['CAMERAID'].values[i])
#...and create a marker
dcmap.polygon_marker(location=[lat,lon],popup=text,
fill_color='#B80000',num_sides=4,radius=6)
Presumably we just made markers en masse. Let's see if we did what we thought we did.
showmap(dcmap,'index.html')