NCEAS Scientific Computing: Solutions Center Use Case
Convert ESRI Shape Files to KML Files for GoogleEarth / GoogleMaps Display
This use case demonstrates how to convert ESRI Shape Files into Google Earth / Google Maps-compatible KML format, then display the KML files over a Google base map. As a bonus, the case also demonstrates the display of a raster image as Google Earth overlay.
To download a .zip file containing the R script and input files for this example, click here.
Researchers often store spatial data in ESRI Shapefile format, the industry standard for storing and exchanging point, line, and polygon geospatial data. And, increasingly, these researchers want to display these datasets as overlays to backgrounds developed within GoogleMaps and/or GoogleEarth.
Three steps add ESRI shape files to a GoogleMaps / GoogleEarth display:
- 1. Convert each Shape File to Keyhole Markup Language (KML) format.
- 2. Navigate to GoogleMaps and create/select a background map of the area surrounding the KML file features; then, Import the KML file into the map environment. OR,
- 3. Start the GoogleEarth application; File/Open the KML file; GoogleEarth will overlay the features on top of a satellite image of the correct geographic location.
Converting Shape Files to KML format
At this writing (August 2009), R provides three methods for creating KML files:
1. Use kmlLine() or kmlPolygon() functions, included in maptools, to convert line and polygon shape files to SEPARATE kml files.
2. Use the writeOGR() function, included in the rgdal package. This is the most convenient and straightforward method.
3. Use the ogr2ogr() function, from the open-source FWTools package, and the R system() command, to convert shape files to KML format.
Note: GoogleMaps / GoogleEarth display spatial data in unprojected geographic coordinate space. Therefore, in order for your spatial data files to display correctly over the Google background, they must be in the same coordinate systems. In Method 1 / Preliminaries, below, we include the R statements that transform spatial files into unprojected geographic space.
Method 1: Using maptools/kmlLine() and kmlPolygon(): Generate KML files for Line and Polygon shape files
Convert Line and Polygon shape files to KML files using the kmlLine() and kmlPolygon() methods supplied by the maptools package. The advantage of this approach is that it relies solely on 'native' R package routines, available on all R-compatible computer platforms. A disadvantage is that the maptools package does not currently support conversion of Point shapefiles to KML format. Use the second or third methods, below, to convert Point shapefiles to KML format.
Preliminaries: First, we will create the data sets that we use to demonstrate all three examples
#
# rgdal provides readOGR() and writeOGR()
# maptools/sp provides transformation and kml translators
#
library(rgdal)
library(maptools)
#
# Read a polygon shape file, which is in a projected map transformation
#
MyCountiesDF = readOGR("FiveNWStatesCounties.shp","FiveNWStatesCounties")
str(MyCountiesDF)
#
# Re-project the counties into Geographic (unprojected) space.
# Then, write the new polygons to a NEW shape file.
#
proj4string(MyCountiesDF)= CRS("+proj=longlat +datum=WGS84")
MyCountiesGP = spTransform(MyCountiesDF,CRS("+proj=longlat +datum=WGS84"))
writeOGR(MyCountiesGP,".","MyCountiesGP",driver="ESRI Shapefile")
#
# Read the (unprojected) Shape File that we created in two ways:
# 1) As a SpatialPolygonsDataFrame - to generate KML file
# 2) As a Map object - To generate polygon centroids for use in our next example.
#
MyCountiesDF = readOGR("MyCountiesGP.shp","MyCountiesGP")
#
# Read the shape file ino a maptools/Map object and generate centroids for each polygon.
#
# Finally, create a projected SpatialPointsDataFrame object from the centroid coordinates
# We sill use this in the second and third examples.
#
MapPolysGP = read.shape("MyCountiesGP")
MapPolyCents = get.Pcent(MapPolysGP)
#
ProjString = CRS(proj4string(MyCountiesGP))
CentCoords = matrix(c(MapPolyCents[,1],MapPolyCents[,2]),nrow=length(MapPolyCents[,1]),ncol=2,byrow=FALSE)
ProjectedCentroids = spTransform(SpatialPoints(coordinates(CentCoords),proj4string=ProjString),ProjString)
#
# Create the @data slot for the centroid Spatial Points Data Frame using the County polygon file's attributes
#
CentroidAttrs = data.frame(MapPolysGP$att.data$COUNTYNAME,MapPolysGP$att.data$COUNTYFIPS)
CentroidsSPDF = SpatialPointsDataFrame(CentCoords,CentroidAttrs,coords.nrs=numeric(0),ProjString)
writeOGR(CentroidsSPDF,".","MyCountyCentroidsGeogProj",driver="ESRI Shapefile")
#
# End of preliminaries!
#Now, convert the polygon shape file into a KML file
# 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.
#
CountyPolys = slot(MapPolysGP[MapPolysGP$STATEFIPS %in% c("53"),], "polygons")[[1]]
kmlPolygon(CountyPolys, kmlfile="NWStateCountyPolys.kml", name="FiveNWStateCounties", col="blue", lwd=1,
border=1, kmlname="R Test",
kmldescription="KML File:NW USA State Counties")
#
# Note: August 2009 - A defect in kml_xxx tools writes only the first feature written to the kml file.
# workaround (until defect is fixed): Use writeOGR() instead
#
writeOGR(MapPolysGP, "MyCountiesGP.kml", "MyCountiesGP", driver="KML")Use the same R statements to produce a KML file for a Line shape file which is ALREADY in unprojected geographic coordinate space.
NevadaStreamsGP = readOGR("NevadaStreamsGP.shp","NevadaStreamsGP")
#
# Note: August 2009 - Because of a defect in kml_xxx tools, only the first feature is written to kml file.
#
kmlLine(NevadaStreamsGP, kmlfile="NevadaStreams.kml", name="NevadaStreams", col="green", lwd=1,
kmlname="R Test",kmldescription="KML File: Nevada Streams")
#
# workaround (until defect is fixed): Use writeOGR() instead
#
writeOGR(NevadaStreamsGP, "NevadaStreams.kml", "NevadaStreams", driver="KML")
#Click here to see the vector (counties and centroids) shape files displayed in GoogleEarth
Using maptools/kmlOverlay(): Overlay Raster Image over a GoogleEarth background
For instances when analysts wish to display a raster image (e.g., satellite image or simulation model result) as an overlay to GoogleMaps or GoogleEarth display, maptools provides the kmlOverlay() metnod. This routine accepts a byte/integer raster image in PNG (Portable Network Graphics) format, and generates the KML 'wrapper' file that facilitates display of the image within GoogleEarth or GoogleMaps. The following example converts a floating point digital elevation model image to byte/integer format, then creates the KML wrapper.
#
# rgdal provides readOGR() and writeOGR()
# maptools/sp provides transformation and kml translators
#
# 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("NevadaSiteDEMRsmp.img")
rng =range(DEMImg$band1)
#
# Translate floating point to integer-byte
#
DEMImg$band1 = as.integer((DEMImg$band1 / rng[2]) * 256)
writeGDAL(DEMImg,"DemByteImg.Tif",driver="GTiff") # normalized DEM data
#
# Create a PNG version of our byte /Tif image
# writeGDAL can not create a PNG, so we use
# an 'indirect' method that works...
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,"NevadaDemImg.kml","NevadaDemImg.png",name="First R Dem Kml")
}Click here to see the raster (terrain) image and vector (streams) shape files displayed in GoogleEarth
Method 2: Use writeOGR() function from rgdal package:
A versatile method for converting ESRI Shape Files to KML files uses the writeOGR() function, supplied by the R rgdal package. writeOGR() accepts as input 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 GoogleMaps and/or GoogleEarth.
Note: writeOGR() works only with Shape Files that 1) are in unprojected, geometric coordinate spate, and 2) whose spatial coordinates section contain only x (longitude) and y (latitude) components. Shape Files containing the z (altitude) component will be rejected.
#
# Demonstrate conversion of point and polygon shape files into KML format using writeOGR()
# load the appropriate R package libraries
#
library(rgdal)
#
# Read polygon and point shape files using the same method
#
MyCountiesDF = readOGR("MyCountiesGP.shp","MyCountiesGeogP")
MyCentroidsDF = readOGR("MyCountyCentroidsGP.shp","MyCountyCentroidsGP")
#
# Create KML files for the point, line, and polygon data objects. Just by specifying the KML file driver.
#
writeOGR(MyCountiesDF, "MyCountiesDF.kml", "MyCountiesDF", driver="KML")
writeOGR(NevadaStreamsGP, "NevadaStreams.kml", "NevadaStreams", driver="KML")
writeOGR(MyCentroidsDF, "MyCentroidsDF.kml", "MyCentroidsDF", driver="KML")
#
# End of example.
#Method 3: Use ogr2ogr utility from FWTools/ OGR package: Access a 'standalone' utility using R system() function
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 ogr2ogr, demonstrated here.
While the previous two examples on this page are 'more than adequate' methods for KML file creation, this example demonstrates the R interface to the substantial body of geospatial software tools that provide capabilities not (yet) included in the R Geospatial tools.
# Advantage of accessing fwtools library from R: Many powerful file conversion utilties available, on several platforms.
# Also, you can work outside of R's restrictive memory space.
# Disadvantage: In some cases, execution speed is sacrificed.
#
# Using the ogr2ogr function from the fwtools package to translate
# our shape files to KML format.
# For this example to work, the host computer must have:
#
# 1) The fwtools file conversion utilities 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.2\bin')
#
sOgr2OgrArgs = sprintf("-dsco NameField=%s","OurTestKMLRegion")
sCommandString = sprintf(" -f KML %s %s %s","MyCountiesGPOgr2ogr.kml","MyCountiesGP.shp",sOgr2OgrArgs)
system(paste('"ogr2ogr"',sCommandString))
Learning More:
The CRAN (R) Task View summarizes the R Analytical Environment's spatial data analysis and management tools.
This edition of the R Newsletter contains an in-depth article describing the R spatial data classes; parts of the article appear here.
This ESRI White Paper describes the Shape File format.
Search for information (including PROJ4 string) about your desired Map Projection: EPSG spatialreference.org Web Site
Alternate PROJ4 information: remotesensing.org Geotiff projection Web Site
This Web Site describes the standard .dbf format which is a component of the ESRI Shape File.