¿leer y escribir geotiffs más grandes que la memoria?

Question

El formato geotiff es un formato de archivo versátil y fácil de entender. La única limitación es que cuando intento hacer cálculos raster prefiero usar GDAL (en lugar de ArcPy) que viene con la limitación de que lees los datos en la memoria del ordenador.

Mi objetivo es realizar cálculos de rasterización como, por ejemplo, establecer todos los valores por debajo de un umbral en un valor NoData.

Alguien compartió una biblioteca python llamada DASK conmigo. ¿Alguien puede explicar cómo abrir, calcular y escribir archivos geotiff más grandes que tu memoria?

El problema es cómo GDAL almacena la matriz en la memoria del ordenador

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"

Este enlace es útil, pero no puedo encontrar una manera de leer la matriz de la geotiff en dask. ¿Algún consejo o truco?

Answer 1

Es posible que pueda realizar el tratamiento con la función gdal_calc.py que procesa bloques (no archivos enteros). Por ejemplo, este comando debería hacer lo mismo que hace su programa, utilizando el comando de Numpy donde mando:

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

(Esto es un comando envuelto en dos líneas. Los usuarios de Windows a través de OSGeo4W Shell tendrán que sustituir gdal_calc.py con gdal_calc y reemplazar el carácter de línea continua \ con ^ )

Sin embargo, si necesita hacer algo más específico, puede vea cómo el código fuente de gdal_calc.py funciona .