Roger of course, the problem was, as you suggested, that I did not adjust the cell dimensions.
Thanks again Karl |---------+----------------------------> | | [EMAIL PROTECTED]| | | no | | | | | | 09/06/2006 21:08 | | | Please respond to| | | Roger.Bivand | | | | |---------+----------------------------> >------------------------------------------------------------------------------------------------------------------------------| | | | To: [EMAIL PROTECTED] | | cc: [EMAIL PROTECTED], r-sig-geo@stat.math.ethz.ch | | Subject: Re: [R-sig-Geo] selection of data subsets in spatial classes | >------------------------------------------------------------------------------------------------------------------------------| Karl, When you reduce the cell sizes, are you increasing the cell dims proportionally, otherwise you will get what you are saying? Roger On Fri, 9 Jun 2006, Edzer J. Pebesma wrote: > Karl, I cannot reproduce your problem: > > library(sp) > grd <- GridTopology(c(617400, 6190060), c(.5,.5), c(400, 300)) > grdx = SpatialGridDataFrame(grd, data = data.frame(z= 1:120000)) > spplot(grdx, col.regions=bpy.colors()) > > shows exactly what I expected. > -- > Edzer > > [EMAIL PROTECTED] wrote: > > 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 > > > -- 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