Display R Spatial Objects With Google Maps/Google Earth

This Use Case creates Google Maps / Google Earth-compatible Keyhole Markup Language (KML) files from R Spatial objects containing data from vector ESRI Shape Files and raster GeoTiff files. The case demonstrates three techniques for writing R Spatial objects into KML files, and one technique for reading KML files into R Spatial objects. Also revealed: reading polygon KML files into R Spatial objects.

This case also introduces the functionality of the recently released RGoogleMaps package. We demonstrate use of this package to download a static GoogleMaps map image and use it as the background for R spatial object plots. Click here to step directly to this section of the Use Case.

Download data and R Code for this Use Case

Project Requirements:

Within R, read vector ESRI Shape Files and a raster GeoTIFF file into R Spatial objects, then create Keyhole Markup Language (KML) files for display in GoogleMaps / Google Earth.

To add R Spatial objects to GoogleMaps / GoogleEarth sessions:
  • Read each ESRI Shape File into R Spatial objects, then translate and write them into KML files.
  • Within GoogleMaps create/select a background map containing the objects in each file. Then, import the KML file into the map. OR,
  • Within GoogleEarth, simply open the KML file; GoogleEarth will read the KML file's spatial metadata and overlay the features on top of a satellite image of the underlying location.
Creating KML files for R Spatial objects: Three techniques
  1. kmlLine()or kmlPolygon()functions, included in maptools package, to convert line and polygon objects to SEPARATE kml files.
  2. writeOGR() function, included in the rgdal package. This is the most convenient and straightforward method.
  3. The R system() command invokes the ogr2ogr() method,  provided by the open-source FWTools  package, to create KML files from Shape Files (Note: requires installation of FWTools library on host computer).
Technique 1: maptools package methods: kmlLine(), kmlPolygon(), and kmlOverlay()

This method uses R maptools package functions; currently, maptools has no methods that convert Point objects to KML format. Instead, use writeOGR, used here and in Technique 2, to write Point objects to a KML file..

Note: To most accurately align R Spatial object features with corresponding Google Maps/Google Earth features, the R objects should be transformed into the map projection used by Google. Google Earth uses a Cylindrical Latitude/Longitude projection; Google Maps uses a variation of the Spherical Mercator projection. These projections are represented by EPSG and Proj.4 codes used by R spatial coordinate transformation methods. The Technique 1 script for this case gives the codes for both projections used by Google, and transforms the Spatial objects into the Google Earth-compatible coordinate system.

First, create the R Spatial objects used in all three examples:

# Generate KML files using three methods: 
# Method 1: KmlPolygon() / KmlLine() functions in Maptools package
# Method 2: ReadOGR() / WriteOGR() functions from rgdal package
# Method 3: ogr2ogr utility from FWTools package, executed using the R system() function # Bonus Method: Reading a KML file into an R Spatial object
# rgdal provides readOGR() and writeOGR()
# maptools/sp provides transformation and kml translators
# PBSmapping provides polygon centroid creation function

library(rgdal) # loads 'dependant' sp package
library(maptools)

# Remember to set the current working directory to the folder into which you unpacked this archive.

# setwd("C:/Projects/GeospatialUseCasesMar2009/CreateKMLFiles/dl") #example

# Re-project the counties into Geographic (unprojected) space.
# Then, write the new polygons to a NEW shape file.

# Read a polygon shape file, which is in a Albers Equal Area map transformation

NwCountiesAE <- readOGR("FiveNWStatesCounties.shp","FiveNWStatesCounties")

