After finding a street tree map for NYC and an “eco impact” map of San Francisco, I decided to make one for the City of Vancouver. The city has released government data online for public access and use. It’s fantastic to see that the city has worked towards providing these resources for citizen engagement.
I joined up with a good friend Laura Jickling to create a website infographic on tumblr. On it you’ll find the interactive map. We decided we start with a neighborhood as a pilot and visualize a smaller batch of the data set.
Screenshot of map:
View the code below:
######################################################
Mount Pleasant Street Trees Data
Date: Nov. 2015
By: Laura Jickling & Audrey Schalhoub
Desc: This script geocodes and scatters street tree data for Vancouver's Mount Pleasant neighbourhood. The data was acquired from the Vancouver Open Data Catalogue.
######################################################
-------------------------------------------------
--------------- Install Libararies ---------------
-------------------------------------------------
install.packages("GISTools") install.packages('rgdal') install.packages('curl') install.packages("devtools") devtools::install_github("corynissen/geocodeHERE")
-------------------------------------------------
---------------- Load Libararies -----------------
-------------------------------------------------
library(GISTools) library(rgdal) library(geocodeHERE) library(curl)
-------------------------------------------------
------------------- Acquire ---------------------
-------------------------------------------------
access the street tree data from local computer
fname = 'C:/Users/laura/Desktop/trees/csv_street_trees/StreetTrees_MountPleasant.csv'
Read data as csv
data = read.csv(fname, header=T)
inspect your data
print(head(data))
-------------------------------------------------
-------------- Parse: Geocoder -------------------
-------------------------------------------------
create a new column with full address for block that tree is facing:
data$block = paste(data$ON_STREET_BLOCK, data$ON_STREET, sep=" ") data$tree_block = paste(data$block, "Vancouver, BC, Canada", sep=" , ")
write out the file
ofile ='C:/Users/laura/Desktop/trees/mp_trees5.csv' write.csv(data, ofile)
Geocode the events using Nokia's Here API
Create an empty vector for lat and lon coordinates
lat = c(); lon = c()
loop through the addresses
for(i in 1: length(data$tree_block)){
store the address at index "i" as a character
address = as.character(data$tree_block[i])
append the latitude of the geocoded address to the lat vector
lat = c(lat,geocodeHERE_simple(address, App_id = "FyE0xI0xlMfHknNvSuAK", App_code="R9v4al_ox3ND2pDKLiGHzA")$Latitude)
append the longitude of the geocoded address to the lon vector
lon = c(lon,geocodeHERE_simple(address, App_id = "FyE0xI0xlMfHknNvSuAK", App_code="R9v4al_ox3ND2pDKLiGHzA")$Longitude) }
add the lat lon coordinates to the dataframe
data$lat = lat data$lon = lon
write out the file
ofile ='C:/Users/laura/Desktop/trees/mp_trees_block_geo.csv' write.csv(data, ofile)
access the street tree data from local computer
fname = 'C:/Users/laura/Desktop/trees/mp_trees_block_geo.csv'
Read data as csv
data = read.csv(fname, header=T)
#Plot the data and look for outliers plot(data$lon, data$lat)
#filter out the outliers block_filter = subset(data, (lat <= 49.272) & (lat >= 49.25) & (lon <= -123.076) & (lon >= -123.12) & is.na(lon) == FALSE )
#Plot the data to see if there are still outliers plot(block_filter$lon, block_filter$lat)
#----------------------------------------------------------------------#
#-----------------------Filter further---------------------------------#
--- Convert Data to Shapefile & GeoJson to use in QGIS or ArcMap ---
#----------------------------------------------------------------------#
store coordinates in dataframe
coords_crime = data.frame(block_filter$lon_offset, block_filter$lat_offset)
create spatialPointsDataFrame
data_shp = SpatialPointsDataFrame(coords = coords_crime, data = block_filter)
set the projection to wgs84
projection_wgs84 = CRS("+proj=longlat +datum=WGS84") proj4string(data_shp) = projection_wgs84
set an output folder for the spatial files:
geofolder2 = 'C:/Users/laura/Desktop/trees/geofolder2/'
Join the folderpath to each file name:
opointsshp = paste(geofolder2, "trees","mp_blocks_test",".shp", sep="") print(opoints_shp) # print to see file name and location opointsgeojson = paste(geofolder2, "trees","mp_blocks_test",".geojson", sep="") print(opoints_geojson)# print to see file name and location
write the file to a shp
writeOGR(check_exists = FALSE, data_shp, opoints_shp, layer="data_shp", driver="ESRI Shapefile")
write file to geojson
writeOGR(check_exists = FALSE, data_shp, opoints_geojson, layer="data_shp", driver="GeoJSON")
-------------------------------------------------
--------------------- Mine ----------------------
-------------------------------------------------
store the file name of the tree shapefile we filtered on ArcMap
mp_filter2 = 'C:/Users/laura/Desktop/trees/geofolder2/trees_blocks_take2.shp'
read in the data
filter= readOGR(mp_filter2, layer="trees_blocks_take2")
#convert shapefile to data frame filter2.df = as(filter, "data.frame")
Classify: Group the tree types in Classes by Genus
get a frequency distribution of the genuses:
top_trees = data.frame(table(filter2.df$GENUS_N)) top_trees = top_trees[order(top_trees$Freq),] print(top_trees)
#Create groups for the top 10 Genuses and assign common name
#Maple Maple='ACER'
#Cherry or Plum Cherry_or_Plum='PRUNUS'
#Pear Pear='PYRUS'
#Oak Oak='QUERCUS'
#Linden Linden='TILIA'
#cHESNUT Chesnut='AESCULUS'
#Hornbeam Hornbeam= 'CARPINUS'
#Beech Beech= 'FAGUS'
#Ash Ash='FRAXINUS'
#Other Other=c('LIQUIDAMBAR','CRATAEGUS','ULMUS','CERCIDIPHYLLUM','HIBISCUS','SORBUS', 'MALUS','MAGNOLIA','SYRINGA','ROBINIA','BETULA','STYRAX','PARROTIA','CORNUS','GLEDITSIA','CHAMAECYPARIS','STEWARTIA','CERCIS','PLANTANUS','THUJA','LIRIODENDRON','GINKO','NYSSA','METASEQUOIA','CLADRASTIS','DAVIDIA','RHUS','PINUS','PICEA','ILEX','PSEUDOTSUGA','OSTRYIA','ABIES','SEQUOIADENDRON','CEDRUS','LABURNUM','POPULUS','ALNUS','ZELKOVA','TSUGA', 'SOPHORA','JUGLANS','CATALPA','CASTANEA','GINKGO','PLATANUS')
give class id numbers to each group:
filter2.df$cid = 9999 for(i in 1:length(filter.df$GENUS_N)){ if(filter2.df$GENUS_N[i] %in% Maple){ filter2.df$cid[i] = 0 }else if(filter2.df$GENUS_N[i] %in% Cherry_or_Plum){ filter2.df$cid[i] = 1 }else if(filter2.df$GENUS_N[i] %in% Pear){ filter2.df$cid[i] = 2 }else if(filter2.df$GENUS_N[i] %in% Oak){ filter2.df$cid[i] = 3 }else if(filter2.df$GENUS_N[i] %in% Linden){ filter2.df$cid[i] = 4 }else if(filter2.df$GENUS_N[i] %in% Chesnut){ filter2.df$cid[i] = 5 }else if(filter2.df$GENUS_N[i] %in% Hornbeam){ filter2.df$cid[i] = 6 }else if(filter2.df$GENUS_N[i] %in% Beech){ filter2.df$cid[i] = 7 }else if(filter2.df$GENUS_N[i] %in% Ash){ filter2.df$cid[i] = 8 }else if(filter2.df$GENUS_N[i] %in% Other){ filter2.df$cid[i] = 9 } }
#----------------------------------------------------------------#
#------------------------------Jitter----------------------------#
--- handle overlapping points and allign them with streets ---
#----------------------------------------------------------------#
offset points along the blockface:
#Create group for East Roads east=c('E 1ST AV', 'E 2ND AV', 'GREAT NORTHERN WAY', 'ATHLETES WAY', 'WALTER HARDWICK AV', 'E 3RD AV', 'E 4TH AV', 'E 5TH AV', 'E 6TH AV', 'E 7TH AV', 'E 8TH AV', 'E BROADWAY', 'E 10TH AV', 'E 11TH AV', 'E 12TH AV', 'E 13TH AV', 'E 14TH AV', 'E 15TH AV', 'E 16TH AV', 'SWITCHMEN ST', 'CENTRAL ST','TERMINAL AV')
#create group for west roads west= c('W 1ST AV','W 2ND AV', 'W 3RD AV', 'W 4TH AV', 'W 5TH AV', 'W 6TH AV', 'W 7TH AV', 'W 8TH AV', 'W BROADWAY', 'W 10TH AV', 'W 11TH AV', 'W 12TH AV', 'W 13TH AV', 'W 14TH AV', 'W 15TH AV', 'W 16TH AV')
#Create group for Notrh-South Roads north_south=c('CAMBIE ST', 'YUKON ST', 'ALBERTA ST', 'COLUMBIA ST', 'MANITOBA ST', 'SALT ST', 'ONTARIO ST', 'QUEBEC ST', 'MAIN ST', 'WATSON ST', 'SOPHIA ST', 'PRINCE EDWARD ST', 'ST. GEORGE ST', 'GUELPH ST', 'CAROLINA ST', 'BRUNSWICK ST', 'SCOTIA ST', 'PRINCE ALBERT ST', 'FRASER ST', 'ST. CATHERINES ST', 'WINDSOR ST', 'GLEN DRIVE', 'CLARK DRIVE', 'KEITH DRIVE')
#Create group for irregular streets other_streets=c('WYLIE ST', 'COOK ST', 'CROWE ST', 'LORNE ST', 'THORNTON ST', 'KINGSWAY')
#Offset lon and lat for each group filter2.df$lat_offset2 = filter2.df$lat filter2.df$lon_offset2 = filter2.df$lon
for(i in 1:length(filter2.df$lat)){ if (filter2.df$ON_STREET[i] %in% east){ filter2.df$lon_offset2[i] = filter2.df$lon[i] + runif(1,0.0008,0.002) filter2.df$lat_offset2[i] = filter2.df$lat[i] + 0.0001 }else if(filter2.df$ON_STREET[i] %in% west){ filter2.df$lon_offset2[i] = filter2.df$lon[i] - runif(1,0.0008,0.002) filter2.df$lat_offset2[i] = filter2.df$lat[i] + 0.0001 }else if(filter2.df$ON_STREET[i] %in% north_south){ filter2.df$lon_offset2[i] = filter2.df$lon[i] - 0.0001 filter2.df$lat_offset2[i] = filter2.df$lat[i] - runif(1,0.-0.0004,0.0004) }else if(filter2.df$ON_STREET[i] %in% other_streets){ filter2.df$lon_offset2[i] = filter2.df$lon[i] + runif(1,0.0001,0.0003) filter2.df$lat_offset2[i] = filter2.df$lat[i] + runif(1,0.0001,0.0003) } }
#Save the file as a geojson to use in CartoDB#
store coordinates in dataframe
coords_trees4 = data.frame(filter2.df$lon_offset2, filter2.df$lat_offset2)
create spatialPointsDataFrame
data_shp4 = SpatialPointsDataFrame(coords = coords_trees4, data = filter2.df)
set the projection to wgs84
projection_wgs84 = CRS("+proj=longlat +datum=WGS84") proj4string(data_shp4) = projection_wgs84
Join the folderpath to each file name:
opoints_geojson = paste(geofolder2, "trees","blocks_please33",".geojson", sep="") print(opoints_geojson)# print to see file name and location
write file to geojson
writeOGR(check_exists = FALSE, data_shp4, opoints_geojson, layer="data_shp4", driver="GeoJSON")