Installation of plotKML and first steps

This tutorial reviews plotKML package functionality with examples of inputs and outputs. The plotKML (R package) project summary page you can find here. To learn more about the Global Soil Information Facilities (GSIF), visit the main project page. This tutorial only explains how to quickly plot various spatial and spatio-temporal data in Google Earth. To use the advanced functionality, refer to the plotKML package documentation.

Download the tutorial as R script. Watch this tutorial on Youtube.

Installation and first steps

KML and Google Earth

Google Earth is probably the biggest publicly accessible GIS in the world, and certainly the most popular browser for geographic data. All displays in Google Earth are controlled by KML files, which are written in the Keyhole Markup Language developed by Keyhole, Inc. KML is an XML-based language for managing the display of three-dimensional geospatial data. This format is now used in several geographical browsers (Google Maps, ArcGIS Explorer, World Wind, KDE Marble and similar) and is one of the most used OGC standards for spatial data.

The KML files can be edited using an ASCII editor (as with HTML), but they can also be edited directly in e.g. Google Earth. KML files are very often distributed as KMZ files, which are just zipped KML files with a *.kmz extension. This interactive KML sampler will give you some good idea about what is all possible via KML.

Although there are several browsers for KML files, the examples shown in this tutorial are based on Google Earth. To install Google Earth, run the Google Desktop installer that you can obtain from the Google's website. To start a KML file, just double-click it and the map will be displayed using the default settings. Other standard Google's background layers, such as roads, borders, places and similar geographic features, can be turned on of off using the Layers panel. There is also a commercial Plus and Pro versions of Google Earth, but for purpose of exercises in this guide, the free version is more than enough.

An important note: loading your data from R to Google Earth can be time consuming. First, you need to understand that there is a difference between loading the vector and raster maps to Google Earth. Typically, it is relatively easy to load vector data such as points or lines to Google Earth, and somewhat more complicated to do the same with raster maps. Also note that, because Google Earth works exclusively only with Latitude/Longitude projection system WGS84 ellipsoid), all vector/raster maps that are not in this system need to be first reprojected before they can be exported. Be aware that there are many other packages that allow importing GIS layers such as shape files, GPX tracks or GeoTiFFs to Google Earth. This tutorial focuses only on writing KML files from R via the plotKML package, but there are also several alternative open source solutions for parsing KML files (e.g. pyKML).

Also note that there are some serious limitations in what you can and cannot do with scientific visualizations produced using Google Earth (i.e. screenshots). Even though most of the screen can be covered with your own data and R plots in Google Earth window, Google still might require that you sign some copyright agreement or request a permission to use their software for commercial applications. For more info please refer to the Google geoguidelines.

Package installation and main functionality

plotKML can be obtained from CRAN. The package is still being developed and filtered for bugs (until version >1.0), so to obtain the most recent development version of the package run:

> sessionInfo()
R version 2.14.1 (2011-12-22)
Platform: x86_64-pc-mingw32/x64 (64-bit)
> install.packages("plotKML", repos=c("http://R-Forge.R-project.org"))
trying URL 'http://R-Forge.R-project.org/bin/windows/contrib/2.15/plotKML_0.3-3.zip'
Content type 'application/zip' length 1996905 bytes (1.9 Mb)
opened URL
downloaded 1.9 Mb

Note that, when installing package from source i.e. from R-forge using R CMD, dependencies need to be taken care of before starting the package.

If the installation was successfully, you can start the package by typing:

> library(plotKML)
plotKML version 0.3-3 (2013-02-04)
URL: http://plotkml.r-forge.r-project.org/

Note that plotKML has many dependencies. It largely builds on the sp and raster-related packages, but it tries to serve a number of packages that either generate space-time data or use them as an input to spatial analysis. The drawback of using many packages is that the package takes more time to load, and it is more sensitive to changes in the packages that it extends and depends on.

plotKML is largely based on the XML package. It basically consists of a number of low-level and wrapper functions to write space-time classes from R to KML. To learn more about KML consider reading Wernecke's (2009) KML Hanbook.

External applications and environmental settings

plotKML uses also a number of optional external packages that need to be installed independently: FWTools / GDAL utilities, SAGA GIS, ImageMagick, Python, and (of course) Google Earth. Before visualizing any data in plotKML, we advise you to first install any additional R and external software package that can speed up the processing and help you achieve the maximum result. These packages are NOT an requirement — plotKML is suppose to be able to run most of operations without them, nevertheless, we recommended that you considet installing and registering on your machine all four external software.

You can check if R has successfully located the external packages by using the plotKML.env() command. For example, if you are running Windows OS, then the output of finding the paths would look something like this:

> paths()
Loading required package: animation
Version: ImageMagick 6.6.8-10 2011-03-26 Q16 http://www.imagemagick.org
Located FWTools from the Registry Hive: "C:\PROGRA~2\FWTOOL~1.7"
Located Python from the Registry Hive: "C:\Python27\"
                                       gdalwarp
1 "C:\\PROGRA~2\\FWTOOL~1.7\\bin\\gdalwarp.exe"
                                 gdal_translate
1 "C:\\PROGRA~2\\FWTOOL~1.7\\bin\\GDAL_T~1.EXE"
                                                  convert
1 "C:\\Program Files\\ImageMagick-6.6.8-Q16\\convert.exe"
                      python saga_cmd
1 "C:\\Python27\\python.exe"
Warning message:
In paths() :
  Could not locate SAGA GIS! Install program and add it to the Windows registry. See http://www.saga-gis.org/en/ for more info.

This means that plotKML has managed to locate all external applications (except for SAGA GIS), which means that the package can be used up to 100% of its functionality. Note that, plotKML has to re-locate external applications every time new R session is started. To permanently preserve the plotKML environmental settings, you either need to save your session, or change settings in the Profile.site.

To access some standard parameter from plotKML, you can point to the plotKML options. For example, to obtain the default coordinate system (WGS84) used by Google Earth use:

> get("ref_CRS", envir = plotKML.opts)
[1] "+proj=longlat +datum=WGS84"

To further customize the plotKML options (on-load), consider putting: library(plotKML); plotKML.env(…, show.env = FALSE) in your /etc/Rprofile.site.

