Difference between revisions of "Export maps to GE"

(Replacing page with '{{#meta: REFRESH | 1;url=http://plotkml.r-forge.r-project.org/tutorial.php}} #REDIRECTS to [http://plotkml.r-forge.r-project.org/tutorial.php plotKML tutorial]')
 
Line 1: Line 1:
{{Drop|T}}his article gives some basic instructions to export various types of GIS layers to Google Earth. It also includes some basic instructions to resample and reproject maps, including steps to read back KML files to R. To obtain and install software used in this exercise, please go to section [[Software|software]]. The complete script used to generate the plots shown down-below, is available [http://spatial-analyst.net/scripts/KML_export.R here].
+
{{#meta: REFRESH | 1;url=http://plotkml.r-forge.r-project.org/tutorial.php}}
  
All displays in '''Google Earth''' are controlled by the KML files, which are written in the [http://earth.google.com/kml/kml_tut.html Keyhole Markup Language] developed by Keyhole, Inc. KML is an XML-based language for managing the display of three-dimensional geospatial data and is, in fact, used in several geographical browsers (Google Maps, Google Mobile, ArcGIS Explorer and World Wind). The KML file specifies a set of standard features (placemarks, images, polygons, 3D models, textual descriptions, etc.) for display in Google Earth. Each place always has a longitude and a latitude. Other data can make the view more specific, such as tilt, heading, altitude, which together define a camera view. The KML datasets can be edited using an ASCII editor (as with HTML), but they can be edited also directly in Google Earth. KML files are very often distributed as KMZ files, which are zipped KML files with a *.kmz extension. This [http://kml-samples.googlecode.com/svn/trunk/interactive/index.html interactive KML sampler] will give you some good idea about what is all possible in Google Earth.
+
#REDIRECTS to [http://plotkml.r-forge.r-project.org/tutorial.php plotKML tutorial]
 
 
To install Google Earth, run the Google Desktop installer that you can obtain from the [http://earth.google.com/download-earth.html 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. To load your own GIS data to Google Earth, there are few possibilities. 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 [http://earth-info.nga.mil/GandG/wgs84/ WGS84] ellipsoid), all vector/raster maps need to be first reprojected before they can be exported to KML format. More about importing the data to Google Earth can be found via the [http://earth.google.com/userguide/v4/ug_importdata.html Google Earth User Guide].
 
 
 
{{Plainlist|align=right|background=lightyellow|heading=Sample script to exchange maps to KML format.|
 
* [[Image:icon_R.png|18px]] [http://spatial-analyst.net/scripts/KML_export.R KML_export.R] : R script to export maps to KML format.
 
}}
 
 
 
 
 
 
 
 
 
= Exporting vector maps =
 
 
 
Export of vector layers from a GIS to KML is commonly straight forward. This is mainly because the structure of point or line features is rather simple, and each spatial element carries spatial coordinates without any topological relations. Export of polygon maps can be somewhat more complicated, especially if the polygons overlap or if they describe 3D features.
 
 
 
== Point and line features  ==
 
 
 
To export point or line features to KML, use the '''writeOGR''' method that is available in R package rgdal. Export can be achieved in three steps, e.g.:
 
 
 
1. Load the rgdal package for GIS data exchange:
 
 
 
<pre>&gt; require(c("rgdal","gstat","lattice","RASAGA","maptools","akima")) </pre>
 
 
 
2. Reproject the original map from local coordinates (in this case [http://spatialreference.org/ref/epsg/28992/ '''Amersfoort / RD New''']) to the longlat system (where necessary):
 
 
 
<pre>
 
> data(meuse)
 
> NL_RD <- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
 
+  +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel
 
+  +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812
 
+  +units=m +no_defs"
 
> coordinates(meuse) <- ~x+y
 
> proj4string(meuse) <- CRS(NL_RD)
 
> meuse.ll <- spTransform(meuse, CRS("+proj=longlat +datum=WGS84"))
 
</pre>
 
 
 
3. Export the point map using the "KML" OGR driver:
 
 
 
<pre>
 
> writeOGR(meuse.ll["lead"], "meuse_lead.kml", "lead", "KML")
 
</pre>
 
 
A more sophisticated way to generate a KML is to directly write to a KML file using loops. This way we have a full control of the visualization parameters. For example, to produce a bubble-type of plot in Google Earth with actual numbers attached as labels to a point map, we can do:
 
 
 
<pre>
 
varname <- "lead"  # variable name
 
maxvar <- max(meuse.ll[varname]@data)  # maximum value
 
filename <- file(paste(varname, "_bubble.kml", sep=""), "w",  blocking=FALSE)
 
write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>", filename)
 
write("<kml xmlns=\"http://earth.google.com/kml/2.2\">", filename, append = TRUE)
 
write("<Document>", filename, append = TRUE)
 
write(paste("<name>", varname, "</name>", sep=" "), filename, append = TRUE)
 
write("<open>1</open>", filename, append = TRUE)
 
for (i in 1:length(meuse.ll@data[[1]])) {
 
  write(paste('  <Style id="','pnt', i,'">',sep=""), filename, append = TRUE)
 
  write("    <LabelStyle>", filename, append = TRUE)     
 
  write("    <scale>0.7</scale>", filename, append = TRUE)
 
  write("    </LabelStyle>", filename, append = TRUE)
 
  write(" <IconStyle>", filename, append = TRUE)
 
  write("  <color>ff0000ff</color>", filename, append = TRUE)
 
#  write("  <colorMode>random</colorMode>", filename, append = TRUE)
 
  write(paste(" <scale>", meuse.ll[i,varname]@data[[1]]/maxvar*2+0.3, "</scale>", sep=""),
 
+    filename, append = TRUE)
 
  write(" <Icon>", filename, append = TRUE)
 
  write(" <href>http://maps.google.com/mapfiles/kml/shapes/donut.png</href>",
 
+    filename, append = TRUE)
 
  write(" </Icon>", filename, append = TRUE)
 
  write(" </IconStyle>", filename, append = TRUE)
 
  write(" </Style>", filename, append = TRUE)
 
}
 
write("<Folder>", filename, append = TRUE)
 
write(paste("<name>Donut icon for", varname,"</name>"), filename, append = TRUE)
 
for (i in 1:length(meuse.ll@data[[1]])) {
 
  write("  <Placemark>", filename, append = TRUE)
 
  write(paste("  <name>", meuse.ll[i,varname]@data[[1]],"</name>", sep=""), filename, append = TRUE)
 
  write(paste("  <styleUrl>#pnt",i,"</styleUrl>", sep=""), filename, append = TRUE)
 
  write("    <Point>", filename, append = TRUE)
 
  write(paste("      <coordinates>",coordinates(meuse.ll)[[i,1]],",",
 
+        coordinates(meuse.ll)[[i,2]],",10</coordinates>", sep=""), filename, append = TRUE)
 
  write("    </Point>", filename, append = TRUE)
 
  write("    </Placemark>", filename, append = TRUE)
 
}
 
write("</Folder>", filename, append = TRUE)
 
write("</Document>", filename, append = TRUE)
 
write("</kml>", filename, append = TRUE)
 
close(filename)
 
</pre>
 
 
This will produce the following plot:
 
 
 
[[Image:Fig meuse bubble plot GE.jpg|350px|Bubble plot in Google Earth.]]
 
 
 
== Ploting 3D points ==
 
 
 
If measurement locations refer to some specific altitude, we can also plot them as 3D points. This is an excellent opportunity to visualize e.g. soil profile data, fish acoustics data and similar.
 
 
 
<pre>
 
# we do not really have 3D data but we can simulate them:
 
> meuse.3D <- data.frame(meuse.ll)
 
# str(meuse.3D)
 
> meuse.3D <- rbind(meuse.3D[,c("x","y","lead")], data.frame(x=meuse.3D$x, y=meuse.3D$y, lead=rep(0, length(meuse.ll$lead))))
 
</pre>
 
 
 
Now we can plot the values of lead as 3D points (in this case we attach values of lead as heights):
 
 
 
<pre>
 
> varname <- "lead"
 
> filename <- file(paste(varname, "_3D.kml", sep=""), "w", blocking=FALSE)
 
> maxvar <- max(log1p(meuse.3D[,varname]), na.rm=TRUE)  # maximum value
 
> write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>", filename)
 
> write("<kml xmlns=\"http://earth.google.com/kml/2.2\">", filename, append = TRUE)
 
> write("<Document>", filename, append = TRUE)
 
> write(paste("<name>", varname, "in ppm</name>"), filename, append = TRUE)
 
> write("<open>1</open>", filename, append = TRUE)
 
> for (i in 1:length(meuse.3D[[1]])) {
 
+  write(paste('  <Style id="','pnt', i,'">',sep=""), filename, append = TRUE)
 
+  write("    <LabelStyle>", filename, append = TRUE) 
 
+  write("    <scale>0.7</scale>", filename, append = TRUE)
 
+  write("    </LabelStyle>", filename, append = TRUE)
 
+  write(" <IconStyle>", filename, append = TRUE)
 
+  write("  <color>ff0008cf</color>", filename, append = TRUE) # obtain the color converter at www.sgrillo.net
 
# write("  <colorMode>random</colorMode>", filename, append = TRUE)
 
+  write(paste(" <scale>", log1p(meuse.3D[i,varname])/maxvar*1.2+0.1,  "</scale>", sep=""), filename, append = TRUE)
 
+  write(" <Icon>", filename, append = TRUE)
 
+  write(" <href>http://maps.google.com/mapfiles/kml/shapes/donut.png</href>",  filename, append = TRUE)
 
+  write(" </Icon>", filename, append = TRUE)
 
+  write(" </IconStyle>", filename, append = TRUE)
 
+  write(" </Style>", filename, append = TRUE)
 
+  }
 
> write("<Folder>", filename, append = TRUE)
 
> write(paste("<name>", i,"</name>", sep=""), filename, append = TRUE)
 
> for (i in 1:length(meuse.3D[[1]])) {
 
+  Lon <- meuse.3D$x[i]
 
+  Lat <- meuse.3D$y[i]
 
+  Depth <- meuse.3D[i,varname]/2
 
+  write("  <Placemark>", filename, append = TRUE)
 
#  write(paste("  <name>", round(meuse.3D[i,varname], 0), "</name>", sep=""), filename, append = TRUE)
 
+  write(paste("  <styleUrl>#pnt",i,"</styleUrl>", sep=""), filename, append = TRUE)
 
+  write("    <Point>", filename, append = TRUE)
 
+  write("    <extrude>1</extrude>", filename, append = TRUE)
 
+  write("    <altitudeMode>absolute</altitudeMode>", filename, append = TRUE)
 
+  write(paste("      <coordinates>", Lon, ",", Lat, ",", Depth, "</coordinates>", sep=""), filename, append = TRUE)
 
+  write("    </Point>", filename, append = TRUE)
 
+  write("    </Placemark>", filename, append = TRUE)
 
+ }
 
> write("</Folder>", filename, append = TRUE)
 
> write("</Document>", filename, append = TRUE)
 
> write("</kml>", filename, append = TRUE)
 
> close(filename)
 
</pre>
 
 
 
[[Image:Fig_ge_meuse_lead_3D.jpg|350px|Plotting 3D points in Google Earth.]]
 
 
 
== Polygon maps ==
 
 
 
Polygon maps can also be exported using the writeOGR command, as implemented in the package '''rgdal'''. In the case of the meuse dataset, we first need to prepare a polygon map:
 
 
 
<pre>
 
> data(meuse.grid)
 
> coordinates(meuse.grid) <- ~x+y
 
> gridded(meuse.grid) <- TRUE
 
> proj4string(meuse.grid) = CRS("+init=epsg:28992")
 
# raster to polygon conversion;
 
> writeGDAL(meuse.grid["soil"], "meuse_soil.sdat", "SAGA", mvFlag=-99999)
 
> rsaga.geoprocessor(lib="shapes_grid", module=6,
 
+  param=list(GRID="meuse_soil.sgrd", POLYGONS="meuse_soil.shp", CLASS_ALL=1))
 
> soil <- readShapePoly("meuse_soil.shp", proj4string=CRS(NL_RD), force_ring=T)
 
> soil.ll <- spTransform(soil, CRS("+proj=longlat +datum=WGS84"))
 
> writeOGR(soil.ll["NAME"], "meuse_soil.kml", "soil", "KML")
 
</pre>
 
 
 
The actual plot in Google Earth can later be adjusted by manually editing the colors:
 
 
 
[[Image:Fig_meuse_polygon_plot_GE.jpg|350px|Exported polygon map.]]
 
 
 
 
 
= Exporting raster maps =
 
 
 
Rasters can not be exported that easy to KML. Google Earth does not allow import of GIS raster formats, but only input of images that can then be draped over a terrain ([http://code.google.com/apis/kml/documentation/kml_tut.html#ground_overlays ground overlay]). The images need to be exclusively in one of the following formats: JPG, BMP, GIF, TIFF, TGA and PNG. Typically, export of raster maps to KML consists of two steps: (1) reprojection/resampling (if the coordinate system of the raster map is different from the LatLonWGS84 system), and (2) export to KML.
 
 
 
== Reprojecting and resampling maps ==
 
 
 
1. Load the data to R and attach the correct coordinate system:
 
 
 
<pre>
 
> data(meuse.grid)
 
> coordinates(meuse.grid) <- ~x+y
 
> gridded(meuse.grid) <- TRUE
 
> proj4string(meuse.grid) <- CRS(NL_RD)
 
</pre>
 
 
 
2. Reproject the map and determine the grid system of the map in the LatLonWGS84 system. You need to determine five parameters: southern edge, northern edge, western edge, eastern edge and grid cell size in arcdegrees.
 
 
 
<pre>
 
> meuse.grid.ll <- spTransform(meuse.grid, CRS("+proj=longlat +datum=WGS84"))
 
# this makes a point map!
 
> corrf <- (1+cos((meuse.grid.ll@bbox[2,2]+meuse.grid.ll@bbox[2,1])/2*pi/180))/2
 
> geogrd.cell <- corrf*(meuse.grid.ll@bbox[1,2]-meuse.grid.ll@bbox[1,1])/meuse.grid@grid@cells.dim[1]
 
> geoarc <- spsample(meuse.grid.ll, type="regular", cellsize=c(geogrd.cell,geogrd.cell))
 
> gridded(geoarc) <- TRUE
 
> gridparameters(geoarc)
 
  cellcentre.offset    cellsize cells.dim
 
x1          5.721994 0.0004603493        95
 
x2        50.957021 0.0004603493        80
 
> gridparameters(meuse.grid)
 
  cellcentre.offset cellsize cells.dim
 
x            178460      40        78
 
y            329620      40      104
 
</pre>
 
 
 
Note that dimension of the two maps are not the same (they will almost never be the same), but we have at least about the same number of pixels in both maps. The grid cell size in arcdegrees is estimated using the following formula ([http://dx.doi.org/10.1016/S0166-2481(08)00002-0 Hengl and Evans, 2008]):
 
 
 
 
 
<tex>
 
\Delta x_{\mathit{m}} = F \cdot \cos(\varphi) \cdot \Delta x_{\mathit{m}}^0
 
</tex>
 
 
 
 
 
where <tex>\delta x_{\mathit{m}}</tex> is the East/West grid cell size estimated for a given latitude <tex>\varphi</tex>, <tex>\Delta x_{\mathit{m}}^0</tex> is the grid cell size in arcdegrees at equator, <tex>F=111,319</tex> is an empirical constant used to convert from degrees to meters.
 
 
 
Reprojection of gridded map is also available via the [http://cran.r-project.org/web/packages/raster/ raster] package, in which case we only need to specify target coordinate system and the package will determine the cell size / bounding box for us:
 
 
 
<pre>
 
> library(raster)
 
> dist.r <- raster(meuse.grid["dist"])
 
> dist.r.ll <- projectRaster(dist.r, crs="+proj=longlat +datum=WGS84")
 
</pre>
 
 
 
We can look at the maps just to see that everything is fine:
 
 
 
<pre>
 
> par(mfrow=c(1,2))
 
> image(dist.r)
 
> image(dist.r.ll)
 
</pre>
 
 
 
[[Image:Fig_meuse_resample.png|450px|Reprojected and resampled meuse distance to river map.]]
 
 
 
To finally export the map (either the map or the actual plot) to KML ground overlay, we can use:
 
 
 
<pre>
 
> meuse.kml <- GE_SpatialGrid(dist.grid.ll)
 
> col.no <- 48 # number of colors;
 
> col.pal <-  c(grey(0), rev(rainbow(65)[1:col.no]), rep(grey(1), 256-col.no))
 
> z.lim <- c(quantile(dist.grid.ll@data[,1], 0.025, na.rm=TRUE), quantile(dist.grid.ll@data[,1], 0.975, na.rm=TRUE))
 
> breaks.var <- c(-99999, seq(from=z.lim[1], to=z.lim[2], length.out=col.no-1), max(dist.grid.ll@data[,1], na.rm=TRUE)+sd(dist.grid.ll@data[,1],
 
+      na.rm=TRUE)/col.no)> dist.grid.ll$img <- as.integer(cut(dist.grid.ll@data[,1], breaks=breaks.var, include.lowest=TRUE))
 
# write a GIF:
 
> library(caTools)
 
> write.gif(image=t(as.matrix(dist.grid.ll["img"])), filename="meuse_dist_ll.gif", col=col.pal, transparent=0)
 
> kmlOverlay(meuse.kml, kmlfile="meuse_dist.kml", imagefile="meuse_dist_ll.gif", name="Distance from the river")
 
[1] "<?xml version='1.0' encoding='UTF-8'?>"                                                                                                     
 
[2] "<kml xmlns='http://earth.google.com/kml/2.0'>"                                                                                               
 
[3] "<GroundOverlay>"                                                                                                                             
 
[4] "<name>Distance from the river</name>"                                                                                                       
 
[5] "<Icon><href>meuse_dist_ll.gif</href><viewBoundScale>0.75</viewBoundScale></Icon>"                                                           
 
[6] "<LatLonBox><north>50.9958956945964</north><south>50.9530556945964</south><east>5.76863376192979</east><west>5.717946932831</west></LatLonBox>"
 
[7] "</GroundOverlay></kml>"
 
</pre>
 
 
 
Which would also allow us to share this load this data to another GIS (note the size of the map is rather large because it is uncompressed). We could also instead export the actual R plot:
 
 
 
<pre>
 
> tf <- tempfile(tmpdir=getwd())
 
> png(file="meuse.png", width=meuse.kml$width, height=meuse.kml$height, bg="transparent")
 
> par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
 
> image(dist.r.ll, col=bpy.colors(), xlim=meuse.kml$xlim, ylim=meuse.kml$ylim)
 
> kmlOverlay(meuse.kml, kmlfile="meuse_dist_plot.kml", imagefile="meuse.png", name="Distance from the river (plot)")
 
[1] "<?xml version='1.0' encoding='UTF-8'?>"                                                                                                     
 
[2] "<kml xmlns='http://earth.google.com/kml/2.0'>"                                                                                               
 
[3] "<GroundOverlay>"                                                                                                                             
 
[4] "<name>Distance from the river (plot)</name>"                                                                                                 
 
[5] "<Icon><href>meuse.png</href><viewBoundScale>0.75</viewBoundScale></Icon>"                                                                   
 
[6] "<LatLonBox><north>50.9958956945964</north><south>50.9530556945964</south><east>5.76863376192979</east><west>5.717946932831</west></LatLonBox>"
 
[7] "</GroundOverlay></kml>"                                                                                                                     
 
> dev.off()
 
null device
 
          1
 
</pre>
 
 
 
This would produce an plot like this:
 
 
 
[[Image:Fig_GE_meuse.jpg|450px|Google Earth visualization of an R plot.]]
 
 
 
== Ground overlays ==
 
 
 
A ground overlay (with a legend) has the following structure:
 
 
 
<pre>
 
<?xml version="1.0" encoding="UTF-8"?>
 
<kml xmlns="http://earth.google.com/kml/2.1">
 
<Document>
 
    <name>Raster map example (meuse)</name>
 
    <GroundOverlay>
 
        <name>Map name</name>
 
        <description>Description of how was map produced.</description>
 
        <Icon>
 
                <href>meuse_dist.jpg</href>
 
            </Icon>
 
            <LatLonBox>
 
                <north>50.99362</north>
 
                <south>50.95679</south>
 
                <east>5.765532</east>
 
                <west>5.721764</west>
 
            </LatLonBox>
 
        </GroundOverlay>
 
        <ScreenOverlay>
 
            <name>Legend</name>
 
            <Icon>
 
                <href>meuse_dist_legend.jpg</href>
 
            </Icon>
 
            <overlayXY x="0" y="1" xunits="fraction" yunits="fraction"/>
 
            <screenXY x="0" y="1" xunits="fraction" yunits="fraction"/>
 
            <rotationXY x="0" y="0" xunits="fraction" yunits="fraction"/>
 
            <size x="0" y="0" xunits="fraction" yunits="fraction"/>
 
        </ScreenOverlay>
 
</Document>
 
</kml>
 
</pre>
 
 
 
== Rescaling maps in SAGA GIS  ==
 
 
 
This section shows some basic instructions on how to import data into a GIS and how to automatically resample data layers to some arbitrary grid cell size. NOTE: aggregation of maps (upscaling) to general scales is often only a technical issue (unless you want to aggregate classes such as land cover); disaggregation of maps (downscaling) is more tricky and can lead to artifacts and decrease in map quality.
 
 
 
The fastest way to re-scale gridded maps is to use the grid_tools library of SAGA GIS:
 
<pre>&gt; library(RSAGA)
 
&gt; rsaga.get.usage("grid_tools", 0)
 
SAGA CMD 2.0.3
 
library path:  C:/Progra~1/saga_vc/modules
 
library name:  grid_tools
 
module name&nbsp;:  Resampling
 
Usage: 0 -INPUT &lt;str&gt; [-GRID &lt;str&gt;] [-METHOD &lt;num&gt;] [-KEEP_TYPE] ...
 
  -INPUT:&lt;str&gt;                  Grid
 
        Grid (input)
 
  -GRID:&lt;str&gt;                  Resampled Grid
 
        Data Object (optional output)
 
  -METHOD:&lt;num&gt;                Target Grid
 
        Choice
 
        Available Choices:
 
        [0] Specify dimensions
 
        [1] Create new grid in existing project
 
        [2] Overwrite existing grid
 
  -KEEP_TYPE                    Preserve Data Type
 
        Boolean
 
  -DIMENSIONS_CELLSIZE:&lt;str&gt;    Cell Size
 
        Floating point
 
        Minimum: 0.000000
 
...
 
  -GRID_GRID:&lt;str&gt;              Target Grid
 
        Grid (input)
 
  -SCALE_UP_METHOD:&lt;num&gt;        Interpolation Method
 
        Choice
 
        Available Choices:
 
        [0] Nearest Neighbor
 
        [1] Bilinear Interpolation
 
        [2] Inverse Distance Interpolation
 
        [3] Bicubic Spline Interpolation
 
        [4] B-Spline Interpolation
 
        [5] Mean Value
 
  -SCALE_DOWN_METHOD:&lt;num&gt;      Interpolation Method
 
        Choice
 
        Available Choices:
 
        [0] Nearest Neighbor
 
        [1] Bilinear Interpolation
 
        [2] Inverse Distance Interpolation
 
        [3] Bicubic Spline Interpolation
 
        [4] B-Spline Interpolation
 
</pre>
 
Before we can resample the <tt>dist_meuse</tt> map, we need to export it to the SAGA format:
 
 
 
<pre>
 
> writeGDAL(meuse.grid,"dist_meuse.sdat", mvFlag=-99999)
 
</pre>
 
 
 
To resample the map from 40 to 80 cell size, we need to set only one parameter -- the grid cell size:
 
<pre>&gt; rsaga.geoprocessor(lib="grid_tools", module=0, param=list(INPUT="dist_meuse.sgrd",
 
+        GRID="dist_meuse_80m.sgrd", METHOD=0, DIMENSIONS_CELLSIZE=80, SCALE_UP_METHOD=1))
 
</pre>
 
This gives the following image:
 
 
 
[[Image:Fig rescaling SAGA.jpg|450px|Rescaling maps in SAGA GIS.]]
 
 
 
Note also that the most automated approach to resample maps from a local coordinate system to LatLonWGS84 is to use the proj4 module in SAGA (note that proj4 string needs to be send with double quotes!). In this case SAGA will estimate the new grid for us, so that we only need to define the correct proj4 string of the input map:
 
 
<pre>
 
> rsaga.geoprocessor(lib="pj_proj4", 2, param=list(SOURCE_PROJ=paste('"', NL_RD, '"', sep=""), TARGET_PROJ="\"+proj=longlat +datum=WGS84\"",
 
+    SOURCE="meuse_dist.sgrd", TARGET="meuse_dist_ll.sgrd", TARGET_TYPE=0, INTERPOLATION=2))
 
# export to PNG:
 
rsaga.geoprocessor(lib="io_grid_image", 0, param=list(GRID="meuse_dist_ll.sgrd", FILE="meuse.png"))
 
</pre>
 
 
 
<br>
 
[[Image:Fig proj4 grid resampling in SAGA.jpg|450px|Resampling of maps in SAGA GIS using the Proj.4 module (Command Line Arguments, Grid).]]
 
 
 
Another excellent package to reproject rasters to LatLonWGS84 system is the gdalwarp utility available via the [[Software#FWTools|FWTools]].
 
 
 
== Time-series ==
 
 
 
Time-series is a range of maps with the same grid topology/definition, but where each map refers to a different time period. Note that also a list of remote sensing image bands can be exported as a time-series. Maptools package is not able to generate GE KMLs using space-time data, but we can instead write directly a KML file from R. For example, imagine that we have a list of maps produced by geostatistical simulations ([http://asdar-book.org Bivand et al., 2008]; p.230):
 
 
 
<pre>
 
> s1.fit <- fit.variogram(variogram(I(soil == 1) ~ 1, meuse),vgm(1, "Sph", 800, 1))
 
# geostatistical simulations of soil type "1":
 
> s1.sim <- krige(I(soil == 1)~1, meuse, meuse.grid, s1.fit, nsim=10, indicators=TRUE, nmax=40)
 
> spplot(s1.sim)
 
</pre>
 
 
 
We would like to visualize these as a (time-)series of ground overlays. We proceed with the same steps as described previously --- first, we reproject the maps to the WGS84 coordinate system:
 
 
 
<pre>
 
> s1.sim.ll <- projectRaster(stack(s1.sim), crs="+proj=longlat +datum=WGS84")
 
> s1.sim.grid <- as(s1.sim.ll, "SpatialGridDataFrame")
 
# export each map as PNG:
 
> grids.kml <- GE_SpatialGrid(s1.sim.grid)
 
> for(i in 1:ncol(s1.sim.ll@data@values)) {
 
+    png(file=paste("S1_sim", i,".png",sep=""), width=grids.kml$width, height=grids.kml$height, bg="transparent")
 
+    par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
 
+    image(s1.sim.ll[[i]], col=bpy.colors(), xlim=grids.kml$xlim, ylim=grids.kml$ylim)
 
+    dev.off()
 
> }
 
</pre>
 
 
 
We can now generate a KML where all maps will be separated using an equal time-span:
 
 
 
<pre>
 
> filename = "soil_sims.kml"
 
> write('<?xml version="1.0" encoding="UTF-8"?>', filename)
 
> write('<kml xmlns="http://earth.google.com/kml/2.2">', filename, append=TRUE)
 
> write('<Folder>', filename, append=TRUE)
 
> write('  <name>Simulations for S1</name>', filename, append=TRUE)
 
> write('  <open>1</open>', filename, append=TRUE)
 
> for(i in 1:ncol(s1.sim.ll@data@values)) {
 
>  write('  <GroundOverlay>', filename, append=TRUE)
 
>  write(paste('  <name>sim',i,'</name>',sep=""), filename, append=TRUE)
 
>  write('    <TimeSpan>', filename, append=TRUE)
 
>  write(paste('  <begin>', i,'</begin>',sep=""), filename, append=TRUE)
 
>  write(paste('  <end>', i+1,'</end>',sep=""), filename, append=TRUE)
 
>  write('    </TimeSpan>', filename, append=TRUE)
 
>  write('    <color>99ffffff</color>', filename, append=TRUE)
 
>  write('    <Icon>', filename, append=TRUE)
 
>  write(paste('  <href>',getwd(),'/S1_sim', i,'.png</href>',sep=""), filename, append=TRUE)
 
>  write('    <viewBoundScale>0.75</viewBoundScale>', filename, append=TRUE)
 
>  write('    </Icon>', filename, append=TRUE)
 
>  write('    <altitude>50</altitude>', filename, append=TRUE)
 
>  write('    <altitudeMode>relativeToGround</altitudeMode>', filename, append=TRUE)
 
>  write('  <LatLonBox>', filename, append=TRUE)
 
>  write(paste('    <north>',grids.kml$ylim[[2]],'</north>',sep=""), filename, append=TRUE)
 
>  write(paste('    <south>',grids.kml$ylim[[1]],'</south>',sep=""), filename, append=TRUE)
 
>  write(paste('    <east>',grids.kml$xlim[[2]],'</east>',sep=""), filename, append=TRUE)
 
>  write(paste('    <west>',grids.kml$xlim[[1]],'</west>',sep=""), filename, append=TRUE)
 
>  write('  </LatLonBox>', filename, append=TRUE)
 
>  write(' </GroundOverlay>', filename, append=TRUE)
 
> }
 
> write('</Folder>', filename, append=TRUE)
 
> write('</kml>', filename, append=TRUE)
 
</pre>
 
 
 
Note that the dates shown in Google Earth time-bar are artificial:
 
 
 
[[Image:Fig_meuse_simulations_GE.jpg|350px|Visualization of multiple simulations in GE.]]
 
 
 
= Reading KML files back to R  =
 
 
 
In principle, everything we export from R to KML we can also read back to R. Unfortunately, this is not an automatic procedure ([http://www.gdal.org/ogr/drv_kml.html GDAL] can only be used to write to KML) and the code often requires editing, for simple reason that KML does not (yet) have standard spatial data formats. The most recent version of the package [http://cran.r-project.org/web/packages/maptools/ maptools] contains methods (e.g. <tt>getKMLcoordinates</tt>) to read KML files to R and coerce them to the sp formats. KML could carry a lot of non-spatial information, or spatial information that is not supported by most of GIS (e.g. viewing angle, transparency). In addition, you could have more maps within the same KML file, which all makes it difficult to import a KML file to a GIS.
 
 
 
Via maptools, you can try importing some KML files (especially polygon maps) via e.g.:
 
 
 
<pre>
 
> meuse.coords <- getKMLcoordinates("meuse_lead.kml", ignoreAltitude=TRUE)
 
> str(meuse.coords) 
 
</pre>
 
 
 
but this has problems getting both coordinates and values for point files.
 
 
 
To read the point map <tt>meuse_lead.kml</tt> we have exported previously using the writeOGR method, we can use the [http://cran.r-project.org/web/packages/XML/ '''XML'''] package:
 
 
 
<pre>
 
&gt; library(XML)
 
&gt; meuse_lead.kml &lt;- xmlTreeParse("meuse_lead.kml")
 
&gt; lengthp &lt;- length(meuse_lead.kml$doc[[1]][[1]][[1]])-1
 
&gt; lead_sp &lt;- data.frame(Longitude=rep(0,lengthp), Latitude=rep(0,lengthp), Var=rep(0,lengthp))
 
&gt; for(j in 1:lengthp) {
 
&gt;    LatLon &lt;- unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][2][[1]][[1]][[1]])$value
 
&gt;    Var &lt;- unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][1][[1]][[1]])$value
 
&gt;      lead_sp$Longitude[[j]] &lt;- as.numeric(matrix(unlist(strsplit(LatLon, split=",")), ncol=2)[1])
 
&gt;      lead_sp$Latitude[[j]] &lt;- as.numeric(matrix(unlist(strsplit(LatLon, split=",")), ncol=2)[2])
 
&gt;      lead_sp$Var[[j]] &lt;- as.numeric(matrix(unlist(strsplit(strsplit(Var, split="&lt;i&gt;")[[1]][2], split="&lt;/i&gt;")), ncol=2)[1])
 
&gt; }
 
&gt; coordinates(lead_sp) &lt;-~Longitude+Latitude
 
&gt; proj4string(lead_sp) &lt;- CRS("+proj=longlat +ellps=WGS84")
 
&gt; bubble(lead_sp, "Var")
 
</pre>
 
 
 
Note that it will take time until you actually locate where in the KML file are the coordinates of points and attribute values, after that it is relatively easy to automate creation of a <tt>SpatialPointsDataFrame</tt>. This code could be shorten by using the <tt>xmlGetAttr()</tt>, <tt>xmlChildren()</tt> and <tt>xmlValue()</tt> methods. You might also consider using this [http://arcscripts.esri.com/details.asp?dbid=14988 ArcView script] to read KML files (points, lines, polygons) and generate shape files directly.
 
 
 
<br>
 
 
 
----
 
 
 
<rating>
 
Rate this article:
 
1 Poor and in error
 
2 Far from ready
 
3 Needs improvements
 
4 Minor corrections
 
5 Very useful!!
 
</rating>
 

Latest revision as of 10:05, 8 September 2012


  1. REDIRECTS to plotKML tutorial