I have a GPKG file with a few vector layers. I have a few processing steps in my script and every now and then a new layer is added to the file. This works fine for vector layers using something like
sf::st_write(vectorlayer, "output.gpkg", "layername", update = TRUE)
but I have no idea how to add raster layers. I can use raster::writeRaster to create a new GPKG, but I'm not able to extend an existing one. Any ideas what I'm missing here?
Short answer: you can't. Geopackage's default raster capabilities are fairly constrained - see https://gdal.org/drivers/raster/gpkg.html. The design of these constraints is very much oriented towards handling imagery - byte data type only (integer values 0-255), no more than four bands which are assumed to be RGB/A etc.
There is also a GPKG extension for handling tiled gridded coverage data e.g. elevation data, but its still constrained in its own way - you can only write a single band of data to a single GPKG, and it uses either 16 bit PNG tiles for integer data or 32 bit TIFF tiles with optional LZW compression for floating-point data. This makes sense; many-banded rasters together with sets of embedded overview tiles could easily blow out to enormous file sizes and the internal table structure would become inefficient and likely unstable.
Geopackages may be SQLite database under the hood, but they were really intended as a snapshot/transfer format, not a space for working datasets as you describe. You might want to consider NetCDF/HDF5 for this kind of task.
Not sure about raster but I did it in terra with this code
writeRaster(
tif,
"test.gpkg",
filetype = "GPKG",
gdal = c("APPEND_SUBDATASET=YES", "RASTER_TABLE=layer_name"),
overwrite = TRUE
)
