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('.') [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?