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?= <bfalevl...@gmail.com> > To: r-sig-geo@r-project.org > Subject: Re: [R-sig-Geo] Spatial correlation between two 'sf' kriging > objects > Message-ID: <043f13cc-20e1-a768-3419-03301f8d8...@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]] >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo