selection of data subsets in spatial classes
Roger of course, the problem was, as you suggested, that I did not adjust the cell dimensions. Thanks again Karl |---------+----------------------------> | | Roger.Bivand at nhh.| | | no | | | | | | 09/06/2006 21:08 | | | Please respond to| | | Roger.Bivand | | | | |---------+----------------------------> >------------------------------------------------------------------------------------------------------------------------------| | | | To: e.pebesma at geo.uu.nl | | cc: karl.sommer at dpi.vic.gov.au, r-sig-geo at stat.math.ethz.ch | | Subject: Re: [R-sig-Geo] selection of data subsets in spatial classes | >------------------------------------------------------------------------------------------------------------------------------| 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