¿Dividir una trama en trozos más pequeños con GDAL?

Question

Tengo un raster (USGS DEM en realidad) y tengo que dividirlo en trozos más pequeños como la imagen de abajo muestra. Eso se logró en ArcGIS 10.0 utilizando la herramienta Split Raster. Me gustaría un método FOSS para hacer esto. He mirado en GDAL, pensando que seguramente lo haría (de alguna manera con gdal_translate), pero no puedo encontrar nada. En última instancia, me gustaría ser capaz de tomar la trama y decir qué tan grande (4KM por 4KM trozos) me gustaría que se divide en.

Answer 1

gdal_translate funcionará utilizando las opciones -srcwin o -projwin.

-srcwin xoff yoff xsize ysize: Selecciona una subventana de la imagen de origen para copiarla basándose en la ubicación del píxel/línea.

-projwin ulx uly lrx lry: Selecciona una subventana de la imagen de origen para copiarla (como -srcwin), pero con las esquinas georreferenciadas. georreferenciadas.

Tendría que obtener las ubicaciones de píxel/línea o las coordenadas de esquina y luego hacer un bucle sobre los valores con gdal_translate. Algo como el python rápido y sucio a continuación funcionará si el uso de valores de píxeles y -srcwin es adecuado para usted, será un poco más de trabajo para resolver con coordenadas.

import os, gdal
from gdalconst import *

width = 512
height = 512
tilesize = 64

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(tilesize)+", " \
            +str(tilesize)+" utm.tif utm_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)

Answer 2

Hay un script python incluido específicamente para retiling rasters, gdal_retile :

gdal_retile.py [-v] [-co NAME=VALUE]* [-of out_format] [-ps pixelWidth pixelHeight]
               [-overlap val_in_pixel]
               [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
                      CInt16/CInt32/CFloat32/CFloat64}]'
               [ -tileIndex tileIndexName [-tileIndexField tileIndexFieldName]]
               [ -csv fileName [-csvDelim delimiter]]
               [-s_srs srs_def]  [-pyramidOnly]
               [-r {near/bilinear/cubic/cubicspline/lanczos}]
               -levels numberoflevels
               [-useDirForEachRow]
               -targetDir TileDirectory input_files

Por ejemplo

gdal_retile.py -ps 512 512 -targetDir C:\example\dir some_dem.tif

Answer 3

Mi solución, basada en la de @wwnick lee las dimensiones raster del propio archivo, y cubre toda la imagen haciendo más pequeños los tiles de los bordes si es necesario:

import os, sys
from osgeo import gdal

dset = gdal.Open(sys.argv[1])

width = dset.RasterXSize
height = dset.RasterYSize

print width, 'x', height

tilesize = 5000

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        w = min(i+tilesize, width) - i
        h = min(j+tilesize, height) - j
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(w)+", " \
            +str(h)+" " + sys.argv[1] + " " + sys.argv[2] + "_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)

Answer 4

Puede utilizar r.azulejo de GRASS GIS. r.tile genera un mapa raster separado para cada tile con nombres de mapa numerados basados en el prefijo definido por el usuario. Se puede definir el ancho de los mosaicos (columnas) y la altura de los mosaicos (filas).

Utilización de la sesión sobre la hierba Python API sólo se necesitan unas pocas líneas de código Python para llamar a la funcionalidad de r.tile desde fuera, es decir, para escribir un script independiente. Uso de r.externo y r.external.out tampoco se produce duplicación de datos durante el paso de procesamiento de GRASS GIS.

Pseudocódigo:

1. inicializar grass-session
2. definir el formato de salida con r.external.out
3. vincular el archivo de entrada con r.external
4. ejecuta r.tile, que genera los mosaicos en el formato definido anteriormente
5. unlink r.external.out
6. cerrar la sesión sobre el césped

Answer 5

Para @Aaron que preguntó:

Espero encontrar una versión gdalwarp de la respuesta de @wwnick que utilice la opción -multi para mejorar las operaciones multinúcleo y multihilo.

Ligera exención de responsabilidad

Esto utiliza gdalwarp aunque no estoy totalmente convencido de que vaya a aumentar mucho el rendimiento. Hasta ahora he estado limitado por la E/S - ejecutar este script en un raster grande cortándolo en muchas partes más pequeñas no parece intensivo para la CPU, así que asumo que el cuello de botella es la escritura en disco. Si estás planeando reproyectar simultáneamente los mosaicos o algo similar, entonces esto podría cambiar. Hay consejos de ajuste aquí . Una breve partida no me produjo ninguna mejora, y la CPU nunca pareció ser el factor limitante.

Descargo de responsabilidad aparte, aquí hay un script que utilizará gdalwarp para dividir una trama en varios mosaicos más pequeños. Es posible que se produzcan algunas pérdidas debido a la división en pisos, pero esto puede solucionarse eligiendo el número de mosaicos que desee. Será n+1 donde n es el número por el que se divide para obtener el tile_width y tile_height variables.

import subprocess
import gdal
import sys

def gdalwarp(*args):
    return subprocess.check_call(['gdalwarp'] + list(args))

src_path = sys.argv[1]
ds = gdal.Open(src_path)

try:
    out_base = sys.argv[2]
except IndexError:
    out_base = '/tmp/test_'

gt = ds.GetGeoTransform()

width_px = ds.RasterXSize
height_px = ds.RasterYSize

# Get coords for lower left corner
xmin = int(gt[0])
xmax = int(gt[0] + (gt[1] * width_px))

# get coords for upper right corner
if gt[5] > 0:
    ymin = int(gt[3] - (gt[5] * height_px))
else:
    ymin = int(gt[3] + (gt[5] * height_px))

ymax = int(gt[3])

# split height and width into four - i.e. this will produce 25 tiles
tile_width = (xmax - xmin) // 4
tile_height = (ymax - ymin) // 4

for x in range(xmin, xmax, tile_width):
    for y in range(ymin, ymax, tile_height):
        gdalwarp('-te', str(x), str(y), str(x + tile_width),
                 str(y + tile_height), '-multi', '-wo', 'NUM_THREADS=ALL_CPUS',
                 '-wm', '500', src_path, out_base + '{}_{}.tif'.format(x, y))