selection of data subsets in spatial classes

after some fiddling with the GridTopology() method I eventually succeeded
in making it work and now can interpolate onto the grid.

Tremendous, thanks.

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

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

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