In principle, all KML parsing functions in plotKML are implemented via the corresponding S4 classes. Consequently, all kml_layer methods require predefined structure in the data i.e. the data first needs to be converted to the corresponding classes.

Writing sp, raster and derived classes to KML

Point, line, polygon and gridded objects

plotKML largely mimics the existing plotting functionality available for spatial data in R; e.g. the spplot functionality available via the sp package (Bivand et al., 2008) and raster plotting functionality available via the raster package (Hijmans, 2012). For example, in the sp package, to plot points we would run:

> data(eberg)
> coordinates(eberg) <- ~X+Y
> proj4string(eberg) <- CRS("+init=epsg:31467")
> eberg <- eberg[runif(nrow(eberg))<.2,]
> bubble(eberg["CLYMHT_A"])

or:

> spplot(eberg.ll["CLYMHT_A"], edge.col="black", alpha=0.8, cex=seq(.3,3,length=5))

which in plotKML looks very similar:

> plotKML(eberg["CLYMHT_A"])
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg__CLYMHT_A__.kml


operating instructions

Various functions are available as part of the map.

mouse/touch operation
  • moving by grabbing the map with a mouse-click you can move the map around
  • overview map using the + button in the bottom right of the map you can expand an overview map
  • zooming in and out using the + and - buttons in the top left of the map you can obtain more or less detail in the map
  • switching themes or maps clicking on the ≡ icon on the right-hand side of the map you can view and select available maps and themes
  • retrieving information the map may contain elements that contain more information, by clicking these a popup will show this information
  • fullscreen display using the ✈ button the map can be maximized to fullscreen display, use the ✕ button to return to page display.
keyboard operation

Keyboard operation becomes available after activating the map using the tab key (the map will show a focus indicator ring).

  • moving using the arrow keys you can move the map
  • overview map using the + button in the bottom right of the map you can expand an overview map
  • zooming in and out using the + and - buttons in the top left of the map or by using the + and - keys you can obtain more or less detail in the map
  • switching themes or maps clicking on the ≡ icon on the right-hand side of the map you can view and select available maps and themes
  • retrieving information the map may contain elements that contain more information, using the i key you can activate a cursor that may be moved using the arrow keys, pressing the enter will execute an information retrieval. press the i or the escape key to return to navigation mode
  • fullscreen display using the ✈ button the map can be maximized to fullscreen display, use the ✕ button to return to page display.

It's possible that some of the functions or buttons describe above have been disabled by the page author or the administrator

 

Points of Interest
id symbol latitude longitude description
KML file KML track KML track: eberg.ll


the only difference is that the plotKML command opens the default KML-viewing software (in most of the times Google Earth) and not just an R plot window. Note that, unlike with sp and other spatial packages, before writing object to KML, plotKML tries to reproject any spatial object to the referent coordinate system. For practical purposes, the file and folder names in the KML file will be generated from the object name. The visual elements can be further customized by triggering some of the aesthetics parameters such as:

  • colour - specifies the color that is used to fill each spatial element (black if not specified otherwise);
  • shape - specifies icons to be used to display each spatial element;
  • width - determines the width of lines (1 by default),
  • labels - specifies text to be added as labels to each spatial element;
  • altitude - adds altitude to each spatial element (0 if not specified otherwise),
  • size - specifies size of icons in the display;
  • balloon - specifies whether to attach attribute data to each placemark;
Bubble-type plots in R and a similar plot produced using the plotKML (shown in Google Earth).

This results in a slightly richer aesthetics that the bubble plot in R. To remove the labels and fill in circles using e.g. yellow color only, we can set:

> plotKML(eberg["CLYMHT_A"], colour_scale=rep("#FFFF00", 2), points_names="")
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg__CLYMHT_A__.kml
Bubble-type plots in R after turning off labels and using a single colour display.

Note that, by default, colors need to be submitted in the hex format. A number of standard built-in GIS palettes are also available.

Lines are by default plotted without any attributes:

> data(eberg_contours)
> plotKML(eberg_contours)
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_contours.kml

To plot the attribute values use e.g.:

> plotKML(eberg_contours, colour=Z, altitude=Z)
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_contours.kml
Contour lines with actual altitudes attached.

When plotting SpatialPolygons (as in the case of points), both attributes and labels will be plotted by default:

> data(eberg_zones)
> plotKML(eberg_zones["ZONES"])
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_zones__ZONES__.kml
Polygon map with labels.

To a polygon map plot, you can also attach some attribute information vai the altitude parameter, which will produce visualizations similar to the ones produced by the GraphEarth software:

> length(eberg_zones)
[1] 11
> plotKML(eberg_zones["ZONES"], altitude=runif(length(eberg_zones))*500)
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_zones__ZONES__.kml
Polygon map with attributes attached as altitude.

To plot RasterLayer, SpatialPixelsDataFrame, and/or SpatialGridDataFrame-class objects, you again need to specify the target layer:

> data(eberg_grid)
> gridded(eberg_grid) <- ~x+y
> proj4string(eberg_grid) <- CRS("+init=epsg:31467")
> data(SAGA_pal)
> plotKML(eberg_grid["TWISRT6"], colour_scale = SAGA_pal[[1]])
KML file header opened for parsing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_grid__TWISRT6__.kml
Compressing to KMZ...

This figure basically mimics SAGA GIS (Böhner et al., 2006) because it shows a DEM parameter derived in SAGA GIS (SAGA Topographic Wetness Index), and it uses the default SAGA GIS legend.

SAGA GIS-derived Topographic Wetness Index visualized in Google Earth.

To plot a factor-type map, you will also need attach an appropriate color legend. To plot for example land cover classes we can use:

> eberg_grid$LNCCOR6 <- as.factor(paste(eberg_grid$LNCCOR6))
> levels(eberg_grid$LNCCOR6)
[1] "112 - Discontinuous urban fabric"
[2] "211 - Non-irrigated arable land"
[3] "231 - Pastures"
[4] "242 - Complex cultivation patterns"
[5] "311 - Broad-leaved forest"
[6] "312 - Coniferous forest"
[7] "313 - Mixed forest"

As you can notice, these are the classes from the Corine Land Cover 2006 map for Europe. We can use the built-in legend for the land cover classes:

> data(worldgrids_pal)
> pal = as.character(worldgrids_pal["corine2k"][[1]][c(1,11,13,14,16,17,18)])
> pal
"#E5481B" "#EEEB9A" "#E0DC45" "#F2D442" "#6BBB4A" "#00A046" "#46AF48"
> plotKML(eberg_grid["LNCCOR6"], colour_scale=pal)
KML file header opened for parsing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg_grid__LNCCOR6__.kml
CORINE Land Cover classes plotted in Google Earth with legend added as a screen overlay.

Field photographs

plotKML can also be used to visualize field photographs of e.g. soil profiles, trees, land cover and similar. This possible via the class SpatialPhotoOverlay. If you have a photograph that can be georeferenced, you can upload it to WikiMedia portal, and then create an object of class SpatialPhotoOverlay and visualize it from R by running:

> imagename = "Soil_monolith.jpg"
> x1 <- getWikiMedia.ImageInfo(imagename)
> sm <- spPhoto(filename = x1$url$url, exif.info = x1$metadata)
> str(sm)
Formal class 'SpatialPhotoOverlay' [package "plotKML"] with 5 slots
  ..@ filename    : chr "http://upload.wikimedia.org/wikipedia/commons/3/3d/Soil_monolith.jpg"
  ..@ pixmap      :Formal class 'pixmapRGB' [package "pixmap"] with 8 slots
  .. .. ..@ red     : num[0 , 0 ]
  .. .. ..@ green   : num[0 , 0 ]
  .. .. ..@ blue    : num[0 , 0 ]
  .. .. ..@ channels: chr [1:3] "red" "green" "blue"
  .. .. ..@ size    : int [1:2] 0 0
  .. .. ..@ cellres : num [1:2] NaN\ NaN\
  .. .. ..@ bbox    : num [1:4] 0 0 0 0
  .. .. ..@ bbcent  : logi FALSE\
  ..@ exif.info   :List of 21
  .. ..$ ImageWidth               : num 1013
  .. ..$ ImageLength              : chr "4952"
  .. ..$ Compression              : chr "1"
  .. ..$ PhotometricInterpretation: chr "2"
  .. ..$ ImageDescription         : chr "image description                "
  .. ..$ Orientation              : chr "1"
  .. ..$ SamplesPerPixel          : chr "3"
  .. ..$ XResolution              : chr "1000000/10000"
  .. ..$ YResolution              : chr "1000000/10000"
  .. ..$ ResolutionUnit           : chr "2"
  .. ..$ Software                 : chr "Adobe Photoshop CS5 Windows"
  .. ..$ DateTime                 : chr "2011-05-24T11:39:04Z"
  .. ..$ ColorSpace               : chr "1"
  .. ..$ MEDIAWIKI_EXIF_VERSION   : chr "1"
  .. ..$ GPSLongitude             : chr "5.6657"
  .. ..$ GPSLatitude              : chr "51.9871"
  .. ..$ ImageHeight              : num 4640
  .. ..$ GPSAltitude              : num 0
  .. ..$ ExposureTime             : chr ""
  .. ..$ FocalLength              : chr "50 mm"
  .. ..$ Flash                    : chr "No Flash"
  ..@ PhotoOverlay:List of 11
  .. ..$ rotation : num 0
  .. ..$ leftFov  : num -6.55
  .. ..$ rightFov : num 6.55
  .. ..$ bottomFov: num -30
  .. ..$ topFov   : num 30
  .. ..$ near     : num 50
  .. ..$ shape    : Factor w/ 1 level "rectangle": 1
  .. ..$ range    : num 1000
  .. ..$ tilt     : num 90
  .. ..$ heading  : num 0
  .. ..$ roll     : num 0
  ..@ sp          :Formal class 'SpatialPoints' [package "sp"] with 3 slots
  .. .. ..@ coords     : num [1, 1:3] 5.67 51.99 0
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL\
  .. .. .. .. ..$ : chr [1:3] "lon" "lat" "alt"
  .. .. ..@ bbox       : num [1:3, 1:2] 5.67 51.99 0 5.67 51.99 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:3] "lon" "lat" "alt"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

To visualize this object we again use the plotKML command:

> plotKML(sm)
KML file header opened for parsing...
Closing  sm.kml
Soil profile photograph with all EXIF information obtained via the WikiMedia API.

Soil profile data

Via the aqp package, we can construct SoilProfile-class objects that can also be visualized via plotKML as three dimensional monoliths. Because soil profiles are 3D samples of soil body, plotKML will attach elevation to different depths to be able to visualize soil profile data in Google Earth. This is an example:

> require(aqp)
Loading required package: aqp
This is aqp 1.5
> lon = 3.90; lat = 7.50; id = "ISRIC:NG0017"; FAO1988 = "LXp"
> top = c(0, 18, 36, 65, 87, 127)
> bottom = c(18, 36, 65, 87, 127, 181)
> ORCDRC = c(18.4, 4.4, 3.6, 3.6, 3.2, 1.2)
> hue = c("7.5YR", "7.5YR", "2.5YR", "5YR", "5YR", "10YR")
> value = c(3, 4, 5, 5, 5, 7); chroma = c(2, 4, 6, 8, 4, 3)
> prof1 <- join(data.frame(id, top, bottom, ORCDRC, hue, value, chroma), data.frame(id, lon, lat, FAO1988), type='inner')
Joining by: id
> prof1$soil_color <- with(prof1, munsell2rgb(hue, value, chroma))
Notice: converting hue to character
> depths(prof1) <- id ~ top + bottom
> site(prof1) <- ~ lon + lat + FAO1988
> coordinates(prof1) <- ~ lon + lat
> proj4string(prof1) <- CRS("+proj=longlat +datum=WGS84")
> prof1
Object of class SoilProfileCollection
Number of profiles: 1
Horizon attributes:
id top bottom ORCDRC   hue value chroma soil_color
1 ISRIC:NG0017   0     18   18.4 7.5YR     3      2  #584537FF
2 ISRIC:NG0017  18     36    4.4 7.5YR     4      4  #7E5A3BFF
3 ISRIC:NG0017  36     65    3.6 2.5YR     5      6  #A96C4FFF
4 ISRIC:NG0017  65     87    3.6   5YR     5      8  #B06A32FF
5 ISRIC:NG0017  87    127    3.2   5YR     5      4  #9A7359FF
6 ISRIC:NG0017 127    181    1.2  10YR     7      3  #C4AC8CFF
Sampling site attributes:
id FAO1988
1 ISRIC:NG0017     LXp
Spatial Data:
min max
lon 3.9 3.9
lat 7.5 7.5
CRS arguments:
 +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> plotKML(prof1, var.name="ORCDRC", color.name="soil_color")
