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=500 +y_0=1000
+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(10, 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
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=500 +y_0=1000 +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@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
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
