Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
4  
Throw us some scraps here. A little code to make a raster and what you expect of it? –  Spacedman Dec 11 '12 at 14:30
add comment

3 Answers

up vote 8 down vote accepted

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.

share|improve this answer
    
thanks mate, I'm also not sure what to want. I try to develop a workflow to produce identical raster maps (in terms of resolution, coorinates, crs, ...) from various spacial point data. rasterize is not an option due to spatial point scattering. Need to use iwd or similar. You are right, grid is rectengular--what I need as a master template probably is a SpatialPointsDataFrame. which can be gridded. S. –  sunpyg Dec 11 '12 at 16:09
    
You'll probably want a master raster (ha!) and to then create all subsequent grids concurrent to it. –  Radar Dec 11 '12 at 19:38
2  
Agreed. If things are gridded, raster is almost always the answer. Either that, or wading through about 20 pages of Applied Spatial Data Analysis in R to figure out how to do SpatialPixels. I gave that choice to a PhD student today, in fact. He chose... wisely! –  Spacedman Dec 11 '12 at 22:18
add comment

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). )
share|improve this answer
    
Thanks for the very instructive statements, which I need to digest. Hopefully enabling me to solve my problem or post a question that is more explicit! –  sunpyg Dec 12 '12 at 7:34
    
Also @Spacedman: I more precisely defined my question background here in the forum. –  sunpyg Dec 12 '12 at 13:02
add comment

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)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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