R Spatial Vignette

This is an R vignette to introduce spatial data analysis. Spatial data can be stored as and comes in many formats. The first part of the vignette will introduce how spatial data can be visualized in web-based platforms through Google Visualisation API, the use of basemaps, selecting areas, and plotting spatial data into a web map.

The next part of the vignette will give an overview of how to use R to load geographic data of vector (1) points and (2) polygons into a spatial dataframe, then analyze in R and a webmap.

Section 1: Create a map from R for web-based platforms

R Spatial packages for introduction to Basemaps and Webmaps

library(sp)  # classes for spatial data
suppressPackageStartupMessages(library(rgeos)) # needed for maptools
suppressPackageStartupMessages(library(maptools)) # has spatial data
library(raster) #needed for dismo
library(dismo) # retrieving base maps from Google with gmap function
library(RgoogleMaps) # web based map abilities through Google server
suppressPackageStartupMessages(library(googleVis)) # interactive web based maps with data frames

## Creating a generic function for 'toJSON' from package 'jsonlite' in package 'googleVis'

Create a basemap in R; zoom/select areas and plot.

# BasemapJapan <- gmap("Japan")
# plot(BasemapJapan)

Change the basemap, options: ‘roadmap’, ‘satellite’, ‘hybrid’, ‘terrain’

# SatJapan <- gmap("Japan", type = "satellite", exp = 1, filename = "Japan.gmap")
# plot(SatJapan)

Manually choose a spatial area of interest:

# select.area <- drawExtent()
# SelectArea <- gmap(select.area)
# plot(SelectArea)

Visualise data in a web browser using Google Visualisation API.

Note: gvisGeoMap needs Chrome set as default browser, example from R Documentation.

# WorldPopulation=data.frame(Country=Population$Country,
#         Population.in.millions=round(Population$Population/1e6,0),
#         Rank=paste(Population$Country, "Rank:", Population$Rank))
# PopMap <- gvisGeoMap(WorldPopulation, "Country", "Population.in.millions", "Rank",
#         options=list(dataMode="regions", width=600, height=300))
# plot(PopMap)

Section 2 Load spatial data vectors in dataframes and visualize as map

Install R packages for spatial vectors and list projections

library(sp) # classes and methods for spatial data: points, lines, polygons and grids
suppressPackageStartupMessages(library(rgdal)) # GIS functionality

CRS.WGS84 = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
CRS.UTM   = CRS("+proj=utm +zone=54 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# safecast.file   = "data/Safecast-2011.Rdata"
# load(safecast.file)
# sc.full <- sc.full[sample(1:nrow(sc.full), 2e4, replace=FALSE),]

# sc.time          = as.Date(sc.full$`Captured Time`)
# begin.date       = "2011-05-01"
# decay.date       = "2011-08-01"
# timeFrame = sc.time > as.Date(begin.date) & sc.time < decay.date
# sc   = sc.full[timeFrame,]  # only keep the data for the valid times
# time = sc.time[timeFrame]  # recall to only keep the times the 'valid' dates
# 
# 
# safe.coords     = cbind(sc$Longitude,sc$Latitude)
# cpm             = sc$Value
# safe.spdf       = SpatialPoints(coords = safe.coords, proj4string=CRS.WGS84)
# 
# safe.data = data.frame(time=time, rad = cpm, radLog = log(cpm) )

Create a new spatial data frame with the spatial/temporal subset

# safe.in.spdf = SpatialPointsDataFrame(coords = safe.coords,
#                                       data = safe.data,
#                                     proj4string=CRS.WGS84)
# 
# plot(safe.in.spdf,cex=.1, pch=3)
# 
# writeOGR(safe.in.spdf, ".", paste(basename(safecast.file),"_",decay.date,sep=""), driver="ESRI Shapefile",overwrite_layer=T)


# fn = "rad"
# 
# safe.in.spdf.utm         = spTransform( safe.in.spdf, CRS.UTM )
# safe.in.raster.utm       = raster(safe.in.spdf.utm, resolution=c(50,50))
# safe.in.raster.utm.max   = rasterize(safe.in.spdf.utm, safe.in.raster.utm, field=fn, fun=max)
# 
# safe.in.raster           = raster(safe.in.spdf, resolution=c(50,50))
# safe.in.raster.max       = projectRaster(safe.in.raster.utm.max, safe.in.raster)

# scWeb <- sc[sample(1:nrow(sc), 50, replace=FALSE),]
# 
# scWeb$locationvar = paste(scWeb$Latitude, scWeb$Longitude, sep = ":")
#  
# Webmap <- gvisMap(scWeb, locationvar = "locationvar", "Value", 
#                 options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE,  mapType='satellite', useMapTypeControl=TRUE, width=400,height=800))
# plot(Webmap)

suppressPackageStartupMessages(library(choroplethr))
library(choroplethrAdmin1)
library(ggplot2)
suppressPackageStartupMessages(library(data.table))

data(admin1.map)
japan.map = get_admin1_map("japan")

ggplot(japan.map, aes(long, lat, group=group)) + 
  geom_polygon() + ggtitle("Japan")

data(df_japan_census)
df_japan_census$value=df_japan_census$pop_density_km2_2010
PopJapan = admin1_choropleth("japan", df_japan_census, 
          num_colors = 5, title = '2010 Population Density of Japan', 
          legend = "Population Estimate")

plot(PopJapan)