Create new output raster dataset with CPL_rasterize/gdal_utils('rasterize', ...)?

[johnbaums]

Since GDAL 1.8.0, gdal_rasterize has been able to create the output dataset. This doesn't seem to occur with gdal_utils('rasterize', ...).

Looking at CPL_rasterize, I see that args such as -tr are currently unsupported. Would it be possible for a future version to support raster dataset creation? (I know this is possible with gdalUtils, but would love to have the functionality available in the self-contained tools provided by sf).

[edzer]

Thanks for getting back to this. I think I now repaired gdal_utils('rasterize', ...) to also handle the case where the destination dataset is not given (as file), but defined by options like e.g. -tr and -te. We can now do:

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
demo(nc, ask = FALSE, echo = FALSE)
# Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64-pc-linux-gnu-library/3.6/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
# bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
# epsg (SRID):    4267
# proj4string:    +proj=longlat +datum=NAD27 +no_defs

# gdal_utils('rasterize') with predefined grid:

(x = st_rasterize(nc))
# stars object with 2 dimensions and 1 attribute
# attribute(s):
#      AREA       
#  Min.   :0.042  
#  1st Qu.:0.108  
#  Median :0.142  
#  Mean   :0.145  
#  3rd Qu.:0.181  
#  Max.   :0.241  
#  NA's   :30904  
# dimension(s):
#   from  to   offset      delta                       refsys point values    
# x    1 461 -84.3239  0.0192484 +proj=longlat +datum=NAD2... FALSE   NULL [x]
# y    1 141  36.5896 -0.0192484 +proj=longlat +datum=NAD2... FALSE   NULL [y]
write_stars(x, "out.tif")
datasource
# [1] "/home/edzer/R/x86_64-pc-linux-gnu-library/3.6/sf/gpkg/nc.gpkg"
gdal_utils('rasterize', datasource, "out.tif", options = c("-b", "1", "-a", "AREA"))

# gdal_utils('rasterize') without predefined grid:
e = as.character(st_bbox(nc))
gdal_utils('rasterize', datasource, "out2.tif", options = c("-tr", ".01", ".01", "-a", "AREA", "-te", e))

[johnbaums]

Great, thanks @edzer. Seems to be working well for me. I appreciate the quick fix.