rasterize simple feature geometries

st_rasterize(
  sf,
  template = guess_raster(sf, ...) %||% st_as_stars(st_bbox(sf), values = NA_real_,
    ...),
  file = tempfile(),
  driver = "GTiff",
  options = character(0),
  ...
)

Arguments

sf

object of class sf

template

stars object with desired target geometry

file

temporary file name

driver

driver for temporary file

options

character; options vector for GDALRasterize

...

arguments passed on to st_as_stars

Examples

demo(nc, echo = FALSE, ask = FALSE)
#> Reading layer `nc.gpkg' from data source `/home/runner/work/_temp/Library/sf/gpkg/nc.gpkg' using driver `GPKG' #> Simple feature collection with 100 features and 14 fields #> Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity #> Geometry type: MULTIPOLYGON #> Dimension: XY #> Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 #> Geodetic CRS: NAD27
(x = st_rasterize(nc)) # default grid:
#> stars object with 2 dimensions and 12 attributes #> attribute(s): #> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's #> AREA 0.042 0.108 0.142 0.1451932 0.181 0.241 30904 #> PERIMETER NA NA NA NaN NA NA 65001 #> CNTY_ NA NA NA NaN NA NA 65001 #> CNTY_ID NA NA NA NaN NA NA 65001 #> FIPSNO NA NA NA NaN NA NA 65001 #> CRESS_ID NA NA NA NaN NA NA 65001 #> BIR74 NA NA NA NaN NA NA 65001 #> SID74 NA NA NA NaN NA NA 65001 #> NWBIR74 NA NA NA NaN NA NA 65001 #> BIR79 NA NA NA NaN NA NA 65001 #> SID79 NA NA NA NaN NA NA 65001 #> NWBIR79 NA NA NA NaN NA NA 65001 #> dimension(s): #> from to offset delta refsys point values x/y #> x 1 461 -84.3239 0.0192484 NAD27 FALSE NULL [x] #> y 1 141 36.5896 -0.0192484 NAD27 FALSE NULL [y]
plot(x, axes = TRUE)
# a bit more customized grid: (x = st_rasterize(nc, st_as_stars(st_bbox(nc), nx = 100, ny = 50, values = NA_real_)))
#> stars object with 2 dimensions and 12 attributes #> attribute(s): #> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's #> AREA 0.042 0.108 0.142 0.145536 0.181 0.241 2373 #> PERIMETER NA NA NA NaN NA NA 5000 #> CNTY_ NA NA NA NaN NA NA 5000 #> CNTY_ID NA NA NA NaN NA NA 5000 #> FIPSNO NA NA NA NaN NA NA 5000 #> CRESS_ID NA NA NA NaN NA NA 5000 #> BIR74 NA NA NA NaN NA NA 5000 #> SID74 NA NA NA NaN NA NA 5000 #> NWBIR74 NA NA NA NaN NA NA 5000 #> BIR79 NA NA NA NaN NA NA 5000 #> SID79 NA NA NA NaN NA NA 5000 #> NWBIR79 NA NA NA NaN NA NA 5000 #> dimension(s): #> from to offset delta refsys point values x/y #> x 1 100 -84.3239 0.0886687 NAD27 FALSE NULL [x] #> y 1 50 36.5896 -0.0541531 NAD27 FALSE NULL [y]
plot(x, axes = TRUE)
(ls = st_sf(a = 1:2, st_sfc(st_linestring(rbind(c(0.1, 0), c(1.1, 1))), st_linestring(rbind(c(0, 0.05), c(1, 0.05))))))
#> Simple feature collection with 2 features and 1 field #> Geometry type: LINESTRING #> Dimension: XY #> Bounding box: xmin: 0 ymin: 0 xmax: 1.1 ymax: 1 #> CRS: NA #> a #> 1 1 #> 2 2 #> c..st_sfc.st_linestring.rbind.c.0.1..0...c.1.1..1.....st_linestring.rbind.c.0..........0.05...c.1..0.05...... #> 1 LINESTRING (0.1 0, 1.1 1) #> 2 LINESTRING (0 0.05, 1 0.05)
(grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = NA_real_))
#> stars object with 2 dimensions and 1 attribute #> attribute(s): #> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's #> values NA NA NA NaN NA NA 100 #> dimension(s): #> from to offset delta refsys point values x/y #> x 1 10 0 0.1 NA NA NULL [x] #> y 1 10 1 -0.1 NA NA NULL [y]
# Only the left-top corner is part of the grid cell: sf_extSoftVersion()["GDAL"]
#> GDAL #> "3.0.4"
plot(st_rasterize(ls, grd), axes = TRUE, reset = FALSE) # ALL_TOUCHED=FALSE;
plot(ls, add = TRUE, col = "red")
plot(st_rasterize(ls, grd, options = "ALL_TOUCHED=TRUE"), axes = TRUE, reset = FALSE)
plot(ls, add = TRUE, col = "red")
# add lines to existing 0 values, summing values in case of multiple lines: (grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = 0))
#> stars object with 2 dimensions and 1 attribute #> attribute(s): #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> values 0 0 0 0 0 0 #> dimension(s): #> from to offset delta refsys point values x/y #> x 1 10 0 0.1 NA NA NULL [x] #> y 1 10 1 -0.1 NA NA NULL [y]
r = st_rasterize(ls, grd, options = c("MERGE_ALG=ADD", "ALL_TOUCHED=TRUE")) plot(r, axes = TRUE, reset = FALSE)
plot(ls, add = TRUE, col = "red")