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 |
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 |
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 |
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 |
Free forum by Nabble | Edit this page |