# Reference: For GoogleMaps overlays, project data into EPSG:3857 (Spherical Mercator projection)
# For this Use Case, we will project all Shape Files for GoogleEarth overlay
epsg3857String <- CRS("
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
# For GoogleEarth overlays, reproject data into EPSG:4326 (simple Latitude/Longitude 'geographic' projection)

epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
NwCountiesGMap <- spTransform(NwCountiesAE,epsg4326String)
writeOGR(NwCountiesGMap,".","NwCountiesGMap",driver="ESRI Shapefile")

# Re-read the (reprojected) Shape File with readOGR to verify the projection information

NwCountiesGMap <- readOGR("NwCountiesGMap.shp","NwCountiesGMap")

# Read a point shape file containing centroids for all Northwest counties

MapPolysCent <- readOGR("NwCountiesCentroids.shp","NwCountiesCentroids")

# We pass ONLY the Polygons component of the SPDF to kmlPolygon, so extract it.
# Create the KML with only Idaho counties (Idaho County FIPS code: 53)
# Extract the Polygons component of the SpatialPointsDataFrame,
# as needed by the kmlToPolygon() function.

# Note: May 2010 - kmlPolygon only writes the first spatial feature
# (with Polygon ID 0) from a Spatial...DataFrame. To capture all of
# the spatial features in this file, need to 'drill down' to the 'polygons' layer)
# CountyPolys <- slot(NwCountiesGP[NwCountiesGP$STATEFIPS %in% c("53"),], "polygons")[[1]]
# kmlPolygon(CountyPolys, kmlfile="NWStateCountyPolys.kml", name="NWStateCountyPolys",
# col="blue", lwd=1,
# border=1, kmlname="R Test",kmldescription="KML File:NW USA State Counties")

sw <- NwCountiesGMap[NwCountiesGMap$STATEFIPS %in% c("53"),] # extract desired SPDF subset
out <- sapply(slot(sw, "polygons"), function(x) {
kmlPolygon(x,
name=as(sw, "data.frame")[slot(x, "ID"), "STATEFIPS"],
col="red", lwd=1.5, border='black',
description=paste("ISO3:", slot(x, "ID")))
}
) # outside of 'sapply' function call
tf <- "MyNewKMLFILE.kml"
kmlFile <- file(tf, "w")

cat(kmlPolygon(kmlname="R Test", kmldescription="Hello")$header,file=kmlFile, sep="\n")
cat(unlist(out["style",]), file=kmlFile, sep="\n")
cat(unlist(out["content",]), file=kmlFile, sep="\n")
cat(kmlPolygon()$footer, file=kmlFile, sep="\n")
close(kmlFile)

# Alternate method: use the GDAL tools writeOGR() to convert and write the shapefile in KML format.

writeOGR(NwCountiesGMap, "NWStateCountyPolysOGR.kml", "NWStateCountyPolysOGR", driver="KML")

Next,  generate KML files for POINT and LINE shape files.

# Creating a KML file for the POINT shape file we read: The only way within R is to use GDAL / OGR utilities  
MapPolysCentGMap <- spTransform(MapPolysCent,epsg4326String) writeOGR(MapPolysCentGMap,".","NwCountyCentGMap",driver="ESRI Shapefile") writeOGR(MapPolysCentGMap, "NwCountyCentroids.kml", "NwCountyCentroids.kml", driver="KML")

# Next, create a KML file for a set of LINES (streams)

NevadaStreamsGP <- readOGR("NevadaStreamsGP.shp","NevadaStreamsGP")
# Same issue: only first feature (with Polygon ID 0) is written, so use GDAL OGR methods

NevadaStreamsGMap <- spTransform(NevadaStreamsGP,epsg4326String)
out <- sapply(slot(NevadaStreamsGMap, "lines"), function(x) {
kmlLine(x,
name=as(sw, "data.frame")[slot(x, "ID"), "GNIS_NAMEhiso"],
col="green", lwd=1,
description=paste("ISO3:", slot(x, "ID")))
}
) # outside of 'sapply' function call
tf <- "MyNewKMLLineFile.kml" ##tempfile()
kmlFile <- file(tf, "w")

cat(kmlLine(kmlname="R Test", kmldescription="Hello")$header,file=kmlFile, sep="\n")
cat(unlist(out["style",]), file=kmlFile, sep="\n")
cat(unlist(out["content",]), file=kmlFile, sep="\n")
cat(kmlLine()$footer, file=kmlFile, sep="\n")
close(kmlFile)

# OR: call to GDAL/writeOGR

writeOGR(NevadaStreamsGMap, "NevadaDesertStreams.kml", "NevadaDesertStreams", driver="KML")

Click here to see the vector (counties and centroids) shape files displayed in GoogleEarth

Overlay a raster image on a Google Earth background: kmlOverlay()

To overlay a raster image (e.g., satellite image or simulation model result) on a Google Maps or Google Earth display, use maptools: kmlOverlay() to prepare the file.

kmlOverlay() reads a byte/integer PNG (Portable Network Graphics) format raster image, and creates the KML 'wrapper' file containing spatial metadata used by Google Earth.

To display the KML-wrapped image in Google Earth: Install Google Earth exectutable on your computer, then double-click the image's KML 'wrapper' file.

To display the KML-wrapped image in Google Maps: The image and wrapper files must be Web-accessible:

  1. Start Google Maps: Go to http://maps.google.com
  2. Within Google Maps, paste KML wrapper file's Web address into the 'Search Maps' box.

For example, try this link:

http://www.nceas.ucsb.edu/files/scicomp/Dloads/ExamplePrograms/NevadaDEMImageOverlay.kml

# Use kmlOverlay() to create a KML 'wrapper' for a Digital Elevation Map (DEM) raster image file
#
# Convert a floating-point DEM elevation map image
# into a KML file complex for display with GoogleEarth.
# Begin with Digital Elevation Map Iamage, floating point elevation values, .IMG format.
# Note: The input DEM image is a .tif file, with elevations in floating point format.
# Steps: 1) Convert image to a TIFF format integer raster with pixre values > 0 < 256
# 2) Convert the TIFF file to PNG image format.
# 3) Create a GoogleEarth-compatible Spatial Grid object from the PNG image.
# 4) Use maptools/KmlOverlay() to create the KML file 'complex' for GoogleEarth display.

DEMImg <- readGDAL("NevadaSiteDemGP.img")
rng =range(DEMImg$band1)

# Translate floating point to integer-byte

DEMImg$band1 <- as.integer((DEMImg$band1 / rng[2]) * 256)
writeGDAL(DEMImg,"DemByteImg.tif",type="Byte",driver="GTiff") # normalized DEM data

# Create a PNG version of our byte /Tif image
# writeGDAL can not create a PNG, so we use
# GDAL library routines to ceate, and then save
# as a .png image, a temporary dataset.

DemTif <- GDAL.open("DemByteImg.tif")
xx <- copyDataset(DemTif,driver="PNG")
saveDataset(xx,"NevadaDEMImg.png")
DemPng <- readGDAL("NevadaDEMImg.png") # na proper 'integer' PNG file

# Preprocessing step: create a special SpatialGrid object
# for display in GoogleEarth

DemPngGK <- GE_SpatialGrid(DemPng)

# Generate the KML 'wrapper' for the png file.
# 'GE-compatible' PNG file is now a complex of three files:
# .png, .kml, and .png.aux.xml files.

kmlOverlay(DemPngGK,"NevadaDEMImageOverlay.kml","NevadaDEMImg.png",name="First R Dem Kml")

Click here to see the raster (terrain) image and vector (streams) shape files displayed in GoogleEarth.

Technique 2: Use rgdal package writeOGR() function:

A versatile method for converting ESRI Shape Files to KML files uses the R rgdal package writeOGR() method.

writeOGR() inputs are: the root name of an ESRI Shapefile, the name of the KML file to be generated, and the name of the file conversion driver (for conversion to KML, this will always be set to "KML") to be used.

writeOGR() automatically detects the shape file type (point/line/polygon) and generates the appropriate KML file, ready for use with Google Maps and/or Google Earth.

Note: writeOGR() works only with shape files that 1) are in geographic coordinate space with x (longitude) and y (latitude) values.

*** Note: Shape Files containing the z (altitude) component will not be processed. ***

	
# Demonstrate conversion of point and polygon shape files into KML format using writeOGR()
# load the appropriate R package libraries
# Read polygon and point shape files using the same method

NwCountiesGMap = readOGR("
NwCountiesGMap.shp","NwCountiesGMap")
NwCentroidsGMap = readOGR("
NwCountyCentGMap.shp","NwCountyCentGMap")

# Create KML files for the point, line, and polygon data objects. Just by specifying the KML file driver.

writeOGR(NwCountiesGMap, "NwCountiesGMap.kml", "NwCountiesGMap", driver="KML")
writeOGR(NevadaStreamsGMap, "NevadaStreamsGMap.kml", "NevadaStreamsGMap", driver="KML")
writeOGR(MyCentroidsGMap, "MyCentroidsGMap.kml", "MyCentroidsGMap", driver="KML")

# End of example.

Technique 3: Use R system() function to call FWTools toolkit ogr2ogr utility

FWTools is an Open Source GIS utility package providing many useful geospatial utilities.

Within FWTools, the Geospatial Data Abstraction Library (GDAL) includes file conversion and manipulation tools, including the ogr2ogr utility, demonstrated here.

This example demonstrates the R interface to the substantial body of geospatial software tools not (yet) included in the R Geospatial tools.

# Use the ogr2ogr function from the fwtools package to translate 
# the 'Northwest Washington Counties' polygon shape file into a KML file.
# For this example to work, the host computer must have:
#
# 1) The fwtools file conversion utilitues installed.
# Download the package at: http://fwtools.maptools.org/

# 2) Add the path to the fwtools ./bin folder to your system 'path' environment variable.
# (on my system, this path is: 'C:/Program Files (x86)/FWTools2.4.7/bin')

sOgr2OgrArgs <-sprintf("-dsco NameField=%s","OurTestKMLRegion")
sCommandString <-sprintf(" -f KML %s %s %s","NwCountiesGeogProjOgr2ogr.kml","NwCountiesGMap.shp",sOgr2OgrArgs)

# Execute the operating system command
# on computers that do not (yet) host the FWTools package, uncomment the line AFTER installing FWTools.
# Note: This line is commented out to avoid error message ('FWtools not found')
# system(paste('"C:/Program Files (x86)/FWTools2.4.7/bin/ogr2ogr"',sCommandString))

Reading KML files into R Spatial objects:
  1. To use this technique, install the Expat XML parsing library on your computer. Be sure that the library is visible to the R and GDAL executable files via the PATH environment variable, and that your R script includes the 'require(rgdal)' statement. 
  2. Use the ogrInfo() command to verify existence of KML driver; use the readOGR() command to read the KML file (Using the KML driver and Expat XML library) into the appropriate R spatial object. The readOGR() command 'autodetects' the file type using the '.kml. suffix, and loads the appropriate software driver to input the file into the appropriate R Spatial object.
# READING a KML file into R: check for presence of OGR KML driver, then read the file

ogrInfo("NWStateCountyPolysOGR.kml",layer="NWStateCountyPolysOGR")
polysFromKML=readOGR("NWStateCountyPolysOGR.kml",layer="NWStateCountyPolysOGR") class(polysFromKML)
plot(polysFromKML)
Using RGoogleMaps package to get map images for plot backgrounds:

The RGoogleMaps package is an interface between any R script and the Google Map server. It equips the R program to download static (i.e., non-scrolling) map images from Google Maps, using spatial coordinates and other parameters, and use the images as R plot backgrounds. The following script extracts two Google Maps for our study area, then uses them as background images for polygon/point plots:

Click here to see the first map; click here to see the second map.



# Demonstrates use of RgoogleMaps package to grab GoogleMaps map/image
# for user-specified geographic rectangle.

   require(RgoogleMaps)
   require(rgdal)      # required if the 'png' image file format is to be used....
   require(PBSmapping) # for PolySet type, to  pass shape files to  PlotPolysOnStaticMap()
   require(ReadImages) # includes read.jpeg(), used by GetMap to write JPEG-format output image

# Using PBSmapping routine, read polygon Shape File to overlay on the Google image

   polyShpFile="./FiveNWStatesCounties.shp"
   shpPolySet=importShapefile(polyShpFile,projection="LL");
   pointShpFile="./MyCountyCentroidsGP.shp"
   shpPointSet=importShapefile(pointShpFile,projection="LL");

# Specify study area (Northwest USA) lower left-hand and upper right-hand corner coords.

   lats = c(40.9,49.0)
   lons = c(-125.0, -104.0)
   centerPt = c(mean(lats),mean(lons))
   boundBox = qbbox(lat = lats,lon = lons)

   mapFromBbox = GetMap.bbox(boundBox$lonR,boundBox$latR,maptype="terrain",destfile="MapFromBbox.png")
   
# Display the GoogleMaps 'base map', then overlay polygons (outlines ONLY, 'col=0') on map
   
   PlotPolysOnStaticMap(mapFromBbox, shpPolySet, lwd=1.5, col=0, add = F)

# Compute the maximum (most detailed) zoom level that will contain study area.
# Use this with the computed 'center point' of the user points to extract the
# desired map using GetMap() function.

   zoomFact = min(MaxZoom(range(lats),range(lons)))
   
   mapFromCenPt = GetMap(center=centerPt,zoom=zoomFact,maptype="roadmap",destfile="MapFromCentPt.jpg")

# Overlay the point set on top of the second, "roadmap" base map.
# Then, add the polygons (only) displayed on the original map.

   PlotOnStaticMap(mapFromCenterPt, shpPointSet$Y,shpPointSet$X, col="red",pch='*', add = F)
   PlotPolysOnStaticMap(mapFromCenterPt, shpPolySet, lwd=1.5, col=0, add = T)

# August 2010: For this version of RgoogleMaps, here is the method for writing the new map to a file
# (until the destfile=fn argument is implemented).
   
   png("./MapFromBbox.png", 640, 640);
   PlotOnStaticMap(mapFromCenterPt, shpPointSet$Y,shpPointSet$X, col="red",pch='*', add = F)
   PlotPolysOnStaticMap(mapFromCenterPt,shpPolySet, lwd=1.5, col=rgb(0,0,0,0), add = T)
   dev.off() # end of example

Learning More:

The CRAN (R) Task View summarizes the R Analytical Environment's spatial data analysis and management tools.

Search for information (including PROJ.4 string) about your desired Map Projection: EPSG spatialreference.org Web Site

This Use Case was updated August, 2010.

 

scicomp: