[R-sig-Geo] selection of data subsets in spatial classes
Edzer J. Pebesma
e.pebesma at geo.uu.nl
Fri Jun 9 11:02:50 CEST 2006
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
>
More information about the R-sig-Geo
mailing list