Spatial correlation between two 'sf' kriging objects

Dear ?kos,

Thank you very much for your help.

I ran your code but unfortunately it returns the error:

Error in Ops.units(distance_matrix[as_units(point_number), , drop = TRUE],  :
   		both operands of the expression should be "units" objects

Any ideas how to fix this?

Cheers!

Nicola Gambaro
BSc Environmental Geoscience, First Class
Durham University

Message: 1
Date: Fri, 14 Aug 2020 14:01:43 +0200
From: =?UTF-8?Q?Bede-Fazekas_=c3=81kos?= <bfalevlist at gmail.com>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo]  Spatial correlation between two 'sf' kriging
	objects
Message-ID: <043f13cc-20e1-a768-3419-03301f8d8f7b at gmail.com>
Content-Type: text/plain; charset="utf-8"; Format="flowed"

Dear Nicola,

Instead of raster::focal(), you can apply the cor() function in
combination with st_distance() (or st_buffer() and st_within()). E.g.
let's say that 'grid' is a POINT type sf object containing columns
'temperature' and 'yield':
distance_threshold <- 40
distance_matrix <- st_distance(grid)
grid$correlation <- vapply(X = 1:nrow(grid), FUN.VALUE = numeric(1), FUN
= function (point_number) {cor(x = st_set_geometry(grid,
NULL)[distance_matrix[point_number, distance_matrix[point_number, , drop
= TRUE] < distance_threshold , drop = TRUE], "temperature"], y =
st_set_geometry(grid, NULL)[distance_matrix[point_number,
distance_matrix[point_number, , drop = TRUE] < distance_threshold , drop
= TRUE], "yield"])})

Or something like this... I have not tested this code, and am sure that
it is not the most efficient solution.

For large, square grids, raster might be faster than sf. You can convert
your grid to RasterLayer with function rasterFromXYZ() combined with
st_coordinates(), st_set_geometry(, NULL) and cbind.data.frame(). There
might be more straightforward solutions for the conversion...

HTH,
?kos Bede-Fazekas
Hungarian Academy of Sciences


2020.08.13. 20:11 keltez?ssel, Nicola Gambaro ?rta:

I have created two ?sf? kriging objects (point vectors), one for temperature and another for agricultural yields. To make the grid and carry out the point interpolation, I have remained within the ?sf? package.

I would now like to create a spatial local correlation ?raster? between these two variables, as shown on this webpage https://statnmap.com/2018-01-27-spatial-correlation-between-rasters/ <https://statnmap.com/2018-01-27-spatial-correlation-between-rasters/>. However, in that example, they use the ?raster? package and the ?focal? function. I was wondering if there was a way of doing this within ?sf?, i.e. without having to change classes? If not, what is the best way to convert those objects into raster classes?

Here is an excerpt of my kriging code for reference:
library(sf)
sf_data <- st_as_sf(x = data, coords = c("longitude", "latitude"), crs = 4326)
library(gstat)
vgm_utci <- variogram(UTCI~1, sf_data)
utci_fit <- fit.variogram(vgm_utci, vgm("Gau"), fit.kappa = TRUE)
plot(vgm_utci, utci_fit)
istria <- read_sf(?./Istria_Boundary.shp")
istria <- istria$geometry
istria.grid <- istria %>%
   st_make_grid(cellsize = 0.05, what = "centers") %>%
   st_intersection(istria)
library(ggplot2)
ggplot() + geom_sf(data = istria) + geom_sf(data = istria.grid)
library(stars)
utci_krig <- krige(formula = sf_data$UTCI ~ 1, locations = sf_data,
                    newdata = istria.grid, model = utci_fit)


Thank you very much in advance,

Nicola
	[[alternative HTML version deleted]]

	[[alternative HTML version deleted]]