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.

Thanks,
Giovanni

2008/1/27, Roger Bivand <Roger.Bivand@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@bbox["r1","min"],0),
> > round(poly.shp@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@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@bbox["coords.x1","max"]-kbw2000.grd@bbox["coords.x1
> ","min"])/2000
> > height =
> > (kbw2000.grd@bbox["coords.x2","max"]-kbw2000.grd@bbox["coords.x2
> ","min"])/2000
> > geogrd.cell = (kbw2000.ll@bbox["s1","max"]-kbw2000.ll@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@grid@cells.dim["x1"]*geogrd@grid@cells.dim["x2"]
> > geogrd = SpatialGridDataFrame(geogrd@grid, data=data.frame(rep(1,
> > nogrids)), proj4string=kbw2000.ll@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@coords[,1], dat.SP@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@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@nhh.no
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]


