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!
You can do that with
terra::rasterizeYour polygons
Solution, use
terra::rasterizewith argumentsfunandbackgroundYou could also achieve this more "manually" with