[R-sig-Geo] selection of data subsets in spatial classes

[karl.sommer at dpi.vic.gov.au]

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)?

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