[R-sig-Geo] Spatial correlation between two 'sf' kriging objects

[Bede-Fazekas Ákos]

Dear Nicola,
you are right, the class of the result of st_distance is not numeric but 
units (if CRS != NA). You may convert it to numeric with as.numeric().
HTH,
Ákos

2020.08.15. 23:15 keltezéssel, Nicola Gambaro írta:
> 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 using gmail.com>
>> To: r-sig-geo using r-project.org
>> Subject: Re: [R-sig-Geo]  Spatial correlation between two 'sf' kriging
>> 	objects
>> Message-ID: <043f13cc-20e1-a768-3419-03301f8d8f7b using 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 using r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo