[R-sig-Geo] How to find distance between grid cell and point that lies in specified radius?

[Ben Tupper]

Hi,

I think you'll enjoy sf once you get going with it.

I am not familiar with the Cressman scheme, but you maybe able to do
step 1 with raster (either pointDistance() or distanceFromPoints() - I
can't quite recall).  In sf you could something like this..

library(raster)
library(sf)
library(units)
library(dplyr)

crs <- "+proj=poly +lat_0=0 +lon_0=-54 +x_0=5000000 +y_0=10000000
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
r <- raster(nrows=9, ncols=18,
            xmn=5153337, xmx=6570069,
            ymn=7462732,ymx=8060416,
            crs = crs)
r[] <- runif(ncell(r))

xy <- raster::xyFromCell(r, seq_len(ncell(r)))
xy <- lapply(seq_len(nrow(xy)), function(i) sf::st_point(xy[i,]))
xy <- sf::st_sfc(xy, crs = crs)

pts <- sampleRandom(r, 10, na.rm = TRUE, sp = TRUE) %>%
  sf::st_as_sfc()

d <- sf::st_distance(xy, pts)
threshold <- units::set_units(100000, m)
d_within <- d <= threshold


d will be a [ncell x npoints] numeric matrix of distance from cell
locations to points
d_within will be a matrix [ncell x npoints] of logicals flagging which
are within your threshold (I used 100km)

Hope that helps start you off.

Cheers,
Ben



On Wed, Feb 5, 2020 at 9:12 PM Thiago V. dos Santos via R-sig-Geo
<r-sig-geo using r-project.org> wrote:
>
> Hello all,
>
> I am trying to implement a function in R to perform the Cressman scheme (which corrects the values of a gridded field, say precipitation, based on the closest observations available).
>
> It looks like it hasn't been done yet, but please let me know if you are aware of any R package that implements the Cressman scheme.
>
> I have a raster and a spatial points layer (which I will call "stations"). A toy example is provided below:
>
> ```
> library(raster)
>
>
> # create random raster
> r <- raster(nrows=9, ncols=18, xmn=5153337, xmx=6570069, ymn=7462732,
>             ymx=8060416, crs = CRS("+proj=poly +lat_0=0 +lon_0=-54 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
> r[] <- runif(ncell(r))
>
>
> # create spatial points
> pts <- sampleRandom(r, 10, na.rm = TRUE, sp = TRUE)
>
>
> # display everything
> plot(r)
> plot(pts, add=T, col="red", pch=16)
> ```
>
> The starting point of my analysis would be:
>
> 1) for each grid cell in the raster, identify the station(s) that lie *within a radius of influence* (say 30 km);
>
> 2) for each station lying under the radius of influence, calculate a distance-weighted expression that will be used later to correct the raster values. Its formula is:
>
> W = (R2 - r2)/(R2 + r2)
>
> where R = influence radius and r = distance between the station and the grid cell.
>
> 3) Based on the weighted correction factor W calculated above, update the value of the grid cell using the expression W * (gridpoint - station). The actual expression is a bit more sophisticated, but I want to get the general method here.
>
> Based on some research, it looks like this could be done with sf. But since I don't have any experience with that package, would someone more familiar please point out some functions in sf or even other packages that could be helpful to solve my problem?
>
> Thanks in advance,
>  -- Thiago V. dos Santos
>
> Postdoctoral Research Fellow
> Department of Climate and Space Science and Engineering
> University of Michigan
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org