2

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?

  • 2
    No sane raster library would require the entire raster be in memory. This seems more likely to be an issue of incorrect use. Please format your Python code so it could execute. – Vince Feb 28 '17 at 0:31
  • 3
    Read the image block by block or line by line... look at band1.ReadRaster gdal.org/gdal_tutorial.html... I've used this for multi-terrabyte rasters on a line-by-line basis successfully. Just remember that the blocks are in cell coordinates from the upper left and not world coordinates. – Michael Stimson Feb 28 '17 at 0:41
  • done. My raster is a float 64 on 30s degree resolution 21600,43200 – RutgerH Feb 28 '17 at 0:42
  • 1
    But you're writing GDT_Float32, you can't write a double array to a float memory block... you would need to convert each element from double to float before writing or change your datatype to GDT_Float64 on creation. If you have arcpy you could have done all this with a Con(Float statement - assuming you have spatial analyst extension. – Michael Stimson Feb 28 '17 at 0:45
  • Thank you for these comments. I am still struggling how to perform algebra (in python IDE) on these scanlines or blocks (getting lost in GDAL docs), can you point me in the right direction? – RutgerH Mar 8 '17 at 19:47
5

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.

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.