[R-sig-Geo] SpatialPoints to SpatialGridDataFrame problems
Roger Bivand
Roger.Bivand at nhh.no
Sat Apr 2 00:44:48 CEST 2011
On Fri, 1 Apr 2011, Hailey Eckstrand wrote:
> Hello All, I am trying to convert some points to a SpatialGridDataFrame
> and am running into difficulties. When I convert the
> SpatialPointsDataFrame into a gridded object, it gives many warning
> messages and the output is incorrect. It says that the cellsize is
> 500x500 which is incorrect, it is actually 50,000 by 50,000)
Looking more closely shows that something is odd in your data (pts.txt is
your data below):
library(sp)
pts <- read.table("pts.txt", header=TRUE)
coordinates(pts) <- c("x", "y")
cangrid.gtop <- GridTopology(c(-633000, -4422000), c(50000, 50000), c(6,
10))
SG <- SpatialGrid(cangrid.gtop)
plot(pts)
plot(SG, add=TRUE, pch=4)
# go full-screen to see that they don't match
sapply(apply(coordinates(pts), 2, unique), sort)
sapply(sapply(apply(coordinates(pts), 2, unique), sort), function(x)
unique(diff(x)))
shows that they are not on a 50k by 50k grid, but some are offset by 500m.
So gridded() was trying hard to guess what you wanted. If they are meant
to be on the grid you specify, continue constructing a SpatialGrid, and
use an over() method to find out how to associate the grid cells with the
data from the SpatialPointsDataFrame:
plot(as(pts, "Spatial"), axes=TRUE)
text(coordinates(pts), label=pts$i)
oo <- over(pts, SG)
unique(table(oo))
# only one point per SG cell
SGi <- rep(as.integer(NA), 60)
SGi[oo] <- pts$i
SGDF <- SpatialGridDataFrame(SG, data=data.frame(i=SGi))
plot(as(SGDF, "Spatial"), axes=TRUE)
text(coordinates(SGDF), label=SGDF$i)
Here i is just the row number, but you would need to create new variables
for each input variable. Once you've sorted out your data, you can assign
the proj4string values.
Hope this helps,
Roger
>
> my.pts <- SpatialPointsDataFrame(coords=spatial.coords, data=my.df,
> proj4string=CRS(p4s))
> gridded(my.pts) <- T
>
> Warning messages:
> 1: In points2grid(points, tolerance, round) :
> grid has empty column/rows in dimension 1
> 2: In points2grid(points, tolerance, round) :
> grid topology may be corrupt in dimension 1
> 3: In points2grid(points, tolerance, round) :
> grid has empty column/rows in dimension 2
> 4: In points2grid(points, tolerance, round) :
> grid topology may be corrupt in dimension 2
>
> str(my.pts)
> Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
> ..@ data :'data.frame': 35 obs. of 2 variables:
> .. ..$ X1: num [1:35] 140.5 84.2 66.3 74.7 98.2 ...
> .. ..$ X2: num [1:35] 155.8 73.9 57.4 72.6 88.4 ...
> ..@ coords.nrs : num(0)
> ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
> .. .. ..@ cellcentre.offset: Named num [1:2] -633000 -4422000
> .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
> .. .. ..@ cellsize : Named num [1:2] 500 500
> .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
> .. .. ..@ cells.dim : Named int [1:2] 503 902
> .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
> ..@ grid.index : int [1:35] 453405 453505 453605 453706 402904 403005
> 403105 403205 403305 403406 ...
> ..@ coords : num [1:35, 1:2] -532500 -482500 -432500 -382000 -633000
> ...
> .. ..- attr(*, "dimnames")=List of 2
> .. .. ..$ : NULL
> .. .. ..$ : chr [1:2] "x" "y"
> ..@ bbox : num [1:2, 1:2] -633250 -4422250 -381750 -3971250
> .. ..- attr(*, "dimnames")=List of 2
> .. .. ..$ : chr [1:2] "x" "y"
> .. .. ..$ : chr [1:2] "min" "max"
> ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
> .. .. ..@ projargs: chr " +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-110
> +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
> +towgs84=0,0,0"
>
>
> I have also tried to build the GridTopology & glue it together with my data:
>
> gtop.coords.min <- c(min(spatial.coords[,1]), min(spatial.coords[,2]))
> gtop.coords.max <- c(max(my.coords[,1]), max(my.coords[,2]))
> gtop.cell.size <- c(50000, 50000)
> gtop.cells.dim.xc <- length(seq(gtop.coords.min[1], gtop.coords.max[1],
> by=gtop.cell.size[1]))
> gtop.cells.dim.yc <- length(seq(gtop.coords.min[2], gtop.coords.max[2],
> by=gtop.cell.size[2]))
> gtop.cells.dim <- c(gtop.cells.dim.xc, gtop.cells.dim.yc)
> cangrid.gtop <- GridTopology(gtop.coords.min, gtop.cell.size,
> gtop.cells.dim)
> str(cangrid.gtop)
> Formal class 'GridTopology' [package "sp"] with 3 slots
> ..@ cellcentre.offset: num [1:2] -633000 -4422000
> ..@ cellsize : num [1:2] 50000 50000
> ..@ cells.dim : int [1:2] 6 10
>
> SpatialGridDataFrame(grid=cangrid.gtop, data=my.df, proj4string=CRS(p4s))
> Error in validityMethod(object) :
> unequal number of objects in full grid and data slot
>
>
> I have included an image of what the coordinates look like. Here is a
> description of my variables:
> p4s
> "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-110 +k=1 +x_0=0 +y_0=0
> +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
>
> spatial.coords
> x y
> 1 -532500 -4422000
> 2 -482500 -4422000
> 3 -432500 -4422000
> 4 -382000 -4422000
> 5 -633000 -4372000
> 6 -582500 -4372000
> 7 -532500 -4372000
> 8 -482500 -4372000
> 9 -432500 -4372000
> 10 -382000 -4372000
> 11 -583000 -4322000
> 12 -532500 -4322000
> 13 -482500 -4322000
> 14 -432500 -4322000
> 15 -382000 -4322000
> 16 -583000 -4272000
> 17 -532500 -4272000
> 18 -482500 -4272000
> 19 -432500 -4272000
> 20 -382000 -4272000
> 21 -583000 -4222000
> 22 -532500 -4222000
> 23 -482500 -4222000
> 24 -432500 -4222000
> 25 -583000 -4172000
> 26 -532500 -4172000
> 27 -482500 -4172000
> 28 -583000 -4122000
> 29 -532500 -4122000
> 30 -482500 -4122000
> 31 -633000 -4071500
> 32 -583000 -4072000
> 33 -532500 -4072000
> 34 -583000 -4021500
> 35 -633000 -3971500
>
> str(my.df)
> 'data.frame': 35 obs. of 2 variables:
> $ X1: num 140.5 84.2 66.3 74.7 98.2 ...
> $ X2: num 155.8 73.9 57.4 72.6 88.4 ...
>
> Many thanks,
> Hailey
>
--
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