Spatial correlation between two 'sf' kriging objects

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Spatial correlation between two 'sf' kriging objects

Nicola Gambaro
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Spatial correlation between two 'sf' kriging objects

Bede-Fazekas Ákos
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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
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: Spatial correlation between two 'sf' kriging objects

Nicola Gambaro
In reply to this post by Nicola Gambaro
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?= <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-sig-Geo]  Spatial correlation between two 'sf' kriging
> objects
> Message-ID: <[hidden email]>
> 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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


        [[alternative HTML version deleted]]

_______________________________________________
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: 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?= <[hidden email]>
>> To: [hidden email]
>> Subject: Re: [R-sig-Geo]  Spatial correlation between two 'sf' kriging
>> objects
>> Message-ID: <[hidden email]>
>> 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
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo