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

Roger Bivand Roger.Bivand at nhh.no
Fri Jun 9 13:08:28 CEST 2006


Karl,

When you reduce the cell sizes, are you increasing the cell dims 
proportionally, otherwise you will get what you are saying?

Roger


On Fri, 9 Jun 2006, Edzer J. Pebesma wrote:

> Karl, I cannot reproduce your problem:
> 
> library(sp)
> grd <- GridTopology(c(617400, 6190060), c(.5,.5), c(400, 300))
> grdx = SpatialGridDataFrame(grd, data = data.frame(z= 1:120000))
> spplot(grdx, col.regions=bpy.colors())
> 
> shows exactly what I expected.
> --
> Edzer
> 
> karl.sommer at dpi.vic.gov.au wrote:
> > 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
> >>>> _________________________________
> >>>> 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
> >
> > _______________________________________________
> > 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




More information about the R-sig-Geo mailing list