File with complete R code can be found here: Full R Code
We’ll get started by first downloading leaflet, and installing a few packages. Leaflet for R is still under development and is not yet available on CRAN, so we download from Github.
We’ll also need the magrittr package to use the forward pipe operator (%>%) and maptools to import and manage shapefiles, which we’ll get to later.
# First, install leaflet
require(devtools)
devtools::install_github("rstudio/leaflet")
# then install a few extra packages
library(leaflet)
library(magrittr)
library(maptools)We’ll be working with geocoded data from the EPA on pollution and particulate matter from major point sources. Data was available for csv download by state, so I’ve focused here on the point sources for NOx in New York state.
Here’s the EPA page, if you’d like to try this yourself: EPA Emissions Sources Data
First thing first, we need to read in our data and inspect the file structure.
## read in data as csv file
NOx_NY <- read.csv(path.expand("~/Desktop/NOx_NY.csv"), header=TRUE, sep=",")
## review data format
head(NOx_NY)## Facility.Name
## 1 BORALEX NEW YORK LP
## 2 PEARL LEATHER FINISHERS INC
## 3 CEDAR CREEK WPCP
## 4 TRIGEN CENTRL UTILITY PLT - MITCHL FIELD
## 5 GLENWOOD LANDING ENERGY CENTER
## 6 GLOBE METALLURGICAL INC
## Facility.Address Latitude Longitude NAICS.Code
## 1 7019 ST RTE 374 44.8935 -74.0722 221112
## 2 11-21 INDUSTRIAL PARK 43.0275 -74.3539 31611
## 3 3340 MERRICK RD & CEDAR CREEK PK 40.6523 -73.5065 22132
## 4 185 CHARLES LINDBERGH BLVD 40.7260 -73.5897 22121
## 5 SHORE RD 40.8320 -73.6466 221112
## 6 3807 HIGHLAND AVE 43.1229 -79.0418 331410
## Pollutant Source.Type Emissions.in.Tons
## 1 Nitrogen Oxides Electricity Generation via Combustion 131.660855
## 2 Nitrogen Oxides 0.866800
## 3 Nitrogen Oxides Wastewater Treatment Facility 146.327375
## 4 Nitrogen Oxides Electricity Generation via Combustion 271.463614
## 5 Nitrogen Oxides Electricity Generation via Combustion 0.804625
## 6 Nitrogen Oxides 410.124500
## Year
## 1 2011
## 2 2011
## 3 2011
## 4 2011
## 5 2011
## 6 2011Great! We’ll put this aside for use once we review the basic leaflet functions.
The basic leaflet() function can be called with a host of parameters inside it, but using the pipe operator makes this much simpler so we’ll use that method here.
The first thing you’ll need is the addTiles() function. This pulls in a set of map tiles from the internet and forms the basis of your interactive map.
## basic leaflet map
leaflet() %>% addTiles()Note that the default is to reference the OpenStreetMaps tiles, but other sets are available. A list of can be found here:
Let’s pick a more interesting background…
## map attributions also recommended
m = leaflet() %>% addTiles('http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
attribution = 'Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid,IGN, IGP, UPR-EGP, and the GIS User Community')
mThere we go. Note that this code pulls in a tile set, and allows you to set attributions for the tiles in the bottom right. Always be sure to included attribution for whoever made the map tiles!
This code also assigns the leaflet map to “m”. We can now use the pipe operator to put different layers on this basic map. This is useful if you want to manage several different map backgrounds, or a single default map that you want to show with different data layers.
Though the map has a zoom button built in, you can also set the initial starting map center and zoom level explicitly.
m %>% setView(-77.03653,38.897676, zoom = 16)Ok, let’s layer on that EPA data from earlier in the simplest way possible.
m %>% addCircles(data = NOx_NY)## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectivelyCool, huh? Leaflet automatically searches for columns in your dataframe named latitude/longitude, or some variant of those words and plots to those points. The zoom is set to the closest level that contains your entire set of points.
You also have some flexibility with this call. Here we explicitly override any data already in m with %>% to plot our NOx sites. If we already had that data referenced in m, we could just do this, with the same result:
m = leaflet(data = NOx_NY) %>% addTiles() %>% addCircles()
mLet’s do something more interesting and weight those dots by emissions amount.
types <- as.list(as.character(unique(NOx_NY$Source.Type)))
m %>% addTiles() %>% addCircleMarkers(data = subset(NOx_NY, Source.Type == "Airport"),
radius = ~Emissions.in.Tons/80, color = 'red')## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectivelyThis did a few things. We are only now displaying a subset of the original data (Airports) and have weighted the radius by a normalized size. We also can update color, and a host of other attributes.
You can also drop customized icons onto your map by geocode. Say we wanted to highlight the biggest statewide NOx polluter. With a quick calculation and a reference to the icon we want to use…
## Layering Markers onto the map
rowmax <- NOx_NY[which.max(NOx_NY$Emissions.in.Tons),3:4]
m %>% addCircleMarkers(data = NOx_NY,
radius = ~Emissions.in.Tons/80, color = 'red') %>% addMarkers(rowmax[1,2], rowmax[1,1],
icon = JS("L.icon({iconUrl: 'http://www.grinningplanet.com/2004/04-20/industry-camera-thumb.gif', iconSize: [30,30]})"))## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectivelyYou can graph more than just circles, too. Leaflet supports popup markers, lines, area enclosing rectangles, and polygons as well. Here are examples of popups and connecting lines.
## Polylines example, with a line connecting each electricity generation site
linetest <- subset(NOx_NY, Source.Type == "Electricity Generation via Combustion")
m %>% addPolylines(linetest$Longitude, linetest$Latitude, stroke=0.2) %>% addCircles(data=linetest, color="red")## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectively## Example adding a Popup marker
m %>% addPopups(-73.959902,40.809379, "Here's our classroom!")If you want to highlight a region, map a specific route, or do anything beyond simple lines, circles, and rectangles you’ll want to use shapefiles. Everything we’ve done up until now has read long/lat data directly from a data frame. Using maptools (or a handful of other packages), you can also upload shapefiles as a SpatialPointsDataFrame and map more complex polygons onto your base map.
To work through this, I downloaded a shapefile of NY county boundaries: NY Country Shapefile
From here, we can use the maptools package to read in the shapefile data. For what it’s worth, the rgdal package is supposed to be simpler to manage but a binary for OSX Mavericks wasn’t available on CRAN…
## import shapefile for county lines and add projection class (maptools forces you to do this manually)
countylines <- readShapeSpatial("~/Desktop/NYCountyShp/CEN2000nycty.036.shp.07865/cty036.shp")
countylines$proj4string <- "+proj=longlat +ellps=clrk66"After reading in the data and manually setting the projection coordinate system, we now can apply these polygons to the leaflet map.
## add counties as polygons
m %>% addPolygons(data=countylines, fill = TRUE, color = "blue")These can again be tuned to different colors, fill opacities, and border sizes.
m %>% addPolygons(data=countylines,
fill = TRUE,
color = "yellow",
fillOpacity= 0.4,
stroke = FALSE)Now that we know all of the pieces we want to include in our map, bringing them all togetehr is simple with the %>% operator.
m = m %>% addPolygons(data=countylines, fill = TRUE, color = "green", fillOpacity= 0.1, stroke = FALSE)
m = m %>% addPolylines(data=countylines, weight = 2, color = "white")
m = m %>% addCircleMarkers(data = subset(NOx_NY, Source.Type == "Airport"),
radius = ~Emissions.in.Tons/80,
color = 'red')## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectivelym = m %>% addCircleMarkers(data = subset(NOx_NY, Source.Type == "Electricity Generation via Combustion"),
radius = ~Emissions.in.Tons/80,
color = 'blue')## Assuming 'Longitude' and 'Latitude' are longitude and latitude, respectivelym = m %>% addMarkers(rowmax[1,2], rowmax[1,1],
icon = JS("L.icon({iconUrl: 'http://www.grinningplanet.com/2004/04-20/industry-camera-thumb.gif', iconSize: [30,30]})"))
m