KML file opened for writing...
Loading required package: fossil
Loading required package: maps
Writing KML...
Closing  prof1.kml
Soil profile description data shown as a block. The soil colors show the exact colors estimated from the Munsell color chart values.

Space-time classes (time series, trajectories)

Spatial classes from the sp package can be extended to spatio-temporal classes available via the package spacetime (Pebesma, 2012). The next two examples show how to visualize irregular point observations and trajectory-type data (cycling log). See also further examples with the RasterBrickTimeSeries class.

To visualize time series of measurements at meteorological stations you need to format the objects so they are of type STIDF (Space-time Irregular Data Frame, or unstructured spatio-temporal data) class:

> data(HRtemp08)
> HRtemp08$ctime <- as.POSIXct(HRtemp08$DATE, format="%Y-%m-%dT%H:%M:%SZ")
> library(spacetime)
> sp <- SpatialPoints(HRtemp08[,c("Lon","Lat")])
> proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")
> HRtemp08.st <- STIDF(sp, time = HRtemp08$ctime, data = HRtemp08[,c("NAME","TEMP")])
> HRtemp08_jan <- HRtemp08.st[1:500]
> str(HRtemp08_jan)
Formal class 'STIDF' [package "spacetime"] with 4 slots
  ..@ data   :'data.frame':     500 obs. of  2 variables:
  .. ..$ NA\ME: Factor w/ 158 levels "<c3><88>akovec",..: 17 50 81 83 84 82 85 86 88 89 ...
  .. ..$ TEMP: num [1:500] -3.5 -1.25 -4 NA\ 2.58 ...
  ..@ sp     :Formal class 'SpatialPoints' [package "sp"] with 3 slots
  .. .. ..@ coords     : num [1:500, 1:2] 17.2 15.6 17.4 15.5 13.6 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL\
  .. .. .. .. ..$ : chr [1:2] "Lon" "Lat"
  .. .. ..@ bbox       : num [1:2, 1:2] 13.6 42.4 19.4 46.4
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "Lon" "Lat"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
  ..@ time   :An ‘xts’ object on 2008-01-01 01:00:00/2008-01-04 01:00:00 containing:
  Data: int [1:500, 1] 298 433 916 1198 1683 2010 2526 2932 2966 3327 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL\
  ..$ : chr "timeIndex"
  Indexed by objects of class: [POSIXct,POSIXt] TZ:
  xts Attributes:
 NULL\
  ..@ endTime: POSIXct[1:500], format: "2008-01-01 01:00:00" "2008-01-01 01:00:00" ...
> plotKML(HRtemp08_jan[,,"TEMP"], dtime = 24*3600)
KML file header opened for parsing...
Parsing to KML...
Closing  HRtemp08_jan_,_,__TEMP__.kml
Time-series of meteorological measurements.

To visualize spatial trajectories, the data needs to be formatted in the STTDF (Spacetime Trajectory Data Frame) class:

> data(gpxbtour)
> gpxbtour$ctime <- as.POSIXct(gpxbtour$time, format="%Y-%m-%dT%H:%M:%SZ")
> coordinates(gpxbtour) <- ~lon+lat
> proj4string(gpxbtour) <- CRS("+proj=longlat +datum=WGS84")
> library(fossil)
> xy <- as.list(data.frame(t(coordinates(gpxbtour))))
> gpxbtour$dist.km <- sapply(xy, function(x) { deg.dist(long1=x[1], lat1=x[2], long2=xy[[1]][1], lat2=xy[[1]][2]) } )
> library(spacetime)
> library(adehabitat)
Loading required package: ade4
Attaching package: ‘ade4’
...
Attaching package: ‘adehabitat’
> gpx.ltraj <- as.ltraj(coordinates(gpxbtour), gpxbtour$ctime, id = "th")
> gpx.st <- methods/html/as.html">as(gpx.ltraj, "STTDF")
> gpx.st$speed <- gpxbtour$speed
> gpx.st@sp@proj4string <- CRS("+proj=longlat +datum=WGS84")
> str(gpx.st)
Formal class 'STTDF' [package "spacetime"] with 4 slots
..@ data:'data.frame':        3228 obs. of  13 variables:
.. ..$ x        : num [1:3228] 5.75 5.75 5.75 5.75 5.75 ...
.. ..$ y        : num [1:3228] 52 52 52 52 52 ...
.. ..$ date     : POSIXct[1:3228], format: "2010-06-25 10:35:29" "2010-06-25 10:35:34" "2010-06-25 10:35:37" ...
.. ..$ dx       : num [1:3228] 0.000466 0.000282 0.00031 0.000403 0.000309 ...
.. ..$ dy       : num [1:3228] -9.5e-05 -7.0e-05 -5.8e-05 -5.5e-05 -4.9e-05 ...
.. ..$ dist     : num [1:3228] 0.000476 0.000291 0.000315 0.000407 0.000313 ...
.. ..$ dt       : num [1:3228] 5 3 3 4 3 4 3 3 3 3 ...
.. ..$ R2n      : num [1:3228] 0.00 2.26e-07 5.87e-07 1.17e-06 2.21e-06 ...
.. ..$ abs.angle: num [1:3228] -0.201 -0.243 -0.185 -0.136 -0.157 ...
.. ..$ rel.angle: num [1:3228] NA\ -0.0422 0.0584 0.0493 -0.0216 ...
.. ..$ burst    : Factor w/ 1 level "th": 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ id       : Factor w/ 1 level "th": 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ speed    : num [1:3228] 23.6 23.1 24 24.9 26.1 ...
.. ..- attr(*, "id")= chr "th"
.. ..- attr(*, "burst")= chr "th"
..@ traj:List of 1
.. ..$ :Formal class 'STI' [package "spacetime"] with 2 slots
.. .. .. ..@ sp  :Formal class 'SpatialPoints' [package "sp"] with 3 slots
.. .. .. .. .. ..@ coords     : num [1:3228, 1:2] 5.75 5.75 5.75 5.75 5.75 ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : NULL\
.. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. .. .. ..@ bbox       : num [1:2, 1:2] 5.75 51.9 7.6 52.04
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. .. .. .. .. ..$ : chr [1:2] "min" "max"
.. .. .. .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. .. .. .. .. .. ..@ projargs: chr NA\
.. .. .. ..@ time:An ‘xts’ object from 2010-06-25 10:35:29 to 2010-06-25 16:49:56 containing:
Data: int [1:3228, 1] 1 2 3 4 5 6 7 8 9 10 ...
Indexed by objects of class: [POSIXct,POSIXt] TZ:
xts Attributes:
List of 3
.. .. .. .. ..$ tclass        : chr [1:2] "POSIXct" "POSIXt"
.. .. .. .. ..$ tzone         : chr ""
.. .. .. .. ..$ timeIsInterval: logi FALSE\
..@ sp  :Formal class 'SpatialPoints' [package "sp"] with 3 slots
.. .. ..@ coords     : num [1:2, 1:2] 5.75 7.6 51.9 52.04
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : NULL\
.. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
.. .. ..@ bbox       : num [1:2, 1:2] 5.75 51.9 7.6 52.04
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
.. .. .. .. ..$ : chr [1:2] "min" "max"
.. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. .. .. ..@ projargs: chr " +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
..@ time:An ‘xts’ object from 2010-06-25 10:35:29 to 2010-06-25 16:49:56 containing:
Data: int [1:2, 1] 1 2
Indexed by objects of class: [POSIXct,POSIXt] TZ:
xts Attributes:
List of 3
.. ..$ tclass        : chr [1:2] "POSIXct" "POSIXt"
.. ..$ tzone         : chr ""
.. ..$ timeIsInterval: logi FALSE\
> plotKML(gpx.st, colour="speed")
KML file header opened for parsing...
Parsing to KML...
Closing  gpx.st.kml
Movemement trajectory (cycling path) with vertices (GPS logs). In this case the cycling speed has been added as a color variable.

