gdal_calc raster calculation fails on 28 GB GeoTiff

Question

I am trying to do a raster calculations on a 28 GB GeoTiff using gdal_calc.py (3.0.4) but the process fails

gdal_calc.py -A DEM_0.tif --outfile=DEM_1.tif --calc="A+0.33" --format=GTiff --type=Float32 --co="COMPRESS=LZW"

when the size of the output file becomes ~650 MB the command keeps running but the file stops increasing in size. It happens if I run the process on the SSD internally but also on an external SSD. Furthermore, I noticed that the rate of increase of the filesize is only 70MB / min (which is rather slow I think, although gdal_calc only uses one CPU).

After a comment below I also created another geotiff with tilesize 64x64. I ran gdal_calc.py overnight and pretty much the same thing happened - it was still running after more than 7 hours so I killed the process.

shortened gdalinfo:

Size is 146140, 86076

Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND

Band 1 Block=256x256 Type=Float32, ColorInterp=Gray NoData Value=-32767 Overviews: 73070x43038, 36535x21519, 18268x10760, 9134x5380, 4567x2690, 2284x1345, 1142x673, 571x337, 286x169, 143x85

Unit Type: metre

Answer 1

Although it isn't a true gdal_calc answer it solves the issue of working with a 28 GB raster using GRASS through R.

library(sf)
librart(raster)
library(rgrass7)

dem  <- raster('DEM.tif) 
bb   <- st_bbox(dem)

initGRASS(gisBase = "/usr/lib/grass78", 
          home = '/home/user/tmp', # to have sufficient disk space
          gisDbase = '/home/user/tmp', 
          location = "wadden", 
          mapset = "PERMANENT",
          override = TRUE)

execGRASS("g.proj", flags = c("c", "quiet"), 
          proj4 = st_crs(dem)$proj4string)

execGRASS("g.region", flags = c("quiet"), 
          n = as.character(bb["ymax"]), s = as.character(bb["ymin"]), 
          e = as.character(bb["xmax"]), w = as.character(bb["xmin"]), 
          res = paste(res(dem)[1])) 

execGRASS("r.in.gdal", flags=c("o", "overwrite"), parameters=list(input='DEM.tif'), output="dem"))
execGRASS("r.mapcalc", expression="dem_fixed = dem + 0.33")
execGRASS("r.out.gdal", parameters=list(
                          input="dem_fixed", 
                          output="dem_fixed.tif"),
                          format="GTiff",
                          type="Float32",
                          createopt=c("COMPRESS=LZW",  "BIGTIFF=YES")),
                        flags=c("f", "overwrite") )

Answer 2

The issue might have been caused by an older version of gdal_calc.py Versions starting from 3.3 have useful new options, so it is good to try with the newest gdal again. Options like TILED=YES or COMPRESS=DEFLATE can be added with the syntax: --co="TILED=YES" one can define many options like this so I use always --co='COMPRESS=DEFLATE' --co='PREDICTOR=YES' --co='BIGTIFF=IF_SAFER' --co='TILED=YES'