selection of data subsets in spatial classes

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



|---------+---------------------------->
|         |           Roger.Bivand at nhh.|
|         |           no               |
|         |                            |
|         |           09/06/2006 18:11 |
|         |           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 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

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