How to create cohesive Spatial Pixels from Spatial Points Dataset

Question

I have a Spatial Point DF spo (covering an irregular shaped area of interest). The data are not on a regular grid due to crs transformation.

My goal is a raster with predefined resolution and extent of the area of interest ( more spatial point data are to be mapped on this master raster) .

Problems start when I

rasterize(spo, raster(ncol, nrow, extent, crs), spo$param)

I need to adjust nrowand ncol in a way so that I wont get moire patterns of NAs within my area of interest. I can't use a predefined (higher) resolution, since rasterize has no interpolation capabilities.

As a solution to this, I thought I might need some kind of Spatial Pixel DF spi, that covers my whole area of interest (just like meuse.grid in library(raster); data(meuse.grid)), and serves as a master grid. Then, I can use it to interpolate my data, e.g.

idw(param~1,spo,spi)

and by this, get full cover of my area of interest at my chosen resolution. But how can a SpatialPixelsDataFrame be produced from the point data?

So in my view, the question boils down to: How to produce meuse.grid from meuse dataset?

Maybe I'm taking the wrong approach here, so please let me know if more easily can achieved what I'm after, using a different way.

Answer 1

If you have a polygon that defines the boundary of your region of interest, (which you should), then it is straight forward. One approach is to use the polygrid function from geoR, which itself is just a wrapper for SpatialPoints, expand.grid and overlay

Lets assume that you have a polygon that defines your region of interest called called ROI

In this case I will create one from meuse.grid

 data(meuse.grid)
 coordinates(meuse.grid) = ~x+y
 x <- chull(meuse.grid@coords)
 borders <- meuse.grid@coords[c(x,x[1]),]

 ROI <- SpatialPolygons(list(Polygons(list(Polygon(borders)), ID = 'border')))

In reality, to use polygrid you only need the coordinates of the polygon that define your region of interest.

To create 10-m grid covering the area of this ROI you can create a call to polygrid

# get the bounding box for ROI an convert to a list
bboxROI <- apply(bbox(ROI), 1, as.list)
# create a sequence from min(x) to max(x) in each dimension
seqs <- lapply(bboxROI, function(x) seq(x$min, x$max, by= 10))

# rename to xgrid and ygrid
names(seqs) <- c('xgrid','ygrid')

thegrid <- do.call(polygrid,c(seqs, borders = list(ROI@polygons[[1]]@Polygons[[1]]@coords)))