3

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

0

2 Answers 2

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") )
1

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'

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.