Lesson 05 – Simple Maps

---

Environmental data almost always comes with latitude and longitude information: “Where was this data collected?” The most natural thing in the world is to display locations on a map. The ‘maps’ package is very simple to use and an excellent choice if your data is located in North America. The biggest downside of the ‘maps’ package is the 1980’s, Soviet Union era borders used in eastern Europe.

Note: The ‘sp’ package gives R an impressive array of GIS capabilities and can be used to read in, analyze and plot shapefile data. The ‘sp’ package is recommended if you intend to do serious geospatial work.

The ‘maps’ package includes ‘world’, ‘usa’, ‘state’ and ‘county’ databases with named elements. To figure out the naming scheme you call the map() function and pass in the ‘namesonly’ argument:

library(stringr)   # for string manipulation (str_~ functions)
library(maps)      # for maps

regionNames <- map('world', plot=FALSE, namesonly=TRUE)
# NOTE:  length(regionNames) # 2284 names includes lots  "country:island" names

# Here is a sorted, unique vector of region names up to the first ':' (see ?str_match, ?regex)
countryNames <- sort(unique(str_match(regionNames,'[^:]+')))

# The US states are in the 'state' database
stateNames <- map('state', plot=FALSE, namesonly=TRUE)

# NOTE:  The maps() function uses the default argument 'exact=FALSE' which forces regular
# NOTE:  expression matching against lower case state names in the database.  Thus we 
# NOTE:  can reference Washington state as 'washington', 'WA' or 'wa' or even 'wAsH'.
# NOTE:  But 'NE' will match Nebraska, Nevada and all of the 'New ...' states.

map('usa',col='gray50',lwd=0.5)
map('state','NE',add=TRUE)

You can add point locations to any map by simply calling points(). Here is an example using the fire locations dataset we have used in previous lessons.

url <- "http://smoke.airfire.org/bluesky-daily/output/hysplit-pp/NAM-4km/2014080100/data/fire_locations.csv"
fires <- read.csv(url, stringsAsFactors=FALSE)

# Create a map for the region of interest
map('world',regions=c('canada','usa','mexico','cuba','hawaii','guatemala','honduras','belize','el salvador','nicaragua'),
    xlim=range(fires$longitude),ylim=range(fires$latitude)) # see ?range

# Overlay state boundaries
map('state', add=TRUE, col='gray50', lwd=0.5)

# Plot the fires on top
points(fires$longitude, fires$latitude)

Not bad but I’m sure we can do better.

One of the great things about arguments to the points() function is that they are all vectorized. This means that you can use other values from the fires dataframe to set symbol, size and color. We can also take advantage of the fact that factors are converted to integer where appropriate.

Here is a map with:

• symbol determined by fire ‘type’
• color determined by ‘heat’
• symbol size determined by log2(‘area’)

oldPar <- par()
par(mar=c(0,0,3,0))  # no margins except at the top

# NOTE:  as.factor(fires$type) will return a vector of 1's and 2's which will be used as indices
pchs <- c(15,17)[as.factor(fires$type)]

# NOTE:  set up lo-med-hi colors and then use .bincode() and quantile() to assign bin #s. 
# NOTE:  By default, NA values are assigned 'transparent'.  We will assign them a medium gray.
loMedHiIndex <- .bincode( fires$heat, quantile(fires$heat,c(0,.33,.67,1.0),na.rm=TRUE) )
cols <- c('orange','red','firebrick')[loMedHiIndex]
cols[is.na(fires$heat)] <- 'gray50'
cols <- adjustcolor(cols,0.1)  # 10% opacity

# Symbol size
cexs <- log2(fires$area) / 8

# Create a map for the region of interest
map('world',
    regions=c('canada','usa','mexico','cuba','hawaii','guatemala','honduras','belize','el salvador','nicaragua'),
    xlim=range(fires$longitude),ylim=range(fires$latitude)) # see ?range

# Overlay state boundaries
map('state',
    add=TRUE, col='gray50', lwd=0.5)

# Plot the fires on top
points(fires$longitude, fires$latitude,
       pch=pchs, col=cols, cex=cexs)

# Add a box
box(lwd=1,lty='solid',col='gray50')

# Add a title and legend
title('Fire Locations on August 1, 2014')
legend('left',
       title="Symbols Used",
       legend=c('WF','RX','warm','hot','HOT','small','med','large'),
       col=adjustcolor(c('gray50','gray50','orange','red','firebrick','gray50','gray50','gray50'),0.7),
       pch=c(15,17,16,16,16,16,16,16),
       pt.cex=c(1,1,1,1,1,1,2,3))

# Restore old graphical parameters
par(oldPar)

---

Task 1: Western Wildfires

• use dplyr’s filter() function to create a wf dataframe
• plot these fires on a map to see what states they are in
• create a vector of state names defining western states
• plot these states with black borders
• add counties in a lighter color
• add wildfires with meaningful col/pch/cex

---

---

Task 2: Projected USA map

The ‘mapproj’ package allows you to perform projections. Try the following:

• use RStudio to install the ‘mapproj’ package
• library(mapproj) and ?mapproject to read up just a little
• use map() with “database=‘state’, proj=‘albers’,par=c(30,40)”
• use points(mapproject(fires$longitude,fires$latitude), …)

---