[R-sig-Geo] kernel density estimation to google earth

Mon Jan 28 10:51:27 CET 2008

On Mon, 28 Jan 2008, G. Allegri wrote:

> Thanks, a really interesting post.
> May I ask a stupid question? Where are the following function documented?
>
> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp
> )[[1]])[[1]])
>
> They are in the sp package but I can't find references but in some posts on
> this ML.

They are deprecated (will be defunct soon); this is why they generate 
warnings. This is because using slot() is much cleaner, easier to 
maintain, and is automatically documented. Say you want the coordinates of 
a "Polygon" object:

getSlots("Polygon")

tells you what they are called, and:

slot(x, "coords")

retrieves them. So an up-to-date version of the above is:

> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp
> )[[1]])[[1]])

slot(slot(slot(poly.shp, "polygons")[[1]], "Polygons")[[1]], "coords")

or equivalently:

poly.shp_1 <- slot(poly.shp, "polygons")[[1]]
poly.shp_1_1 <- slot(poly.shp_1, "Polygons")[[1]]
crds_1_1 <- slot(poly.shp_1_1, "coords")

where all you need are slot() and the names of the slots, rather than very 
long and rather arbitrary names for wrapper functions for slot().

In the code below, note that use of "@" should be avoided in user-level 
code, and that bbox() methods retrieve the bounding box of a Spatial* 
object cleanly.

Roger

>
> Thanks,
> Giovanni
>
> 2008/1/27, Roger Bivand <Roger.Bivand at nhh.no>:
>>
>> On Sun, 27 Jan 2008, Tomislav Hengl wrote:
>>
>>>
>>> # Modified by T. Hengl (http://spatial-analyst.net)
>>> # 27.01.2008
>>
>> Thanks, very clear - the kernel values at the grid of projected points are
>> not an a grid in geographical coordinates, and need to be interpolated
>> (warped) to a grid in geographical coordinates. The alternative is to make
>> the kernel function use Great Circle distances.
>>
>> Roger
>>
>>>
>>> library(rgdal)
>>> library(maptools)
>>> library(splancs)
>>>
>>> # Import the points and study area:
>>>
>>> data.shp <- readOGR("C:/", layer="events")
>>> str(data.shp)
>>>
>>> poly.shp <- readOGR("C:/", layer="hull")
>>> str(poly.shp)
>>>
>>> poly <-
>>> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp
>> )[[1]])[[1]])
>>> grd <-
>>> GridTopology(cellcentre.offset=c(round(poly.shp at bbox["r1","min"],0),
>>> round(poly.shp at bbox["r2","min"],0)), cellsize=c(2000,2000),
>>> cells.dim=c(30,25))
>>>
>>> # Run the 2D kernel smoother:
>>>
>>> kbw2000 <- spkernel2d(data.shp, poly, h0=2000, grd)
>>> hist(kbw2000)
>>>
>>> # Pack and plot a SpatialGridDataFrame:
>>>
>>> kbw2000.grd <- SpatialGridDataFrame(grd, data=data.frame(kbw2000))
>>> proj4string(kbw2000.grd) <- data.shp at proj4string
>>> spplot(kbw2000.grd, col.regions=terrain.colors(16),
>>> sp.layout=list("sp.points",pch="+",cex=1.2,col="black",data.shp))
>>>
>>> # Export to KML
>>> # First, reproject the grid to longlat:
>>>
>>> kbw2000.ll <- spTransform(kbw2000.grd, CRS("+proj=longlat
>> +datum=WGS84"))
>>> str(kbw2000.ll)
>>>
>>> # The cell size you need to determine yourself!!
>>>
>>> width =
>>> (kbw2000.grd at bbox["coords.x1","max"]-kbw2000.grd at bbox["coords.x1
>> ","min"])/2000
>>> height =
>>> (kbw2000.grd at bbox["coords.x2","max"]-kbw2000.grd at bbox["coords.x2
>> ","min"])/2000
>>> geogrd.cell = (kbw2000.ll at bbox["s1","max"]-kbw2000.ll at bbox
>> ["s1","min"])/width
>>>
>>> # Define a new grid:
>>>
>>> geogrd = spsample(kbw2000.ll, type="regular",
>>> cellsize=c(geogrd.cell,geogrd.cell))
>>> gridded(geogrd) = TRUE
>>>
>>> gridparameters(geogrd)
>>> # cellcentre.offset   cellsize cells.dim
>>> # x1          15.90165 0.02636685        30
>>> # x2          47.95541 0.02636685        16
>>>
>>> # This is an empty grid without any topology (only grid nodes are
>> defined)
>>> and coordinate
>>> # system definition. To create topology, we coerce a dummy variable
>> (1s),
>>> then
>>> # specify that the layer has a full topology:
>>>
>>> nogrids = geogrd at grid@cells.dim["x1"]*geogrd at grid@cells.dim["x2"]
>>> geogrd = SpatialGridDataFrame(geogrd at grid, data=data.frame(rep(1,
>>> nogrids)), proj4string=kbw2000.ll at proj4string)
>>>
>>> # and estimate the values of the reprojected map at new grid locations
>>> using the bilinear resampling:
>>> # this can be time-consuming for large grids!!!
>>>
>>> library(gstat)
>>> kbw2000.llgrd = krige(kbw2000~1, kbw2000.ll, geogrd, nmax=4)
>>>
>>> # Optional, convert the original shape to latlong coordinates:
>>>
>>> data.ll <- spTransform(data.shp, CRS("+proj=longlat +datum=WGS84"))
>>> spplot(kbw2000.llgrd["var1.pred"], col.regions=terrain.colors(16),
>>> scales=list(draw=TRUE),
>>> sp.layout=list("sp.points",pch="+",cex=1.2,col="black",data.ll))
>>>
>>> # The final grid map can be exported to KML format using the maptools
>>> package and kmlOverlay method:
>>>
>>> kbw2000.kml = GE_SpatialGrid(kbw2000.llgrd)
>>>
>>> tf <- tempfile()
>>> png(file=paste(tf, ".png", sep=""), width=kbw2000.kml$width,
>>> height=kbw2000.kml$height, bg="transparent")
>>>
>>> par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
>>> image(as.image.SpatialGridDataFrame(kbw2000.llgrd[1]), col=bpy.colors(),
>>> xlim=kbw2000.kml$xlim, ylim=kbw2000.kml$ylim)
>>> plot(data.ll, pch="+", cex=1.2, add=TRUE, bg="transparent")
>>>
>>> kmlOverlay(kbw2000.kml, paste(tf, ".kml", sep=""), paste(tf, ".png",
>> sep=""))
>>> dev.off()
>>>
>>> # see also:
>>> # Hengl, T., 2007. A Practical Guide to Geostatistical Mapping of
>>> Environmental Variables. EUR 22904 EN Scientific and Technical Research
>>> series, Office for Official Publications of the European Communities,
>>> Luxemburg, 143 pp.
>>>
>>>
>>>> Dear List,
>>>>
>>>> I've calculated a kernel density estimation (splancs function) and now
>> I
>>>> want to export the result as a kml-file to open it in the google earth
>>>> viewer, but I get stuck at converting the SpatialGridDataFrame to a
>>>> SpatialPixelDataFrame...
>>>>
>>>> here (http://www.zshare.net/download/689742375ef5fb/) you can download
>>>> some
>>>> testdata and my R-code so far:
>>>>
>>>> ###########
>>>> data.shp <- readOGR("C:/", layer="events")
>>>> prj <- data.shp@ proj4string@ projargs
>>>> dat <- data.shp
>>>> str(dat)
>>>> poly.shp <- readOGR("C:/", layer="hull")
>>>> str(poly.shp)
>>>>
>>>> dat.SP <- as(dat, "SpatialPoints")
>>>> str(dat.SP)
>>>> pp_poi <- as.points(dat.SP at coords[,1], dat.SP at coords[,2])
>>>> poly <-
>>>> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(
>> poly.shp)[[1]])[[1]])
>>>> polymap(poly)
>>>> points(pp_poi)
>>>>
>>>> grd <- GridTopology(cellcentre.offset=c(590511, 396191),
>> cellsize=c(2000,
>>>> 2000),
>>>>   cells.dim=c(30,25))
>>>> kbw2000 <- spkernel2d(pp_poi, poly, h0=2000, grd)
>>>> spplot(SpatialGridDataFrame(grd, data=data.frame(kbw2000)),
>>>> col.regions=terrain.colors(16))
>>>>
>>>> test <- SpatialGridDataFrame(grd, data=data.frame(kbw2000))
>>>> proj4string(test) <- CRS(prj)
>>>> str(test)
>>>> test1 <- spTransform(test, CRS("+proj=longlat +datum=WGS84"))
>>>>
>>>> ## here I got following error message:
>>>> # validityMethod(as(object, superClass)): Geographical CRS given to
>>>> non-conformant data
>>>> test2 <- spsample(test1, type="regular", cellsize=c(2000,2000))
>>>>
>>>> # export the SpatialPixelDataFrame (Code (not testet) from the
>> help-file)
>>>> tf <- tempfile()
>>>> SGxx <- GE_SpatialGrid(test2)
>>>> png(file=paste(tf, ".png", sep=""), width=SGxx$width,
>> height=SGxx$height,
>>>>   bg="transparent")
>>>> par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
>>>> plot(x, xlim=SGxx$xlim, ylim=SGxx$ylim, setParUsrBB=TRUE)
>>>> dev.off()
>>>> kmlOverlay(SGxx, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
>>>> ###########
>>>>
>>>> I appreciate every hint! Thanks.
>>>> Marco
>>>>
>>>> --
>>>> Marco Helbich
>>>> Institute for Urban and Regional Research
>>>> Austrian Academy of Sciences
>>>> Postgasse 7/4/2, A-1010 Vienna, Austria (EU)
>>>> e-mail: marco.helbich(at)oeaw.ac.at
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at stat.math.ethz.ch
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>> --
>> Roger Bivand
>> Economic Geography Section, Department of Economics, Norwegian School of
>> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
>> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: Roger.Bivand at nhh.no
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no