Spatial metadata

plotKML offers a set of simple functions such as spMetadata, that allow preparation, automated generation and export of spatial metadata. It uses the The Federal Geographic Data Committee (FGDC) metadata XML standards to save and parse the metadata into the KML files. These metadata can then be added to the plotKML plot via the metadata argument. This is an example with the Ebergötzen case study:

> eberg.md <- spMetadata(eberg, xml.file=system.file("eberg.xml", package="plotKML"), Target_variable="SNDMHT_A")
Reading the metadata file: C:/Program Files/R/R-2.15.3/library/plotKML/eberg.xml
Generating metadata...
Estimating the bounding box coordinates...
Reprojecting to +proj=longlat +datum=WGS84 ...
> str(eberg.md)
Formal class 'SpatialMetadata' [package "plotKML"] with 4 slots
  ..@ xml        :Classes 'XMLInternalDocument', 'XMLAbstractDocument', 'oldClass' <externalptr>
  ..@ field.names: chr [1:131] "Attribute_accuracy_report" "NA" "Completeness_Report" "Process_Date" ...
  ..@ palette    :Formal class 'sp.palette' [package "plotKML"] with 5 slots
  .. .. ..@ type  : chr "numeric"
  .. .. ..@ bounds: num [1:173] 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 ...
  .. .. ..@ color : chr [1:172] "#2C7BB6" "#2F7DB7" "#327FB8" "#3582BA" ...
  .. .. ..@ names : chr [1:172] "6.75" "7.25" "7.75" "8.25" ...
  .. .. ..@ icons : chr [1:172] "" "" "" "" ...
  ..@ sp         :Formal class 'Spatial' [package "sp"] with 2 slots
  .. .. ..@ bbox       : num [1:2, 1:2] 3569344 5707622 3580992 5718873
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "X" "Y"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. .. .. ..@ projargs: chr "+init=epsg:31467 +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +towgs84"| __truncated__
> plotKML(eberg["CLYMHT_A"], metadata=eberg.md)
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  eberg__CLYMHT_A__.kml
Adding some basic metadata to an object plot.

Writing analysis outputs to KML

plotKML can be used to only to visualize spatial and spatio-temporal classes, but also the results of complex spatial analysis. In this case the package will plot a series of layers and attach screen overlays, to allow direct interpretation of the analysis results. Such visualization is also implemented via specific S4 classes: RasterBrickTimeSeries, SpatialPredictions, SpatialSamplingPattern, RasterBrickSimulations, SpatialVectorsSimulations, SpatialMaxEntOutput and similar. To further customize visualizations consider combining the lower level functions kml_open, kml_layer, kml_close, kml_compress, kml_screen into your own wrapper function.

Time-series of images

Time series of images of the same variable can be imported into R as RasterBrickTimeSeries objects. This following example shows how to visualize a time series of MODIS Land Surface Temperature (LST) images. These first need to reprojected and converted to a RasterBrick class:

> require(raster)
> data(LST)
> gridded(LST) <- ~lon+lat
> proj4string(LST) <- CRS("+proj=longlat +datum=WGS84")
> dates <- sapply(strsplit(names(LST), "LST"), function(x){x[[2]]})
> datesf <- format(as.Date(dates, "%Y_%m_%d"), "%Y-%m-%dT%H:%M:%SZ")
> # begin / end dates +/- 4 days:
> TimeSpan.begin = as.POSIXct(unclass(as.POSIXct(datesf))-4*24*60*60, origin="1970-01-01")
> TimeSpan.end = as.POSIXct(unclass(as.POSIXct(datesf))+4*24*60*60, origin="1970-01-01")
> # pick few climatic stations:
> pnts <- HRtemp08[which(HRtemp08$NAME=="Pazin")[1],]
> pnts <- rbind(pnts, HRtemp08[which(HRtemp08$NAME=="Crni Lug - NP Risnjak")[1],])
> pnts <- rbind(pnts, HRtemp08[which(HRtemp08$NAME=="Cres")[1],])
> coordinates(pnts) <- ~Lon + Lat
> proj4string(pnts) <- CRS("+proj=longlat +datum=WGS84")
> # get the dates from the file names:
> LST_ll <- brick(LST[1:5])
> LST_ll@title = "Time series of MODIS Land Surface Temperature (8-day mosaics) images"

