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

<

[hidden email]> 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

>

[hidden email]
>

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_______________________________________________

R-sig-Geo mailing list

[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo