Dividindo a varredura em pedaços menores usando o GDAL?

18

Eu tenho uma varredura (USGS DEM na verdade) e preciso dividi-la em partes menores, como mostra a imagem abaixo. Isso foi realizado no ArcGIS 10.0 usando a ferramenta Split Raster. Eu gostaria de um método FOSS para fazer isso. Eu olhei para o GDAL, pensando que certamente o faria (de alguma forma com gdal_translate), mas não consigo encontrar nada. Por fim, eu gostaria de poder pegar a varredura e dizer qual o tamanho (pedaços de 4KM por 4KM) que eu gostaria que fosse dividido.

insira a descrição da imagem aqui

Chad Cooper
fonte
Oopen para executar várias traduções gdal ao mesmo tempo que eu uso para extrair uma varredura grande para blocos usando uma rede de pesca, particularmente útil se a entrada e / ou saída for altamente compactada (por exemplo, LZW ou desinflar o GeoTiff ), se nenhum deles estiver altamente compactado, o processo atingirá o pico no acesso ao disco rígido e não será muito mais rápido do que executar um de cada vez. Infelizmente, não é genérico o suficiente para compartilhar devido a convenções rígidas de nomenclatura, mas é algo que deve ser pensado de qualquer maneira. A opção -multi para GDALWarp geralmente causa problemas e usa apenas 2 threads (uma leitura, uma gravação) nem todos disponíveis.
Michael Stimson

Respostas:

18

O gdal_translate funcionará usando as opções -srcwin ou -projwin.

-srcwin xoff yoff xsize ysize: seleciona uma sub-janela da imagem de origem para copiar com base na localização dos pixels / linhas.

-projwin ulx uly lrx lry: seleciona uma sub-janela da imagem de origem para cópia (como -srcwin), mas com os cantos fornecidos nas coordenadas georreferenciadas.

Você precisaria apresentar os locais de pixel / linha ou coordenadas de canto e, em seguida, percorrer os valores com gdal_translate. Algo como o python rápido e sujo abaixo funcionará se o uso de valores de pixel e -srcwin for adequado para você, será um pouco mais trabalhoso para resolver as 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
fonte
1
Olá, quando tento a opção -projwin com uma imagem geotiff, recebo um aviso dizendo "O aviso: Computado -srcwin -3005000 1879300 50 650 fica completamente fora do alcance da varredura. Continuando no entanto" Não sei ao certo onde estou fazendo errado, parece que não. usando suas coordenadas georreferenciadas.
Ncelik 26/10/16
@ncelik isso é provavelmente porque você está usando coordenadas de célula no seu projwin e deveria estar usando srcwin. Se estiver com dificuldades, envie uma nova pergunta com todas as informações relevantes para que possamos fazer sugestões sobre o seu problema específico.
Michael Stimson
15

Minha solução, baseada na do @wwnick, lê as dimensões de varredura do próprio arquivo e cobre toda a imagem, tornando os blocos de borda menores, se necessário:

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
fonte
Eu acho que deveria ser sys.argv [1] onde diz sys.argv [2], certo?
Oskarlin 30/07
3
sys.argv [2] é usado como um prefixo para os arquivos de saída, acredito. Super útil-- obrigado @Ries!
93017 Charlie Hofmann
4

Existe um script python incluído especificamente para rasters de reciclagem , 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 exemplo:

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

mikewatt
fonte
4

Para @Aaron, que perguntou:

Espero encontrar uma versão gdalwarp da resposta do @ wwnick que utilize a opção -multi para operações multicore e multithread avançadas

Isenção de responsabilidade

Isso usa gdalwarp, embora eu não esteja totalmente convencido de que haverá muito ganho de desempenho. Até agora, eu estou vinculado à E / S - executar esse script em uma grande varredura, cortando-o em muitas partes menores, não parece muito com a CPU, portanto, presumo que o gargalo esteja gravando no disco. Se você planeja reprojetar simultaneamente os blocos ou algo semelhante, isso pode mudar. Existem dicas de ajuste aqui . Uma breve jogada não resultou em nenhuma melhoria para mim e a CPU nunca pareceu ser o fator limitante.

Isenção de responsabilidade, aqui está um script que será usado gdalwarppara dividir uma varredura em vários blocos menores. Pode haver alguma perda devido à divisão do piso, mas isso pode ser resolvido escolhendo o número de peças que você deseja. Será n+1onde nestá o número que você divide para obter as variáveis tile_widthe 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
fonte
3

Você pode usar r.tile do GRASS GIS. O r.tile gera um mapa de varredura separado para cada bloco com nomes de mapas numerados com base no prefixo definido pelo usuário. A largura dos ladrilhos (colunas) e a altura dos ladrilhos (linhas) podem ser definidas.

Usando a API Python da sessão de grama, apenas algumas linhas de código Python são necessárias para chamar a funcionalidade r.tile de fora, ou seja, para escrever um script independente. Usando r.external e r.external.out também não ocorre duplicação de dados durante a etapa de processamento do GRASS GIS.

Pseudo-código:

  1. inicializar sessão de grama
  2. defina o formato de saída com r.external.out
  3. vincular arquivo de entrada com r.external
  4. execute r.tile que gera os blocos no formato definido acima
  5. desvincular r.external.out
  6. fechar sessão de grama
markusN
fonte