In the last step, we can construct an object of class RasterBrickTimeSeries that contains variable name, sample points, rasters and time period:

> LST.ts <- methods/html/new.html">new("RasterBrickTimeSeries", variable = "LST", sampled = pnts,
+     rasters = LST_ll, TimeSpan.begin = TimeSpan.begin[1:5], TimeSpan.end = TimeSpan.end[1:5])
> data(SAGA_pal)
> # plot MODIS images in Google Earth:
> plotKML(LST.ts, colour_scale=SAGA_pal[[1]])
KML file opened for writing...
Parsing to KML...
Writing to KML...
Closing  LST.ts.kml

The output plot allows us to visualize inspect changes in temperatures by scrolling the time slide-bar in the upper left corner. By clicking on sample points we can also see how the values in the LST images change for each date of interest.

Time series of MODIS LST images with a sample point showing the temporal changes at a meterological station.

Geostatistical mapping

The objective of geostatistical mapping is to generate spatial predictions of some target variable for a study area of interest. The outputs from geostatistical mapping are usually standard:

  • Observed values (point locations)
  • Geostatistical model (i.e. a list of fitted parameters)
  • Predicted values (gridded map)
  • Results of validation (point locations)

These can be presented in plotKML via the SpatialPredictions class. The following example (based on the GLM-kriging implemented in the GSIF package) shows how to quickly visualize results of geostatistical mapping via plotKML:

> data(meuse)
> coordinates(meuse) <- ~x+y
> proj4string(meuse) <- CRS("+init=epsg:28992")
> data(meuse.grid)
> gridded(meuse.grid) <- ~x+y
> proj4string(meuse.grid) <- CRS("+init=epsg:28992")
> library(GSIF)
GSIF version 0.3-0 (2013-03-07)
URL: http://gsif.r-forge.r-project.org/
> omm <- fit.gstatModel(observations = meuse, formulaString = om~dist, family = gaussian(log), covariates = meuse.grid)
Fitting a GLM...
Fitting a 2D variogram...
Saving an object of class 'gstatModel'...
> methods/html/show.html">show(omm@regModel)
Call:  glm(formula = om ~ dist, family = family, data = rmatrix)
 
Coefficients:
(Intercept)         dist
      2.390       -1.871
 
Degrees of Freedom: 152 Total (i.e. Null);  151 Residual
  (2 observations deleted due to missingness)
Null Deviance:      1791
Residual Deviance: 1076         AIC: 738.6
> om.rk <- predict(omm, predictionLocations = meuse.grid)
Generating predictions using the trend model (RK method)...
[using ordinary kriging]
 
 37% done
100% done
Running 5-fold cross validation using 'krige.cv'...
> methods/html/show.html">show(om.rk)
  Variable           : om
  Minium value       : 1
  Maximum value      : 17
  Size               : 153
  Total area         : 4964800
  Total area (units) : square-m
  Resolution (x)     : 40
  Resolution (y)     : 40
  Resolution (units) : m
  GLM call formula   : om ~ dist
  Family             : gaussian
  Link function      : log
  Vgm model          : Exp
  Nugget (residual)  : 4.24
  Sill (residual)    : 5.87
  Range (residual)   : 1080
  RMSE (validation)  : 2.323
  Var explained      : 53.9%
  Effective bytes    : 812
  Compression method : gzip
> plotKML(om.rk)
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  om.rk.kml
Plotting the whole geostatistical mapping project at once: sampling locations, predictions, cross-validation correlation plot, and variograms for residuals and cross-validation residuals.

Spatial Sampling patterns

Spatial samples can be generated via a number of packages in R (sp, spatstat, spcosa etc). In spcosa (Walvoort et al, 2010), for example, one can also optimize sampling design to maximize the geographical coverage of the sampling locations:

> library(spcosa)
> shpFarmsum <- readOGR(dsn = system.file("maps", package = "spcosa"), layer = "farmsum")
OGR data source with driver: ESRI Shapefile
Source: "C:/Program Files/R/R-2.15.3/library/spcosa/maps", layer: "farmsum"
with 1 features and 4 fields
Feature type: wkbPolygon with 2 dimensions
> # stratify `Farmsum' into 50 strata
> myStratification <- stratify(shpFarmsum, nStrata = 50)
> # sample two sampling units per stratum
> mySamplingPattern <- spsample(myStratification, n = 2)
> # attach the correct proj4 string:
> library(RCurl)
> nl.rd <- getURL("http://spatialreference.org/ref/sr-org/6781/proj4/")
> proj4string(mySamplingPattern@sample) <- CRS(nl.rd)
> # prepare spatial domain (polygons):
> sp.domain <- methods/html/as.html">as(myStratification@cells, "SpatialPolygons")
> sp.domain <- SpatialPolygonsDataFrame(sp.domain, data.frame(ID=as.factor(myStratification@stratumId)), match.ID = FALSE)
> proj4string(sp.domain) <- CRS(nl.rd)
> mySamplingPattern.ssp <- methods/html/new.html">new("SpatialSamplingPattern", method = class(mySamplingPattern), pattern = mySamplingPattern@sample, sp.domain = sp.domain)
> shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
> plotKML(mySamplingPattern.ssp, shape = shape)
Plotting the first variable on the list
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  mySamplingPattern.ssp.kml
Visualizing results of sampling optimization produced using the spcosa package.

Spatial simulations

In some cases, the result of some (geo)statistical model will be realizations i.e. a series of spatial objects with same intrinsic properties. The following two examples show how to visualize such data via plotKML i.e. by putting the outputs of analysis to the RasterBrickSimulations or SpatialVectorsSimulations class objects.

