[R-sig-Geo] selection of data subsets in spatial classes
Roger Bivand
Roger.Bivand at nhh.no
Fri Jun 9 10:11:27 CEST 2006
On Fri, 9 Jun 2006 karl.sommer at dpi.vic.gov.au 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
>
>
>
> |---------+---------------------------->
> | | Roger.Bivand at nhh.|
> | | no |
> | | |
> | | 08/06/2006 18:53 |
> | | Please respond to|
> | | Roger.Bivand |
> | | |
> |---------+---------------------------->
> >------------------------------------------------------------------------------------------------------------------------------|
> | |
> | To: karl.sommer at dpi.vic.gov.au |
> | cc: r-sig-geo at stat.math.ethz.ch |
> | Subject: Re: [R-sig-Geo] selection of data subsets in spatial classes |
> >------------------------------------------------------------------------------------------------------------------------------|
>
>
>
>
> On Thu, 8 Jun 2006 karl.sommer at dpi.vic.gov.au 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
> >
> >
> >
> > |---------+---------------------------->
> > | | Roger.Bivand at nhh.|
> > | | no |
> > | | |
> > | | 07/06/2006 19:17 |
> > | | Please respond to|
> > | | Roger.Bivand |
> > | | |
> > |---------+---------------------------->
> > >
> ------------------------------------------------------------------------------------------------------------------------------|
>
> > |
> |
> > | To: karl.sommer at dpi.vic.gov.au
> |
> > | cc: r-sig-geo at stat.math.ethz.ch
> |
> > | Subject: Re: [R-sig-Geo] selection of data subsets in spatial
> classes |
> > >
> ------------------------------------------------------------------------------------------------------------------------------|
>
> >
> >
> >
> >
> > On Wed, 7 Jun 2006 karl.sommer at dpi.vic.gov.au 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: karl.sommer at dpi.vic.gov.au
> > >
> > > _______________________________________________
> > > 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
> >
> >
> >
> >
> >
> >
>
> --
> 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
>
>
>
>
>
>
--
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
More information about the R-sig-Geo
mailing list