Diviser le raster en petits morceaux à l'aide de GDAL?

18

J'ai un raster (DEM USGS en fait) et je dois le diviser en petits morceaux comme le montre l'image ci-dessous. Cela a été accompli dans ArcGIS 10.0 à l'aide de l'outil Fractionner le raster. Je voudrais une méthode FOSS pour ce faire. J'ai regardé GDAL, pensant sûrement qu'il le ferait (en quelque sorte avec gdal_translate), mais je ne trouve rien. En fin de compte, j'aimerais pouvoir prendre le raster et dire dans quelle taille (morceaux de 4KM par 4KM) je voudrais qu'il soit divisé.

entrez la description de l'image ici

Chad Cooper
la source
J'ai un utilitaire qui utilise un sous-processus.Popen pour exécuter plusieurs traductions gdal en même temps que j'utilise pour extraire un grand raster en tuiles à l'aide d'un résille, particulièrement utile si l'entrée et / ou la sortie est très compressée (par exemple LZW ou dégonfler GeoTiff ), si aucun des deux n'est fortement compressé, le processus atteint un pic sur l'accès au disque dur et n'est pas beaucoup plus rapide que l'exécution un par un. Malheureusement, ce n'est pas assez générique pour partager en raison de conventions de dénomination rigides, mais de quoi réfléchir. L'option -multi pour GDALWarp cause souvent des problèmes et n'utilise que 2 threads (une lecture, une écriture) qui ne sont pas tous disponibles.
Michael Stimson

Réponses:

18

gdal_translate fonctionnera en utilisant les options -srcwin ou -projwin.

-srcwin xoff yoff xsize ysize: sélectionne une sous-fenêtre de l'image source pour la copie en fonction de l'emplacement des pixels / lignes.

-projwin ulx uly lrx lry: sélectionne une sous-fenêtre de l'image source pour la copie (comme -srcwin) mais avec les coins donnés en coordonnées géoréférencées.

Vous devrez trouver les emplacements des pixels / lignes ou les coordonnées des coins, puis parcourir les valeurs avec gdal_translate. Quelque chose comme le python rapide et sale ci-dessous fonctionnera si l'utilisation de valeurs de pixels et -srcwin vous convient, sera un peu plus de travail pour trier les coordonnées.

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
la source
1
Salut quand j'essaie l'option -projwin avec une image géotiff, je reçois un avertissement disant "Avertissement: -srcwin -3005000 1879300 50 650 tombe complètement en dehors de l'étendue du raster. En continuant cependant" Je ne suis pas sûr où je me trompe semble ne pas en utilisant ses coordonnées géoréférencées.
ncelik
@ncelik, c'est probablement parce que vous utilisez des coordonnées de cellule dans votre projet et que vous devriez utiliser srcwin à la place. Si vous rencontrez des difficultés, veuillez poster une nouvelle question avec toutes les informations pertinentes afin que nous puissions faire des suggestions sur votre problème spécifique.
Michael Stimson
15

Ma solution, basée sur celle de @wwnick, lit les dimensions du raster à partir du fichier lui-même et couvre toute l'image en réduisant les tuiles de bord si nécessaire:

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
la source
Je pense que ce devrait être sys.argv [1] où il est écrit sys.argv [2], non?
oskarlin du
3
sys.argv [2] est utilisé comme préfixe pour les fichiers de sortie, je crois. Super utile - merci @Ries!
Charlie Hofmann
4

Il y a un script python fourni spécifiquement pour retiler les 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

par exemple:

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

mikewatt
la source
4

Pour @Aaron qui a demandé:

J'espère trouver une version gdalwarp de la réponse de @ wwnick qui utilise l'option -multi pour des opérations multicœurs et multithread améliorées

Légère exclusion de responsabilité

Cela utilise gdalwarp, même si je ne suis pas totalement convaincu qu'il y aura beaucoup de gain de performances. Jusqu'à présent, je suis lié aux E / S - l'exécution de ce script sur un grand raster le découpant en de nombreuses parties plus petites ne semble pas nécessiter beaucoup de CPU, donc je suppose que le goulot d'étranglement écrit sur le disque. Si vous prévoyez de re-projeter simultanément les tuiles ou quelque chose de similaire, cela pourrait changer. Il y a des conseils de réglage ici . Un bref jeu ne m'a apporté aucune amélioration, et le CPU n'a jamais semblé être le facteur limitant.

Mise à part la mise en garde, voici un script qui servira gdalwarpà diviser un raster en plusieurs tuiles plus petites. Il peut y avoir une perte due à la division du sol, mais cela peut être résolu en choisissant le nombre de carreaux que vous souhaitez. Ce sera n+1nest le nombre que vous divisez par pour obtenir les variables tile_widthet 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
la source
3

Vous pouvez utiliser r.tile de GRASS GIS. r.tile génère une carte raster distincte pour chaque tuile avec des noms de carte numérotés en fonction du préfixe défini par l'utilisateur. La largeur des carreaux (colonnes) et la hauteur des carreaux (rangées) peuvent être définies.

En utilisant l' API Python grass-session, seules quelques lignes de code Python sont nécessaires pour appeler la fonctionnalité r.tile de l'extérieur, c'est-à-dire pour écrire un script autonome. En utilisant r.external et r.external.out, aucune duplication de données ne se produit également pendant l'étape de traitement de GRASS GIS.

Pseudo code:

  1. initialiser grass-session
  2. définir le format de sortie avec r.external.out
  3. lier le fichier d'entrée avec r.external
  4. exécutez r.tile qui génère les tuiles dans le format défini ci-dessus
  5. dissocier r.external.out
  6. close grass-session
markusN
la source