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
