How to find distance between grid cell and point that lies in specified radius?

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

How to find distance between grid cell and point that lies in specified radius?

R-sig-geo mailing list
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
Reply | Threaded
Open this post in threaded view
|

Re: 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
<[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