How to produce spatial grid from raster?

Question

I need a Spatial Grid as a master grid for diverse thematic maps. How do I produce a Spatial Grid from a raster discarding all NA pixels?

Answer 1

You can get all the non-NA coordinates of cells in a raster with:

r = raster(matrix(runif(20),5,4))
r[r>.5]=NA

coordinates(r)[!is.na(values(r)),]
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1

those are the cells that are not NA. You can then probably feed these to SpatialPixels

SpatialPixels(SpatialPoints(coordinates(r)[!is.na(values(r)),]))
Object of class SpatialPixels
Grid topology:
  cellcentre.offset cellsize cells.dim
x             0.125     0.25         4
y             0.100     0.20         4
SpatialPoints:
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1
Coordinate Reference System (CRS) arguments: NA 

Although personally anything on a grid I'd keep as a raster.

I'm not totally sure still what it is you want - SpatialGrid objects are define full rectangular grids, so one without the NA pixels doesn't make sense.

Answer 2

Your two requirements seem to be about different things:

1) Some kind of reliable raster grid template.

2) A sparse grid that doesn't explicitly store missing cells.

sp::GridTopology provides the first one, it's just a description of the grid based on the lower left cell coordinate (cellcentre.offset), the cell spacing (cellsize), and the dimensions of the grid (cells.dim).

The sp::SpatialPixelsDataFrame class allows you to store sparse grids, but on its own it stores much more than the "template" - it also stores every coordinate explicitly. This is because it does two jobs, one it allows you to preserve the original coordinates that come from the grid and are possibly slightly irregular, two it lets you store only those cells that have valid values. (Arguably* these two goals should have been separated, but that's another story).

I don't think the raster package has an explicit analogue to GridTopology, but you can get hold of the components to "roll your own":

library(raster)
r1 <- raster(nrows=108, ncols=21, xmn=0, xmx=10)

## "cellsize"
res(r1)
## [1] 0.4761905 1.6666667

## extreme cell corners (just a different convention to sp's cellcentre)
extent(r1)
class       : Extent 
xmin        : 0 
xmax        : 10 
ymin        : -90 
ymax        : 90 

## we can also use bbox to get the same thing
bbox(r1)
min max
s1   0  10
s2 -90  90

## grid dimensions (including number of attributes/layers as 3rd "dim")
dim(r1)
## [1] 108  21   1

Linking these all together, we can go from raster to sp:

GridTopology(bbox(r1)[,1] + res(r1)/2, res(r1), dim(r1)[2:1])

(Note how the dimensions have to be reversed). Another simpler way is to coerce to SpatialGrid and use sp's getGridTopology, though this is more expensive since the coordinates end up being generated along the way:

getGridTopology(as(r1, "SpatialGrid"))

Those three parts of the raster "topology" are not all necessary, since some are redundant but there isn't otherwise a way to capture them all as one object - except, that raster created above is "empty" and so it can do the job that GridTopology does for sp. I'm not sure on the details of how "empty" it is, but it certainly does not explicitly store the "data" slot and is smaller than it would be with values in it. The raster package does in general do its best to keep memory usage to a minimum, and so with it you might not need to worry about really being "sparse".

That might help explain a little more, I know I'm overlapping Spacedman's answer but it's still not clear exactly what you mean in the question yet.

• (Arguably, since you can store sparse cells by just storing their index rather than explicit coordinates, and originally "slightly irregular" cell coordinates could just be stored as attributes if you really wanted them. Neither sp or raster deal with irregular rasters at all, not even the simple rectilinear case - this is in keeping with most GIS tools, but is unfortunate since it's quite common in formats like NetCDF and is handled by R's good old graphics::image function (though not if the recent rasterImage option is used). )

Answer 3

To transform a RasterLayer to a Spatial* object (Grid or Pixels) you can use the coercion function "as"

library(raster)
r <- raster(matrix(runif(20),5,4))
r[r>.5] <- NA
g <- as(r, 'SpatialGridDataFrame')
p <- as(r, 'SpatialPixels')

plot(r)
points(p)