Replicating GMT's blockmedian/blockmean in R
I have frequently used GMT’s blockmean and blockmedian functions with great success. These method essentially quasi-grid your points which is a particularly useful technique when trying to reduce irregularly spaced data such as for onward point sampling or data reduction. One case I have employed this is where the points were to be used as an input to kriging.
Sometimes I need to make use of these blockmean/blockmedian functions within R when system commands cannot be executed. Consequently, I’ve looked for a way to mimic the GMT functions within R - below is my current effort. This could be improved but provides a useful start if you’ve been looking for something similar… data with columns other than x and y are not handled in the below example.
The below code:
- reads in points
- creates a grid over which the points lie of a given dimension
- assigns unique ID to each cell (in this case, the cell number)
- extracts raster cell value to all points within a cell
- calculates new point positions for all points within a specific cell
- output new point list - blockmean or blockmedian filtered :)
library(sp)
library(raster)
# Create some points
xs=sample(seq(from=1, to=20, by=0.2), size=500, replace=TRUE)
ys=sample(seq(from=1, to=20, by=0.2), size=500, replace=TRUE)
plot(xs,ys)
pnts=SpatialPoints(cbind(xs,ys))
ext=extent(pnts)
# Create grid
nrows=ext[4]-ext[3]
ncols=ext[2]-ext[1]
r = raster(nrows = nrows, ncols = ncols, xmn = ext[1], xmx = ext[2], ymn = ext[3], ymx = ext[4])
# Assign cell number to each cell
vals=1:ncell(r)
r[]=vals
# Extract cell value to points
pnts$cell=extract(r, pnts)
# Calculate median point position
# Assumes dataframe with columns x,y,cell (<<< cell number)
# method - median or mean - otherwise fail
block_xy <-function(df, method='median'){
cells=unique(df$cell)
xys=data.frame()
for (cell_inst in cells){
grp=subset(df, cell==cell_inst)
if (method=='median'){
xs=median(grp$x)
ys=median(grp$y)
} else if (method=='mean'){
xs=mean(grp$x)
ys=mean(grp$y)
}
xys=rbind(xys,c(xs, ys, cell_inst))
}
colnames(xys)<-c('x','y','cell')
return(xys)
}
# Run the filtering
pnts=as.data.frame(pnts) # cell, xs, ys
colnames(pnts)=c('cell','x','y')
med_xys=block_xy(pnts, method='median')
mean_xys=block_xy(pnts, method='mean')
# Quick plot
plot(r)
plot(rasterToPolygons(r), add=TRUE, border='black', lwd=0.5)
points(pnts$x, pnts$y, pch=20, col='blue')
points(med_xys, pch=4, col='green')
points(mean_xys, pch=4, col='red')
Quasi-gridded points. Raster is coloured by cell number, original points are in blue, blockmedian points are in green and blockmean points are in red.