¿Dividir ráster en trozos más pequeños utilizando GDAL?

16

Tengo un raster (USGS DEM en realidad) y necesito dividirlo en partes más pequeñas como muestra la imagen de abajo. Eso se logró en ArcGIS 10.0 usando la herramienta Dividir ráster. Me gustaría un método de software libre 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 la trama y decir qué tan grande (4KM por 4KM trozos) me gustaría dividir en.

    
pregunta Chad Cooper 16.09.2011 - 22:22

5 respuestas

18

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

  

-srcwin xoff yoff xsize ysize:   Selecciona una subventana de la imagen de origen para copiar según la ubicación del 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 georreferenciadas   coordenadas.

Necesitará encontrar las ubicaciones de píxel / línea o coordenadas de esquina y luego recorrer los valores con gdal_translate. Algo así como el pitón rápido y sucio que se muestra a continuación funcionará si el uso de valores de píxel y -srcwin es adecuado para usted, será un poco más de trabajo resolverlos con las 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)
    
respondido por el wwnick 17.09.2011 - 01:13
15

Mi solución, basada en la de @wwnick lee las dimensiones de la trama desde el propio archivo, y cubre toda la imagen haciendo que los mosaicos de borde 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)
    
respondido por el Ries 27.03.2014 - 12:33
3

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

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

Pseudo código:

  1. inicializar sesión de césped
  2. defina el formato de salida con r.external.out
  3. enlace de archivo de entrada con r.externo
  4. ejecutar r.tile que genera los mosaicos en el formato definido anteriormente
  5. desvincular r.external.out
  6. cerrar sesión de hierba
respondido por el markusN 14.10.2018 - 23:34
3

Para @Aaron quien 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 usa gdalwarp , aunque no estoy totalmente convencido de que habrá mucho aumento de rendimiento. Hasta ahora he estado enlazado a la E / S: ejecutar este script en un ráster grande y cortarlo en muchas partes más pequeñas no parece un uso intensivo de 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 me produjo ninguna mejora, y la CPU nunca pareció ser el factor limitante.

Dejando de lado el descargo de responsabilidad, aquí hay un script que usará gdalwarp para 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 puede solucionarse seleccionando el número de fichas que desee. Será n+1 donde n es el número por el que se divide para obtener las variables tile_width y 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))
    
respondido por el RoperMaps 14.10.2018 - 23:22
2

Hay un script de Python incluido específicamente para volver a compilar 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

por ejemplo:

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

    
respondido por el gberard 19.10.2018 - 21:22

Lea otras preguntas en las etiquetas