2

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.

3

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.

  • Thanks, Jeffrey. I get an error when I try to coerce: "Error in is(object, Class) : trying to get slot "className" from an object of a basic class ("function") with no slots" and when I try call the @data: "Error in str(ks_sg@data) : no slot of name "data" for this object of class "SpatialGrid" ". – shea Mar 4 '16 at 16:21
  • 1
    Sorry, I did not notice that it is a SpatialGrid object. It does not have any attributes and will be of limited use until it is a SpatialGridDataFrame. I will modify my answer accordingly. I should note that I ran the example on my laptop with 8GB RAM and ran out of memory. I would point out again that this does not seem like a tractable approach and you should explore staying in a raster class (not sp SpatialPointsDataFrame) environment. – Jeffrey Evans Mar 4 '16 at 16:32
  • Jeffrey, you are right. I am helping someone else who want to use this approach and trying what I can. Thanks again. – shea Mar 4 '16 at 16:35
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"
)
  • Thanks for the suggestion of spsample. I did not write in my question that I had tried rasterToPoints and I ran in to memory issues. I'll think over sampling the data instead. – shea Mar 4 '16 at 16:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.