Export maps to GE

From spatial-analyst.net
Revision as of 13:02, 21 September 2011 by Hengl (Talk | contribs)

Jump to: navigation, search
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. The complete script used to generate the plots shown down-below, is available here.

All displays in Google Earth are controlled by the 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 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 interactive KML sampler will give you some good idea about what is all possible in 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. 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 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 Google Earth User Guide.

Sample script to exchange maps to KML format.



Contents

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:

> require(c("rgdal","gstat","lattice","RASAGA","maptools","akima")) 

2. Reproject the original map from local coordinates (in this case Amersfoort / RD New) to the longlat system (where necessary):

> 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"))

3. Export the point map using the "KML" OGR driver:

> writeOGR(meuse.ll["lead"], "meuse_lead.kml", "lead", "KML") 

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:

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)

This will produce the following plot:

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.

# 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))))

Now we can plot the values of lead as 3D points (in this case we attach values of lead as heights):

> 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)

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:

> 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")

The actual plot in Google Earth can later be adjusted by manually editing the colors:

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 (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:

> data(meuse.grid)
> coordinates(meuse.grid) <- ~x+y
> gridded(meuse.grid) <- TRUE
> proj4string(meuse.grid) <- CRS(NL_RD)

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.

> 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

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 (Hengl and Evans, 2008):



\Delta x_{\mathit{m}} = F \cdot \cos(\varphi) \cdot \Delta x_{\mathit{m}}^0


where \delta x_{\mathit{m}} is the East/West grid cell size estimated for a given latitude \varphi, \Delta x_{\mathit{m}}^0 is the grid cell size in arcdegrees at equator, F=111,319 is an empirical constant used to convert from degrees to meters.

Reprojection of gridded map is also available via the 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:

> library(raster)
> dist.r <- raster(meuse.grid["dist"])
> dist.r.ll <- projectRaster(dist.r, crs="+proj=longlat +datum=WGS84")

We can look at the maps just to see that everything is fine:

> par(mfrow=c(1,2))
> image(dist.r)
> image(dist.r.ll)

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:

> 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>"

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:

> 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

This would produce an plot like this:

Google Earth visualization of an R plot.

Ground overlays

A ground overlay (with a legend) has the following structure:

<?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>

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:

> library(RSAGA)
> 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 :   Resampling
Usage: 0 -INPUT <str> [-GRID <str>] [-METHOD <num>] [-KEEP_TYPE] ...
  -INPUT:<str>                  Grid
        Grid (input)
  -GRID:<str>                   Resampled Grid
        Data Object (optional output)
  -METHOD:<num>                 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:<str>    Cell Size
        Floating point
        Minimum: 0.000000
...
  -GRID_GRID:<str>              Target Grid
        Grid (input)
  -SCALE_UP_METHOD:<num>        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:<num>      Interpolation Method
        Choice
        Available Choices:
        [0] Nearest Neighbor
        [1] Bilinear Interpolation
        [2] Inverse Distance Interpolation
        [3] Bicubic Spline Interpolation
        [4] B-Spline Interpolation 

Before we can resample the dist_meuse map, we need to export it to the SAGA format:

> writeGDAL(meuse.grid,"dist_meuse.sdat", mvFlag=-99999)

To resample the map from 40 to 80 cell size, we need to set only one parameter -- the grid cell size:

> 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))

This gives the following image:

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:

> 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"))


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 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 (Bivand et al., 2008; p.230):

> 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)

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:

> 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()
> }

We can now generate a KML where all maps will be separated using an equal time-span:

> 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)

Note that the dates shown in Google Earth time-bar are artificial:

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 (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 maptools contains methods (e.g. getKMLcoordinates) 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.:

> meuse.coords <- getKMLcoordinates("meuse_lead.kml", ignoreAltitude=TRUE)
> str(meuse.coords)  

but this has problems getting both coordinates and values for point files.

To read the point map meuse_lead.kml we have exported previously using the writeOGR method, we can use the XML package:

> library(XML)
> meuse_lead.kml <- xmlTreeParse("meuse_lead.kml")
> lengthp <- length(meuse_lead.kml$doc[[1]][[1]][[1]])-1
> lead_sp <- data.frame(Longitude=rep(0,lengthp), Latitude=rep(0,lengthp), Var=rep(0,lengthp))
> for(j in 1:lengthp) {
>     LatLon <- unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][2][[1]][[1]][[1]])$value
>     Var <- unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][1][[1]][[1]])$value
>      lead_sp$Longitude[[j]] <- as.numeric(matrix(unlist(strsplit(LatLon, split=",")), ncol=2)[1])
>      lead_sp$Latitude[[j]] <- as.numeric(matrix(unlist(strsplit(LatLon, split=",")), ncol=2)[2])
>      lead_sp$Var[[j]] <- as.numeric(matrix(unlist(strsplit(strsplit(Var, split="<i>")[[1]][2], split="</i>")), ncol=2)[1])
> }
> coordinates(lead_sp) <-~Longitude+Latitude
> proj4string(lead_sp) <- CRS("+proj=longlat +ellps=WGS84")
> bubble(lead_sp, "Var")

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 SpatialPointsDataFrame. This code could be shorten by using the xmlGetAttr(), xmlChildren() and xmlValue() methods. You might also consider using this ArcView script to read KML files (points, lines, polygons) and generate shape files directly.



<rating> Rate this article: 1 Poor and in error 2 Far from ready 3 Needs improvements 4 Minor corrections 5 Very useful!! </rating>

Personal tools