Consider for example the case with measured elevations for the Baranja hill case study. This data set contains 6370 precise observations of elevations, which can be used to construct a Digital Elevation Model (DEM). Instead of deriving just one smooth surface, we can model the land surface as a stochastic feature:

> data(barxyz)
> # define the projection system:
> prj = "+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +units=m +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824"
> coordinates(barxyz) <- ~x+y
> proj4string(barxyz) <- CRS(prj)
> data(bargrid)
> coordinates(bargrid) <- ~x+y
> gridded(bargrid) <- TRUE
> proj4string(bargrid) <- CRS(prj)
> # fit a variogram and generate simulations:
> Z.ovgm <- vgm(psill=1352, model="Mat", range=650, nugget=0, kappa=1.2)
> sel <- runif(length(barxyz$Z))<.2  # Note: this operation can be time consuming
> sims <- krige(Z~1, barxyz[sel,], bargrid, model=Z.ovgm, nmax=20, nsim=10, debug.level=-1)
drawing 10 GLS realisations of beta...
[using conditional Gaussian simulation]
 
  3% done
...
100% done
> # specify the cross-section:
> t1 <- Line(matrix(c(bargrid@bbox[1,1],bargrid@bbox[1,2],5073012,5073012), ncol=2))
> transect <- SpatialLines(list(Lines(list(t1), ID="t")), CRS(prj))
> # glue to a RasterBrickSimulations object:
> bardem_sims <- methods/html/new.html">new("RasterBrickSimulations", variable = "elevations",
+ sampled = transect, realizations = brick(sims))
> # plot the whole project and open in Google Earth:
> data(R_pal)
> plotKML(bardem_sims, colour_scale = R_pal[[4]])
KML file opened for writing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Reprojecting to +proj=longlat +datum=WGS84 ...
Writing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  bardem_sims.kml

The resulting plot can be used to inspect the uncertainty in the multiple equiprobable DEMs. Note also that, by default, plotKML adds a cross-section from the mid point West-East, that allows us to see the probability distributions in vertical dimension.

Multiple equiprobable geostatistical simulations of land surface and associated uncertainty.

Likewise, plotKML can also be used to visualize multiple realizations of vector features via the SpatialVectorsSimulations class. Here is an example with multiple stream networks derived using geostatistical simulations of DEMs, shown previously (see Hengl et al, 2010 for more details):

> data(barstr)
> data(bargrid)
> coordinates(bargrid) <- ~ x+y
> gridded(bargrid) <- TRUE
> cell.size = bargrid@grid@cellsize[1]
> bbox = bargrid@bbox
> gridT = GridTopology(cellcentre.offset=bbox[,1], cellsize=c(cell.size,cell.size),
+    cells.dim=c(round(abs(diff(bbox[1,])/cell.size), 0),
+    ncols=round(abs(diff(bbox[2,])/cell.size), 0)))
> bar_sum <- aggregate(gridT, barstr[1:5])
> plotKML(bar_sum, grid2poly = TRUE)
Warning in grid2poly(obs) :
Operation not recommended for large grids (>>1e4 pixels).
Reprojecting to +proj=longlat +datum=WGS84 ...
KML file header opened for parsing...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
...
Closing  bar_sum.kml
Multiple realizations of equiprobable vector-type features (streams) and the cumulative probability from blue (0) to red (1) colours (gridded).

Species distribution modeling

A popular method to run the species distribution model to estimate potential distribution of a species is the MaxEnt software. Via the package dismo (Hijmans and Elith, 2012), one can run the modelling in R, and then visualize results in Google Earth by using:

> data(bigfoot)
> aea.prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
> coordinates(bigfoot) <- ~Lon+Lat
> proj4string(bigfoot) <- CRS("+proj=latlon +datum=WGS84")
> bigfoot.aea <- spTransform(bigfoot, CRS(aea.prj))
> data(USAWgrids)
> gridded(USAWgrids) <- ~s1+s2
> proj4string(USAWgrids) <- CRS(aea.prj)
> bigfoot.aea <- as.ppp(spTransform(bigfoot, CRS(aea.prj)))
> sel.grids <- c("globedem","nlights03","sdroads","gcarb","twi","globcov")
> library(GSIF)
> jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep='')
> if(file.exists(jar)){
+   bigfoot.smo <- MaxEnt(bigfoot.aea, USAWgrids[sel.grids])
+   icon = "http://plotkml.r-forge.r-project.org/bigfoot.png"
+   data(R_pal)
+   plotKML(bigfoot.smo, colour_scale = R_pal[["bpy_colors"]], shape = icon)
+ }
KML file header opened for parsing...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Closing  bigfoot.smo.kml
Mapping distribution of the bigfoot using the BigFoot Research Organization (BFRO) data.

Note that the result are two maps: one showing the 'ecological' probability of bigfoot's niche, and the second map showing the predicted domain (with values 1 or missing) i.e. areas where the occurrence of the species based on cross-validation is significant.

Working with large data sets

Visualization of large datasets via plotKML and kml methods is not recommended because R is relatively inefficient for working with large rasters. Two robust procedures to export larger areas i.e. larger data sets to KML are (a) by simplifying the spatial structure or raster resolution, and (b) by tiling large grids or vector maps to regular blocks (see tile function in the GSIF package) and then writing them in a loop. The following two sections provide some general guidelines on how to visualize large grids.

Tiling grids and vectors

Consider for example the 5 km resolution grids for the contiguous 48-state part of USA. We can download the grids and load them to R by using:

