read and write geotiffs bigger than memory?

Question

The geotiff format is a versatile and easy to understand file format. The only limitation is when I try to do raster calculations I prefer to use GDAL (over ArcPy) which comes with the limitation that you read the data in the computer's memory.

My objective is to do raster calculations such as setting all values below a threshold to a NoData value.

Somebody shared a python library called DASK with me. Can somebody explain how to open, calculate and write geotiff files bigger than your memory?

The problem is how GDAL stores the array in the computer memory

try:
    from osgeo import ogr, osr, gdal
except:
    sys.exit('ERROR: cannot find GDAL/OGR modules')
# Enable GDAL/OGR exceptions
gdal.UseExceptions()
import os
import numpy as np
import dask.array as da

INPUTPATH = os.path.join('Path to store the geotiffs')

OUTPUTPATH = os.path.join('Path to store the output geotiffs')


def readFile(filename):
    filehandle = gdal.Open(filename)
    band1 = filehandle.GetRasterBand(1)
    geotransform = filehandle.GetGeoTransform()
    geoproj = filehandle.GetProjection()
    Z = band1.ReadAsArray() # Can this be avoided?
    xsize = filehandle.RasterXSize
    ysize = filehandle.RasterYSize
    filehandle = None
    return xsize,ysize,geotransform,geoproj,Z

def writeFile(filename,geotransform,geoprojection,data):
    (x,y) = data.shape
    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    # you can change the dataformat but be sure to be able to store negative values including -9999
    dst_datatype = gdal.GDT_Float32
    dst_ds = driver.Create(filename,y,x,1,dst_datatype, [ 'COMPRESS=LZW' ])
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(geoprojection)
    dst_ds = None
    return 1

files = os.listdir(INPUTPATH)

for oneFile in files:
    print oneFile
    InputBaseName = oneFile.split('.')[0]
    [xsize, ysize, geotransform, geoproj, Z] = readFile(os.path.join(INPUTPATH, oneFile))
    Z[Z<-9000] = -9999
    writefilename = os.path.join(OUTPUTPATH,(InputBaseName+".tif"))
    writeFile(writefilename, geotransform, geoproj, Z)

print "Done"

This link is useful but I can't seem to find a way to read the array from the geotiff in dask. Any tips/tricks?

Answer 1

You might be able to do your processing with the gdal_calc.py utility, which has block-level processing (not whole file processing). For example, this command should do the same thing your program does, using Numpy's where command:

gdal_calc.py -A input1.tif --outfile=result.tif --calc="where(A < -9000, -9999, A)" \
             --co=COMPRESS=LZW --NoDataValue=-9999 --type=Float32

(This is one command wrapped on two lines. Windows users via OSGeo4W Shell will need to replace gdal_calc.py with gdal_calc and replace continue line character \ with ^)

However, if you need to do something more specific, you can see how the source code for gdal_calc.py works.