4

I have a SpatialPointsDataFrame with thousands of datapoints in a relatively small geographic area. The points are grouped around regions of physical measurements and are hence irregularly spaced:

plot(old)

old data

I'd like to resample these data onto a coarser grid that is sparsely populated so as not to overextrapolate observations. Each point in the grid grid should take on the value of the closest point in old --- unless there are no data within (say) a 500m radius, in which case the grid point in new should take on the value of NA.

ggg <-raster(crs = old@proj4string,
               ext = extent(old@bbox),
               resolution = 1000)
gg <- as(pppp, 'SpatialGrid')
plot(gg)

grid for new data

I have played around with creating an empty SpatialGrid object as shown above and using over(), but I get a vector of NAs. I assume that's because the points in the grid do not precisely align with my data.

I am currently trying to see if there's a way to do it with aggregate but no dice so far. Note that I don't want to average all values in each grid cell---instead, I want to use the nearest single observation to each grid cell center / point.

Can anyone point me in the right direction?

2

You could treat this as a (degenerate) interpolation problem, using a neighbourhood of one, try something along the lines of

library(gstat)
i = idw(var~1, old, ggg, nmax = 1, maxdist = 500)
library(sp)
spplot(i[1])

where var refers to the name of the variable you sampled.

  • Thank you. Since posting, I solved my problem. I will post it here for reference. I ended up discovering nncross() in spatstat. However, thanks for the response. As a general question: what's degenerate interpolation? Wikipedia didn't help me there? (And what, if any, is the label given to interpolation methods that are not degenerate)? Thanks. – Kevin Jul 15 '16 at 4:00
  • (just edited) - I meant interpolation with a single nearest observation, which does not really interpolate (basically spatial k-NN, with k = 1, but generating missing value if value is further than 500) – Edzer Pebesma Jul 16 '16 at 9:34
1

I found I was able to solve the problem using the nncross() function in the package spatstat. I don't know if this was the simplest solution, since it required first converting my SpatialPoints to class ppp. In any case here it is:

library('raster')
library('spatstat')

# Make the grid. Use raster() and then convert.
# Not the most direct but it works.
gridRes = 1000 # meters
ggg <-raster(crs = old@proj4string,
             ext = extent(old@bbox),
             resolution = gridRes)
gg <- as(ggg, 'SpatialGrid')
g <- as(gg, 'SpatialPoints')

grd <- ppp(x = g@coords[,1],
           y = g@coords[,2],
           window = owin(xrange=g@bbox[1,], yrange = g@bbox[2,]))
old_ppp <- ppp(x = old@coords[,1],
           y = old@coords[,2],
           window = owin(xrange=old@bbox[1,],
           yrange = old@bbox[2,]))

distz <- nncross(grd, old_ppp)
# distz is a dataframe with a row for each item in grd,
# and each row contains (1) distance and (2) reference to the nearest
# point in old_ppp.
  • 2
    The grumpy old man here wants to edit out all the @ uses here with proj4string(), extent(), coordinates(), xmin(), xmax(), ymin(), ymax(). – mdsumner Jul 15 '16 at 8:41
  • As a new R user this was far from clear the first time I read the rgdal / sp doc. Took ages just to stumble on how to extract information I needed from the datastructures :) Thanks for the pointer; will try to remember this! – Kevin Jul 18 '16 at 1:06

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.