> download.file("http://worldgrids.org/rda/USgrids.rda", "USgrids.rda")
trying URL 'http://worldgrids.org/rda/USgrids.rda'
Content type 'text/plain' length 7952764 bytes (7.6 Mb)
opened URL
downloaded 7.6 Mb
> load("USgrids.rda")
> str(USgrids)
'data.frame':   312366 obs. of  23 variables:
 $ mukey   : Factor w/ 9570 levels "657753","657754",..: 8568 8568 8568 8564 8568 8568 8568 8568 8529 8569 ...
 $ suborder: Factor w/ 69 levels "","Albolls","Aqualfs",..: 39 39 39 65 39 39 39 39 66 39 ...
 $ pct     : num  61 61 61 57 61 61 61 61 77 60 ...
 $ anthroms: Factor w/ 21 levels "0","11","12",..: 9 10 9 9 10 10 3 3 3 10 ...
 $ biocl15 : int  46 45 45 45 44 43 43 42 42 45 ...
 $ biocl5  : num  22.2 22.3 22.1 22.2 22.3 ...
 $ biocl6  : num  -0.588 -0.769 -0.205 -0.378 -0.572 ...
 $ globcov : Factor w/ 22 levels "14","20","30",..: 5 5 5 5 5 5 5 5 3 5 ...
 $ globedem: int  31 58 11 22 33 37 37 46 44 40 ...
 $ LAIm    : int  97 99 62 97 100 102 100 99 100 92 ...
 $ LSTDm   : num  12.9 13.4 12.1 12.9 13.4 13.8 14.2 14.5 14.4 11.9 ...
 $ LSTDs   : num  418 436 375 418 436 ...
 $ LSTNm   : num  3.9 3.6 3.8 3.4 3.3 3.3 3.2 3.2 3.3 4.3 ...
 $ LSTNs   : num  317 322 305 311 316 ...
 $ PCEVI1  : num  4275 4670 NA\ 4484 4655 ...
 $ PCEVI2  : num  1018 1077 NA\ 1132 979 ...
 $ PCEVI3  : num  -180 -185 NA\ -272 -144 ...
 $ PCEVI4  : num  -63 -81.9 NA\ -86.4 -65.3 -36.4 4 56.5 55.1 -59.5 ...
 $ PRECm   : int  33 35 31 33 35 36 38 41 43 31 ...
 $ slope   : num  0.00351 0.0032 0.00303 0.0033 0.00336 ...
 $ wildness: int  0 0 0 0 0 0 0 0 0 0 ...
 $ x       : num  -1952500 -1947500 -1957500 -1952500 -1947500 ...
 $ y       : num  3162500 3162500 3157500 3157500 3157500 ...
> gridded(USgrids) <- ~x+y
> aea.prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
> proj4string(USgrids) = aea.prj
> zlims = range(USgrids@data[,"LSTDm"], na.rm=TRUE)
> zlims
[1] -2.8 40.4

This is a relatively large data set with 312,366 pixels. To visualize these grids using plotKML method would be very inefficient. Instead, we can first tile one layer in the grids into 5 decimal degree tiles. First, we need to reproject the layer of interest (LSTDm: MODIS-estimated Land Surface Temperature day time):

> library(GSIF)
> LSTDm <- reproject(USgrids["LSTDm"], program="FWTools")
> LSTDm <- methods/html/as.html">as(LSTDm, "SpatialPixelsDataFrame")
> LSTDm <- tile(USgrids["LSTDm"], block.x=5e5)
Generating 78 tiles...
Returning a list of tiles with 0 percent overlap
Subseting object of class "SpatialPixelsDataFrame"...

This resulted in a list of 78 tiles. To plot this list we use:

> plotKML(LSTDm.l, z.lim = zlims, colour_scale=SAGA_pal[[1]])
KML file opened for writing...
Writing to KML...
Writing to KML...
...
Tiling large grids using tile function from the GSIF package.

The whole process can still take few minutes, but this is still more efficient than if we would run:

> plotKML(USgrids["LSTDm"], colour_scale=SAGA_pal[[1]])

Important: when working with objects of type SpatialLinesDataFrame, SpatialPolygonsDataFrame and or RasterLayer, the tile function looks for FWTools binary files ogr2ogr and gdalwarp. FWTools is a separate program and must be installed separately otherwise the program will not run.

Writing objects in loops

Another possibility to write large grids / large vectors to KML files is the use loops. For example, to write a polygon map using loops we would run something like this:

> data(eberg_zones)
> ## this one requires ogr2ogr function:
> pol.lst <- tile(eberg_zones, block.x=5000)
Generating 4 tiles...
Returning a list of tiles with 0 percent overlap
Clipling lines using 'ogr2ogr'...
> kml_open("pol_lst.kml")
KML file opened for writing...
> for(i in 1:length(pol.lst)){
+   kml_layer.SpatialPolygons(obj=pol.lst[[i]], colour = ZONES)
+ }
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
Reprojecting to +proj=longlat +datum=WGS84 ...
Parsing to KML...
> kml_close("pol_lst.kml")
Closing  pol_lst.kml

In fact, the lower level functions kml_open, kml_layer, kml_close, kml_compress, kml_screen can be combined into new functions, so that the whole process can be automated. For more info refer to the plotKML manual.

References

  1. Bivand, R.S., Pebesma, E.J., and Gómez-Rubio, V., (2008) Applied Spatial Data Analysis with R. Springer, 378 p.
  2. Böhner, J., McCloy, K. R. and Strobl, J. (Eds), (2006) SAGA — Analysis and Modelling Applications. Göttinger Geographische Abhandlungen, Heft 115. Verlag Erich Goltze GmbH, Göttingen, 117 pp.
  3. Hengl, T., Heuvelink, G. B. M., van Loon, E. E., (2010) On the uncertainty of stream networks derived from elevation data: the error propagation approach. Hydrology and Earth System Sciences, 14:1153-1165.
  4. Hijmans, R.J., Elith, J., (2012) Species distribution modeling with R. CRAN, Vignette for the dismo package, 72 p.
  5. Hijmans, R.J., (2012) Introduction to the 'raster' package. Contributed R Achive Network (CRAN), p. 26.
  6. Pebesma, E., (2012) Classes and Methods for Spatio-Temporal Data in R. Journal of Statistical Software, 51(7): 30.
  7. Walvoort, D.J.J, Brus, D.J., de Gruijter, J.J., (2010) An R package for spatial coverage sampling and random sampling from compact geographical strata by k-means. Computers & Geosciences, 36(10): 1261-1267.
  8. Wernecke, J., (2009) The KML Hanbook: Geographic Visualization for the Web. Addison-Wesley, p. 331.

wiki tutorial plotkml
 stars  from 1 votes (Details)