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@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@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