Yes, the grid coordinates reflect the c(0,5, 0,5) [1,] 617400.0 6190061.0 [2,] 617400.5 6190061.0 [3,] 617401.0 6190061.0 [4,] 617401.5 6190061.0 [5,] 617400.0 6190060.5 [6,] 617400.5 6190060.5
But for some reason, when using the smaller cellsize my spplot of the interpolated field all of a sudden only shows a subsection of the original plot. There are no error messages. Karl |---------+----------------------------> | | [EMAIL PROTECTED]| | | no | | | | | | 09/06/2006 18:11 | | | Please respond to| | | Roger.Bivand | | | | |---------+----------------------------> >------------------------------------------------------------------------------------------------------------------------------| | | | To: [EMAIL PROTECTED] | | cc: r-sig-geo@stat.math.ethz.ch | | Subject: Re: [R-sig-Geo] selection of data subsets in spatial classes | >------------------------------------------------------------------------------------------------------------------------------| On Fri, 9 Jun 2006 [EMAIL PROTECTED] wrote: > > after some fiddling with the GridTopology() method I eventually succeeded > in making it work and now can interpolate onto the grid. > > Tremendous, thanks. Good > > One problem I encountered when specifying a number smaller than 1 for the > cellsize argument eg. c(0.5, 0.5) instead of c(1,1), I get erroneous > results. > > grd <- GridTopology(c(617400, 6190060), c(1,1), c(400, 300)) > I think this is a matter of screen representation: coordinates(GridTopology(c(617400, 6190060), c(0.5,0.5), c(4, 3))) looks wrong, but print(coordinates(GridTopology(c(617400, 6190060), c(0.5,0.5), c(4, 3))), digits=10) shows that the coordinate values were been rounded for display (is that what you are seeing?) Roger > Cheers > Karl > > > > |---------+----------------------------> > | | [EMAIL PROTECTED]| > | | no | > | | | > | | 08/06/2006 18:53 | > | | Please respond to| > | | Roger.Bivand | > | | | > |---------+----------------------------> > > ------------------------------------------------------------------------------------------------------------------------------| > | | > | To: [EMAIL PROTECTED] | > | cc: r-sig-geo@stat.math.ethz.ch | > | Subject: Re: [R-sig-Geo] selection of data subsets in spatial classes | > > ------------------------------------------------------------------------------------------------------------------------------| > > > > > On Thu, 8 Jun 2006 [EMAIL PROTECTED] wrote: > > > > > thanks Roger Bivand, this is really powerful and much less cumbersome > than > > I feared. > > > > I ran an inverse distance interpolation on some of the screened data. In > > order to do so I created a destination file from scratch with a 2 x 2 m > > grid. The destination file is of rectangular shape and is larger than the > > field for which a I have the aerial image data which is a polygon. As a > > result the kriging procedure extends beyond the field. > > > > Is there a way to either create a sampling grid that is the same shape as > > the field, ie from a shape file, or is there a way to subtract the two > > layers such that only the field shape is retained (something similar to > > mapcalc in GRASS)? > > My immediate guess would be that the overlay method will get you there: > > library(sp) > plot(1:100, 1:100, type="n") > poly <- locator() > poly <- do.call("cbind", poly) > poly <- rbind(poly, poly[1,]) > lines(poly) > field <- SpatialPolygons(list(Polygons(list(Polygon(poly)), ID="field"))) > grd <- GridTopology(c(0.5,0.5), c(1,1), c(100, 100)) > SG <- SpatialGrid(grd) > plot(SG, add=TRUE) > t1 <- overlay(SG, field) > SGDF <- SpatialGridDataFrame(grd, data=AttributeList(list(field=t1))) > SPDF <- as(SGDF, "SpatialPixelsDataFrame") > plot(SPDF, col="red", add=TRUE) > > and then use the SPDF points for prediction. If you've read in the field > polygon, you're already mostly set up. It will also work with multiple > fields, when t1 will know which field the grid point belongs to. > > Using GridTopology to set up prediction grids is quite neat too, it > provides an object to encapsulate what you are doing, but ensures that you > register the grid to the grid centres. > > Hope this helps, > > Roger > > > > > here is my R code > > > > library(rgdal) > > # load rgdal library and import ESRI raster file (aerial image) > > t1 <- readGDAL( "shi94a9_bnd.bil" , half.cell=c(0.5, 0.5), silent = > FALSE) > > > > image(t1, col=grey(1:9/10)) > > summary(t1) > > > > t1a <- as(t1, "SpatialPixelsDataFrame") > > summary(t1a) > > > > t1b <- t1a[t1a$band5 > 0,] # this eliminates values outside the field > > boundaries which are 0; inter-row pixels had been eliminated previously > > > > library(gstat) > > t1b.idw <- krige(band5~1, t1b, grid.pts) > > spplot(t1b.idw["var1.pred"],main="idw PCD map", > col.regions=bpy.colors(64)) > > > > # the grid.pts file was created as follows > > # > > utm.n.min<- 6190068 > > utm.n.max<- 6190326 > > utm.e.min<- 617405 > > utm.e.max<- 617812 > > > > spacing<-2 # average linear spacing between points > > # > > xseq<-seq(utm.e.min, utm.e.max, by=spacing) > > yseq<-seq(utm.n.min, utm.n.max, by=spacing) > > > > sample.pts<-expand.grid(xseq, yseq) > > names(sample.pts)<-c("x", "y") > > > > ## include the variable in a new object yld.pts > > grid.pts <- cbind(sample.pts) > > coordinates(grid.pts)=~x+y > > gridded(grid.pts)=TRUE > > > > Cheers > > > > Karl > > > > > > > > |---------+----------------------------> > > | | [EMAIL PROTECTED]| > > | | no | > > | | | > > | | 07/06/2006 19:17 | > > | | Please respond to| > > | | Roger.Bivand | > > | | | > > |---------+----------------------------> > > > > ------------------------------------------------------------------------------------------------------------------------------| > > > | > | > > | To: [EMAIL PROTECTED] > | > > | cc: r-sig-geo@stat.math.ethz.ch > | > > | Subject: Re: [R-sig-Geo] selection of data subsets in spatial > classes | > > > > ------------------------------------------------------------------------------------------------------------------------------| > > > > > > > > > > > On Wed, 7 Jun 2006 [EMAIL PROTECTED] wrote: > > > > > Hello list > > > > > > I have been trying to figure out if it is possible to select a subset > of > > > data in a SpatialGridDataFrame without prior conversion to a data.frame > > > using the method as.data.frame(). > > > > > > I have imported an ESRI " *. bil " file containing different spectral > > bands > > > using rgdal. The file imported as a SpatialGridDataFrame. The data > > > originates from an aerial photograph of a row crop and I would like to > > > select pixels from within a row as opposed to those from the inter-row > > > space for further interpolation. Inter-row pixels have a different > > > signature and therefore may be potentially screened out. However, I > have > > > been unable to access the SpatialGridDataFrame directly with, for > > example, > > > a subset(x, band5 < lowerLimit, select c(a, b) ) call for selecting > > values > > > according to given criteria. > > > > This is what the SpatialPixelsDataFrame class is for (unless you want to > > change resolution too): > > > > library(rgdal) > > SP27GTIF <- readGDAL(system.file("pictures/SP27GTIF.TIF", package = > > "rgdal")[1]) > > summary(SP27GTIF) > > image(SP27GTIF, col=grey(1:9/10)) > > SP27GTIF_a <- as(SP27GTIF, "SpatialPixelsDataFrame") > > summary(SP27GTIF_a) > > SP27GTIF_b <- SP27GTIF_a[SP27GTIF_a$band1 < 100,] > > summary(SP27GTIF_b) > > image(SP27GTIF_b) > > > > because SpatialPixelsDataFrame objects have their coordinates recorded > > explicitly, and "know where they are" on a full grid. > > > > fullgrid(SP27GTIF_b) <- TRUE > > summary(SP27GTIF_b) > > > > will drop them again, inserting NAs where band1 >= 100. > > > > This class is borrowed from the Terralib "CELL" data structure, part > > raster, part DB record, only recording cells with data. > > > > Roger > > > > > > > > Am I using the wrong strategy, should I just demote the spatial class > > using > > > method as.data.frame(), do the manipulation and then promote the > screened > > > values back to a spatial class? > > > > > > Regards > > > > > > Karl > > > _________________________________ > > > Karl J Sommer, > > > Department of Primary Industries, > > > PO Box 905 > > > Mildura, VIC, Australia 3502 > > > > > > Tel: +61 (0)3 5051 4390 > > > Fax +61 (0)3 5051 4534 > > > > > > Email: [EMAIL PROTECTED] > > > > > > _______________________________________________ > > > 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: [EMAIL PROTECTED] > > > > > > > > > > > > > > -- > 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: [EMAIL PROTECTED] > > > > > > -- 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: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo