Rasterize overlapping sf objects in R and specify which values to use for pixels
I have an sf dataframe with multiple polygons in R (each row is a square in a slightly different location with a different value associated with it). There are also some gaps between the polygons. I would like to rasterize the sf dataframe using stars::st_rasterize (using a second sf dataframe as the extent of the new raster) but am running into 2 issues.
When there are multiple overlapping sf objects in a raster pixel, I would like to specify that I want the raster pixel to take the lowest possible value from the sf objects.
In the rasterizing process, I would like to fill in the gaps between the sf objects using interpolation based on the values of the nearest sf objects. In other words, I don't want the raster to have NA values within the extent of the second sf df.
Here are some dummy sf dataframes:
a) the polygon that will be the extent of the raster
### load libraries
library(tidyverse)
library(sf)
library(stars)
### set vertices
polygon_vertices <- matrix(c(-1489008, 1369848,
-1488191, 1370128,
-1488132, 1369347,
-1488985, 1369243,
-1489008, 1369848),
ncol = 2,
byrow = TRUE)
### make it a polygon
polygon_extent <- st_polygon(list(polygon_vertices))
### make it an sf object and set the crs
polygon_extent <- polygon_extent %>%
st_sfc(crs = 6350)
b) the sf dataframe with the overlapping sf objects
### set the vertices for the rows of the sf df
p1_vertices <- matrix(c(-1488355, 1369737,
-1488355, 1370072,
-1488237, 1370112,
-1488190, 1370112,
-1488162, 1369737,
-1488355, 1369737),
ncol = 2, byrow = TRUE)
p2_vertices <- matrix(c(-1488355, 1369737,
-1488355, 1370072,
-1488237, 1370112,
-1488190, 1370112,
-1488162, 1369737,
-1488355, 1369737),
ncol = 2, byrow = TRUE)
p3_vertices <- matrix(c(-1488253, 1369484,
-1488253, 1369859,
-1488171, 1369859,
-1488143, 1369484,
-1488253, 1369484),
ncol = 2, byrow = TRUE)
p4_vertices <- matrix(c(-1488292, 1369795,
-1488292, 1370093,
-1488191, 1370128,
-1488166, 1369795,
-1488292, 1369795),
ncol = 2, byrow = TRUE)
p5_vertices <- matrix(c(-1488260, 1369634,
-1488154, 1369634,
-1488132, 1369347,
-1488260, 1369332,
-1488260, 1369634),
ncol = 2, byrow = TRUE)
p6_vertices <- matrix(c(-1488255, 1369734,
-1488630, 1369734,
-1488630, 1369977,
-1488255, 1370106,
-1488255, 1369734),
ncol = 2, byrow = TRUE)
p7_vertices <- matrix(c(-1488240, 1369419,
-1488615, 1369419,
-1488615, 1369794,
-1488240, 1369794,
-1488240, 1369419),
ncol = 2, byrow = TRUE)
p8_vertices <- matrix(c(-1488212, 1370120,
-1488191, 1370128,
-1488190, 1370120,
-1488212, 1370120),
ncol = 2, byrow = TRUE)
p9_vertices <- matrix(c(-1488970, 1369519,
-1488995, 1369519,
-1489008, 1369848,
-1488970, 1369861,
-1488970, 1369519),
ncol = 2, byrow = TRUE)
p10_vertices <- matrix(c(-1488970, 1369519,
-1488995, 1369519,
-1489008, 1369848,
-1488970, 1369861,
-1488970, 1369519),
ncol = 2, byrow = TRUE)
p11_vertices <- matrix(c(-1488518, 1369823,
-1488518, 1370016,
-1488191, 1370128,
-1488168, 1369823,
-1488518, 1369823),
ncol = 2, byrow = TRUE)
p12_vertices <- matrix(c(-1488452, 1369986,
-1488452, 1370039,
-1488191, 1370128,
-1488180, 1369986,
-1488452, 1369986),
ncol = 2, byrow = TRUE)
p13_vertices <- matrix(c(-1488705, 1369618,
-1488999, 1369618,
-1489008, 1369848,
-1488705, 1369952,
-1488705, 1369618),
ncol = 2, byrow = TRUE)
p14_vertices <- matrix(c(-1488804, 1369607,
-1488999, 1369607,
-1489008, 1369848,
-1488804, 1369918,
-1488804, 1369607),
ncol = 2, byrow = TRUE)
p15_vertices <- matrix(c(-1488921, 1369704,
-1489002, 1369704,
-1489008, 1369848,
-1488921, 1369878,
-1488921, 1369704),
ncol = 2, byrow = TRUE)
p16_vertices <- matrix(c(-1488299, 1369901,
-1488674, 1369901,
-1488674, 1369962,
-1488299, 1370091,
-1488299, 1369901),
ncol = 2, byrow = TRUE)
p17_vertices <- matrix(c(-1488904, 1369779,
-1489005, 1369779,
-1489008, 1369848,
-1488904, 1369884,
-1488904, 1369779),
ncol = 2, byrow = TRUE)
p18_vertices <- matrix(c(-1488525, 1369926,
-1488780, 1369926,
-1488525, 1370013,
-1488525, 1369926),
ncol = 2, byrow = TRUE)
### combine all of them into an sf df
sf_boxes <- st_sf(pix_value = c(1:18), # Unique identifier for each polygon
geometry = st_sfc(st_polygon(list(p1_vertices)),
st_polygon(list(p2_vertices)),
st_polygon(list(p3_vertices)),
st_polygon(list(p4_vertices)),
st_polygon(list(p5_vertices)),
st_polygon(list(p6_vertices)),
st_polygon(list(p7_vertices)),
st_polygon(list(p8_vertices)),
st_polygon(list(p9_vertices)),
st_polygon(list(p10_vertices)),
st_polygon(list(p11_vertices)),
st_polygon(list(p12_vertices)),
st_polygon(list(p13_vertices)),
st_polygon(list(p14_vertices)),
st_polygon(list(p15_vertices)),
st_polygon(list(p16_vertices)),
st_polygon(list(p17_vertices)),
st_polygon(list(p18_vertices))),
crs = 6350)
So here's what that looks like:
ggplot() +
geom_sf(data = polygon_extent, color = "gray") +
geom_sf(data = sf_boxes, aes(color = as.character(pix_value),
fill = as.character(pix_value)),
alpha = 0.2) +
NULL
And then I can rasterize it, but I'm having trouble figuring out how to set the code so that 1) each raster pixel with overlapping sf objects from sf_boxes gets the lowest value assigned to it and 2) the missing values within the extent of polygon_extent are filled in.
### rasterize
dummy_raster <- st_rasterize(sf_boxes,
st_as_stars(st_bbox(polygon_extent),
nx = 399,
ny = 460))
plot(dummy_raster)
Any ideas? Thanks!
Solution
You can do that with terra::rasterize
Your polygons
wkt <- c('POLYGON ((-1488355 1369737, -1488355 1370072, -1488237 1370112, -1488190 1370112, -1488162 1369737, -1488355 1369737))', 'POLYGON ((-1488355 1369737, -1488355 1370072, -1488237 1370112, -1488190 1370112, -1488162 1369737, -1488355 1369737))', 'POLYGON ((-1488253 1369484, -1488253 1369859, -1488171 1369859, -1488143 1369484, -1488253 1369484))', 'POLYGON ((-1488292 1369795, -1488292 1370093, -1488191 1370128, -1488166 1369795, -1488292 1369795))', 'POLYGON ((-1488260 1369634, -1488154 1369634, -1488132 1369347, -1488260 1369332, -1488260 1369634))', 'POLYGON ((-1488255 1369734, -1488630 1369734, -1488630 1369977, -1488255 1370106, -1488255 1369734))', 'POLYGON ((-1488240 1369419, -1488615 1369419, -1488615 1369794, -1488240 1369794, -1488240 1369419))', 'POLYGON ((-1488212 1370120, -1488191 1370128, -1488190 1370120, -1488212 1370120))', 'POLYGON ((-1488970 1369519, -1488995 1369519, -1489008 1369848, -1488970 1369861, -1488970 1369519))', 'POLYGON ((-1488970 1369519, -1488995 1369519, -1489008 1369848, -1488970 1369861, -1488970 1369519))', 'POLYGON ((-1488518 1369823, -1488518 1370016, -1488191 1370128, -1488168 1369823, -1488518 1369823))', 'POLYGON ((-1488452 1369986, -1488452 1370039, -1488191 1370128, -1488180 1369986, -1488452 1369986))', 'POLYGON ((-1488705 1369618, -1488999 1369618, -1489008 1369848, -1488705 1369952, -1488705 1369618))', 'POLYGON ((-1488804 1369607, -1488999 1369607, -1489008 1369848, -1488804 1369918, -1488804 1369607))', 'POLYGON ((-1488921 1369704, -1489002 1369704, -1489008 1369848, -1488921 1369878, -1488921 1369704))', 'POLYGON ((-1488299 1369901, -1488674 1369901, -1488674 1369962, -1488299 1370091, -1488299 1369901))', 'POLYGON ((-1488904 1369779, -1489005 1369779, -1489008 1369848, -1488904 1369884, -1488904 1369779))', 'POLYGON ((-1488525 1369926, -1488780 1369926, -1488525 1370013, -1488525 1369926))')
library(terra)
v <- vect(wkt)
values(v) <- data.frame(ID=1:nrow(v))
Solution, use terra::rasterize with arguments fun and background
r <- rast(v, res=2)
x <- rasterize(v, r, "ID", fun="min", background=0)
plot(x)
You could also achieve this more "manually" with
r <- rast(v, res=2)
v <- sort(v, "ID", TRUE)
x <- rasterize(v, r, "ID")
x <- subst(x, NA, 0)
- rlanguageserver: Disable specific linter and avoid line-breaking
- Replace truncated word in a string
- Add label to right y axis in R ggplot
- Popup offset without using markers
- Multi marker depending a variable in a data frame
- Drawing numerous polygons
- Select icons or CircleMarkers based on factor variable
- Why is a map not displayed when I call leaflet in a function
- Display SpatialPolygonsDataFrame on leaflet map
- Leaflet in R not showing variable circle colors
- Overlay two population pyramid plots in one ggplot graph
- How to most efficiently filter a dataframe conditionally of values in another one, in the tidyverse framework?
- R group_by multiple conditionnal aggregation functions dplyr
- Conditionally split strings cells in a dataframe R
- How to have Shiny app render multiple plots
- Disable Shiny button conditionally to a rHandsontable format
- Is there a way to render dynamic headings in Quarto using R?
- What is the time complexity of sample?
- Add geom_text with x bar and round values
- Color-fill ggplot2 maps with a single color
- lapply a ggplot function on multiple inputs with consistent groups
- Is there a way to save text in R without it actively implementing \n?
- How to get number of unique observations using tbl_summary
- Using if_else in a mutate command yields an error about false having the wrong size
- Rotate x-axis value in dendogram using fviz_dend
- Export data frames in list to xlsx with named sheets
- Query all primary roads from PostGis Database filled with OSM Data within a boundary
- Write a list in r to an excel file
- Calculating net assets with missing values: discrepancy in results
- Error in curl::curl_fetch_memory(url, handle = handle) : Couldn't connect to server in R: oauth2.0_token())