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

18

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

ingrese la descripción de la imagen aquí

Chad Cooper
fuente
Tengo una utilidad que utiliza subprocesos. Abrir para ejecutar múltiples traducciones de gdal al mismo tiempo que utilizo para extraer un gran ráster a mosaicos usando una red de pesca, particularmente útil si la entrada y / o salida está muy comprimida (por ejemplo, LZW o desinflar GeoTiff ), si ninguno de los dos está muy comprimido, el proceso alcanza un máximo en el acceso al HDD y no es mucho más rápido que ejecutar uno a la vez. Desafortunadamente, no es lo suficientemente genérico como para compartirlo debido a las rígidas convenciones de nomenclatura, pero de todos modos es motivo de reflexión. La opción -multi para GDALWarp a menudo causa problemas y solo usa 2 hilos (uno de lectura, uno de escritura) no todos disponibles.
Michael Stimson

Respuestas:

18

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

-srcwin xoff yoff xsize ysize: selecciona una subventana de la imagen de origen para copiar en función de la ubicación de píxel / línea.

-projwin ulx uly lrx lry: selecciona una subventana de la imagen de origen para copiar (como -srcwin) pero con las esquinas dadas en coordenadas georreferenciadas.

Tendría que encontrar las ubicaciones de píxeles / líneas o las coordenadas de las esquinas y luego recorrer los valores con gdal_translate. Algo así como el pitón rápido y sucio a continuación funcionará si usar valores de píxeles y -srcwin es adecuado para usted, será un poco más difícil de 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)
wwnick
fuente
1
Hola, cuando intento la opción -projwin con una imagen geotiff, recibo una advertencia que dice "Advertencia: Computado -srcwin -3005000 1879300 50 650 cae completamente fuera de la extensión de la trama. Sin embargo, continúo" No estoy seguro de dónde estoy haciendo mal parece que no utilizando sus coordenadas georreferenciadas.
ncelik
@ncelik probablemente se deba a que está utilizando coordenadas de celda en su projwin y debería utilizar srcwin en su lugar. Si tiene dificultades, publique una nueva pregunta con toda la información relevante para que podamos hacer sugerencias sobre su problema específico.
Michael Stimson
15

Mi solución, basada en la de @wwnick, lee las dimensiones ráster del archivo en sí y cubre toda la imagen haciendo que los mosaicos de bordes sean más pequeños 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)
Ries
fuente
Creo que debería ser sys.argv [1] donde dice sys.argv [2], ¿verdad?
oskarlin
3
Creo que sys.argv [2] se usa como prefijo para los archivos de salida. Super útil, gracias @Ries!
Charlie Hofmann
4

Hay un script de Python incluido específicamente para retilizar rásteres, 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

p.ej:

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

mikewatt
fuente
4

Para @Aaron que preguntó:

Espero encontrar una versión gdalwarp de la respuesta de @ wwnick que utilice la opción -multi para operaciones multinúcleo y multiproceso mejoradas

Descargo de responsabilidad leve

Esto utiliza gdalwarp, aunque no estoy totalmente convencido de que habrá mucho aumento de rendimiento. Hasta ahora he estado vinculado a E / S: ejecutar este script en un gran ráster cortándolo en muchas partes más pequeñas no parece intensivo en la CPU, por lo que supongo que el cuello de botella está escribiendo en el disco. Si planea volver a proyectar simultáneamente los mosaicos o algo similar, esto podría cambiar. Hay consejos de ajuste aquí . Una breve jugada no produjo ninguna mejora para mí, y la CPU nunca pareció ser el factor limitante.

Dejando de lado el descargo de responsabilidad, aquí hay un script que se usará gdalwarppara dividir un ráster en varios mosaicos más pequeños. Puede haber alguna pérdida debido a la división del piso, pero esto se puede solucionar eligiendo la cantidad de mosaicos que desee. Será n+1dónde nestá el número por el que divide para obtener las variables tile_widthy tile_height.

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))
RoperMaps
fuente
3

Puede usar r.tile de GRASS GIS. r.tile genera un mapa ráster separado para cada mosaico con nombres de mapas numerados basados ​​en el prefijo definido por el usuario. Se puede definir el ancho de los mosaicos (columnas) y el alto de los mosaicos (filas).

Al utilizar la API de Python de sesión de hierba, solo se necesitan unas pocas líneas de código de Python para llamar a la funcionalidad r.tile desde el exterior, es decir, para escribir un script independiente. Al usar r.external y r.external.out tampoco se produce duplicación de datos durante el paso de procesamiento de GRASS GIS.

Pseudocódigo:

  1. inicializar sesión de hierba
  2. definir formato de salida con r.external.out
  3. vincular archivo de entrada con r.external
  4. ejecute r.tile que genera los mosaicos en el formato definido anteriormente
  5. desvincular r.external.out
  6. cerrar sesión de hierba
markusN
fuente