R use data in SpatialGrid to create SpatialPoints

Question

I am trying to make a grid of points to overlay on a raster of vegetation of Kansas, U.S.A. My idea is to use the points as a sort of centroid for the raster pixels, which are 30m resolution. When I have the points, or "centroids", I want to use getValues() to capture some attributes of the raster and add the values from getValues() to the points. The steps I have taken so far, from Bivand et al.:

> bbox(ks08)
         min       max
s1  229687.8  889807.8
s2 4094357.4 4434227.4
> cs <- c(30, 30)
> cc <- bbox(ks08)[,1]+ (cs/2)
> cd <- ceiling(diff(t(bbox(ks08)))/cs)
> p4k <- CRS(proj4string(ks08))
> p4k
CRS arguments:
 +proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs 
> ks_grid <- GridTopology(cellcentre.offset = cc, cellsize = cs, cells.dim = cd)
> ks_sg <- SpatialGrid(ks_grid, proj4string = p4k)
> summary(ks_sg)
Object of class SpatialGrid
Coordinates:
         min       max
s1  229687.8  889807.8
s2 4094357.4 4434227.4
Is projected: TRUE 
proj4string :
[+proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
Grid attributes:
   cellcentre.offset cellsize cells.dim
s1          229702.8       30     22004
s2         4094372.4       30     11329
> str(ks_sg)
Formal class 'SpatialGrid' [package "sp"] with 3 slots
  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots    
  .. .. ..@ cellcentre.offset: Named num [1:2] 229703 4094372
  .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
  .. .. ..@ cellsize         : num [1:2] 30 30
  .. .. ..@ cells.dim        : int [1:2] 22004 11329
  ..@ bbox       : num [1:2, 1:2] 229688 4094357 889808 4434227
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "s1" "s2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

When I view ks_sg, part of the output include SpatialPoints and then the output chokes the memory (249,000,000+ points) for a few moments while it prints all the points' coordinates. I don't know what to do next to get those points into a vector format, like SpatialPoints, that would be the equivalent of a ArcGIS shapefile. Making the vector of points is my goal at the moment.

After the reply from Dave-Evans: I did try rasterToPoints before this, but ran into memory issues. Setting up the grid seemed like my best alternative.

After reply from Jeffrey Evans: My ultimate goal is to set up a very large data frame/data.table that has a column for each year of cover. The column will be vegetation classes for the year for every pixel, like a cbind(getValues(ks08),getValues(ks09),...). Each row is a pixel and each column is a year and the cell for a pixel and year would have the vegetation cover at that location and time.

Answer 1

Directly addressing your question, you can use the coordinates from your "SpatialGrid" object to create a "SpatialPointsDataFrame" object. Please note that it is not necessary perform the first coercion (using as), I just wanted to get the data into the same object class as yours.

library(sp)
data(meuse.grid)
coordinates(meuse.grid) = ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) = TRUE
class( meuse.grid <- as(meuse.grid, "SpatialGrid") )

class( meuse.grid <- SpatialPointsDataFrame( coordinates(meuse.grid),
                 data.frame(ID=1:length(meuse.grid))) )
str(meuse.grid@data)

The raster package has made this considerably easier. In fact, I just retooled an old function to deprecate a method, much like yours, in leu of a raster approach. Please consider this approach as it is quite a bit more efficient and will result in the same object class that you are after.

library(raster)
e <- extent(229687.8, 889807.8, 4094357.4, 4434227.4)
prj = "+proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
r <- raster(e, resolution = 30, crs = CRS(prj) )
r[] <- 1:ncell(r)           
r.pts <- rasterToPoints(r)

There really is no need to print the entire object and, in this case, I believe that you intended to look at the structure of the resulting data.frame and not the entire object (which will print everything). To look at the structure of just the data.frame you need to call the @data slot eg., str(ks_sg@data).

Since you are representing data as a systematic point array, I question the choice of coercing to points. I highly doubt that this is necessary and if you flushed out your goals a bit more, some relevant alternative could be provided. If you are wanting to look a spatial summary or aggregate aground each cell you could use a focal approach rather than coercing to points and extracting values.

Answer 2

I wonder do you really need to make a point from every cell? It seems to defeat the purpose of the raster format. Perhaps you could generate a regular grid of points using spsample and then use over to extract the data you are interested in from the raster. spsample will allow you to generate a grid of points that can be sparser than your raster.

But if it needs to be done you could try using rasterToPoints from the raster package (example from rasterToPoints help pages):

r <- raster(nrow=18, ncol=36)
r[] <- runif(ncell(r)) * 10
r[r>8] <- NA
p <- rasterToPoints(r, spatial=T)

With the flag spatial=T the function will return a SpatialPointsDataFrame which can be easily exported to a shapefile with

writeOGR(
    p
    dsn="outputDir",
    layer="grid_of_points",
    driver="ESRI Shapefile"
)