[R-sig-Geo] selection of data subsets in spatial classes
Roger Bivand
Roger.Bivand at nhh.no
Thu Jun 8 10:53:56 CEST 2006
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
More information about the R-sig-Geo
mailing list