Raster neu klassifizieren mit Python, GDAL und Numpy

9

Ich möchte eine Rasterdatei von einem Raster mit 10 Klassen in ein Raster mit 8 Klassen unter Verwendung von Pyhton, GDAL und / oder Numpy umklassifizieren. Die Klassen werden als Ganzzahlen dargestellt. Ich habe versucht, die Schritte aus diesem Beitrag zu befolgen. Raster mit GDAL und Python , dem doc numpy.equal und dem Dokument gdal_calc neu klassifizieren. Allerdings ohne Erfolg.

Die neu zu klassifizierende Rasterdatei hat ganzzahlige Werte zwischen 0 und 11 und enthält auch die Werte 100 und 255. Im Folgenden wird die Neuklassifizierung (von Wert: zu Wert) dargestellt:

nodata: 4, 0: 4, 1: 1, 2: 2, 3: 3, 4: 3, 5: 4, 6: 5, 7: 5, 8: 6, 9: 7, 10: 8, 100: nodata, 255: nodata,

Was ich tun konnte, ist die Rasterdatei auszuwählen, die mit tkinter.FileDialog neu klassifiziert werden soll, und die Rasterinformationen wie Geotransform und Pixelgröße mit reclass = gdal.Open (Raster, GA_ReadOnly) abzurufen.

Wie löse ich das oben genannte?

Es könnte erwähnenswert sein, dass die zu klassifizierenden Raster in einigen Fällen ziemlich groß sein können (500 MB bis 5 GB).

PyMapr
quelle
Es gibt ein weiteres Beispiel auf dem GeoExamples Blog
Bennos
@bennos, habe das Skript im Blog ausprobiert, aber beim Entpacken des Arrays wird ein Speicherfehler zurückgegeben.
PyMapr
Ich schlage vor, Sie besprechen dieses Problem mit Roger Veciana i Rovira, dem Autor des Beitrags, da er seinen Code besser kennt als ich und vielleicht weiß, wie man das Problem
löst
Durch Ändern des Eingabe-Rasters von 16 Bit ohne Vorzeichen auf 8 Bit ohne Vorzeichen wurde das Speicherproblem behoben. Die Neuklassifizierung dauert jedoch ungefähr genauso lange wie die des folgenden Skripts dmh126.
PyMapr

Antworten:

6

Hier haben Sie ein einfaches Python-Skript zur Neuklassifizierung, ich habe es geschrieben und es funktioniert für mich:

from osgeo import gdal

driver = gdal.GetDriverByName('GTiff')
file = gdal.Open('/home/user/workspace/raster.tif')
band = file.GetRasterBand(1)
lista = band.ReadAsArray()

# reclassification
for j in  range(file.RasterXSize):
    for i in  range(file.RasterYSize):
        if lista[i,j] < 200:
            lista[i,j] = 1
        elif 200 < lista[i,j] < 400:
            lista[i,j] = 2
        elif 400 < lista[i,j] < 600:
            lista[i,j] = 3
        elif 600 < lista[i,j] < 800:
            lista[i,j] = 4
        else:
            lista[i,j] = 5

# create new file
file2 = driver.Create( 'raster2.tif', file.RasterXSize , file.RasterYSize , 1)
file2.GetRasterBand(1).WriteArray(lista)

# spatial ref system
proj = file.GetProjection()
georef = file.GetGeoTransform()
file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.FlushCache()

Ändern Sie einfach die Bereiche.

Ich hoffe es wird helfen.

dmh126
quelle
3
Sie sollten file2entweder mit del file2oder schließen file2 = None, um sicherzustellen, dass es auf die Festplatte geschrieben wird. .FlushCache()beeinflusst nur den internen Kachel-Cache von GDAL.
Kersten
@ dmh126, ich habe die Bereiche geändert und das Skript funktioniert. Es ist jedoch nicht sehr schnell (schnell umstritten). Das Eingabe-Raster war ungefähr 120 MB groß und dauerte ungefähr 15 Minuten. Mit Hilfe eines kommerziellen Fernerkundungspakets dauert es Sekunden. Empfehlungen zur Verkürzung der Bearbeitungszeit?
PyMapr
Ich denke, dass Multithreading helfen kann. Sie können versuchen, alle Ihre Kerne zu verwenden, es gibt eine Frage gis.stackexchange.com/questions/162978/…
dmh126
Es macht keinen Sinn, eine Double-for-Schleife zu verwenden, siehe Antwort unten
Mattijn
Richtig, die Neuklassifizierung von Doppelschleifen und Elementen ist die langsamste aller möglichen Möglichkeiten, dies zu tun. Verwenden Sie die leistungsstarken Teile von numpy wie ufuncs: docs.scipy.org/doc/numpy-1.10.1/reference/ufuncs.html .
Sgillies
16

Anstatt die von dmh126 beschriebene Neuklassifizierung als Double-for-Schleife durchzuführen, verwenden Sie Folgendes np.where:

# reclassification    
lista[np.where( lista < 200 )] = 1
lista[np.where((200 < lista) & (lista < 400)) ] = 2
lista[np.where((400 < lista) & (lista < 600)) ] = 3
lista[np.where((600 < lista) & (lista < 800)) ] = 4
lista[np.where( lista > 800 )] = 5

Auf einem Array von 6163 x 3537 Pixel (41,6 MB) erfolgt die Klassifizierung in 1,59 Sekunden, wobei die Verwendung der Double-for-Schleife 12 Minuten und 41 Sekunden dauert. Insgesamt nur eine Beschleunigung von 478x.

Fazit: Verwenden Sie niemals eine Double-for-Schleife numpy

Mattijn
quelle
Vielen Dank für den Hinweis, aber ich denke, das wird zu einem Problem führen, wenn sich die Eingabeklassen mit den Ausgabeklassen überschneiden. Ich möchte nicht, dass mein neuer Wert durch die nächste Regel geändert wird.
Etrimaille
@Gustry - Bin gerade auf dieses Problem gestoßen.
Relima
Überprüfen Sie also meine Antwort unten: gis.stackexchange.com/questions/163007/…
etrimaille
6

Hier ist ein grundlegendes Beispiel mit Rasterio und Numpy:

import rasterio as rio
import numpy as np


with rio.open('~/rasterio/tests/data/rgb1.tif') as src:
    # Read the raster into a (rows, cols, depth) array,
    # dstack this into a (depth, rows, cols) array,
    # the sum along the last axis (~= grayscale)
    grey = np.mean(np.dstack(src.read()), axis=2)

    # Read the file profile
    srcprof = src.profile.copy()

classes = 5
# Breaks is an array of the class breaks: [   0.   51.  102.  153.  204.]
breaks = (np.arange(classes) / float(classes)) * grey.max()

# classify the raster
classified = np.sum(np.dstack([(grey < b) for b in breaks]), axis=2).reshape(1, 400, 400).astype(np.int32)

# Update the file opts to one band
srcprof.update(count=1, nodata=None, dtype=classified.dtype)

with rio.open('/tmp/output.tif', 'w', **srcprof) as dst:
    # Write the output
    dst.write(classified)
Damon
quelle
2

Um die Antwort von @Mattijn zu vervollständigen, denke ich, dass dies zu einem Problem führen wird, wenn sich die Eingabeklassen mit den Ausgabeklassen überschneiden. Ich möchte nicht, dass mein neuer Wert durch die nächste Regel geändert wird.

Ich weiß nicht, ob ich an Geschwindigkeit verliere, aber ich sollte eine tiefe Kopie machen:

list_dest = lista.copy()

list_dest[np.where( lista < 0 )] = 0
list_dest[np.where((0 <= lista) & (lista <= 1)) ] = 1
list_dest[np.where((1 < lista) & (lista <= 5)) ] = 2
list_dest[np.where( 5 < lista )] = 3
Etrimaille
quelle
1

Hier ist ein weiterer Rasterio- Ansatz, den ich mithilfe des Rasterio-Kochbuchs und der Antwort von @ Mattijn zusammen gehackt habe .

import rasterio
import numpy as np

with rasterio.open('input_raster.tif') as src:    
    # Read as numpy array
    array = src.read()
    profile = src.profile

    # Reclassify
    array[np.where(array == 0)] = 4 
    array[np.where(array == 2)] = 1
    # and so on ...  

with rasterio.open('output_raster.tif', 'w', **profile) as dst:
    # Write to disk
    dst.write(array)
Aaron
quelle
0

In einigen Fällen kann die Numpy- Digitalisierung hilfreich sein, um schnell von Bereichen zu Behältern zu gelangen.

import rasterio
import numpy as np

with rasterio.open('my_raster.tif') as src:    
    array = src.read()
    profile = src.profile
    bins = np.array([-1.,-0.7,-0.4, 0.2, 1]) 
    inds = np.digitize(array, bins)

with rasterio.open('output_raster.tif', 'w', **profile) as dst:
    dst.write(inds)
Tactopoda
quelle
0

Mit Raster-RGB-Farbtabellenunterstützung:

import numpy as np
from osgeo import gdal

path_inDs = "/data/OCS_2016.extract.tif"
path_outDs = "/data/OCS_2016.postpython.tif"

driver = gdal.GetDriverByName('GTiff')
file = gdal.Open(path_inDs)

if file is None:
  print ('Could not open image file')
  sys.exit(1)

band = file.GetRasterBand(1)
lista = band.ReadAsArray()


# reclassification
lista[np.where(lista == 31)] = 16

# create new file
file2 = driver.Create(path_outDs, file.RasterXSize , file.RasterYSize , 1, gdal.GPI_RGB)
file2.GetRasterBand(1).WriteArray(lista)

# spatial ref system
proj = file.GetProjection()
georef = file.GetGeoTransform()
meta = file.GetMetadata()
colors = file.GetRasterBand(1).GetRasterColorTable()

file2.SetProjection(proj)
file2.SetGeoTransform(georef)
file2.SetMetadata(meta)
file2.GetRasterBand(1).SetRasterColorTable(colors)

file2.FlushCache()
del file2

Ser
quelle
0

Eine etwas andere Alternative könnte die folgende sein:

import numpy as np
from osgeo import gdal

original = gdal.Open('**PATH**\\origianl_raster.tif')



# read the original file

band = original.GetRasterBand(1) # assuming that the file has only one band
band_array = band.ReadAsArray()



#create a new array with reclassified values

new_array = np.where(band_array == np.nan, 4, 
                np.where(band_array == 0, 4,
                    np.where(band_array == 1, 1,
                        np.where(band_array == 2, 2,
                            np.where(band_array == 3, 3,
                                np.where(band_array == 4, 3, 
                                    np.where(band_array == 5, 4,
                                        np.where(band_array == 6, 5,
                                            np.where(band_array == 7, 5,
                                                np.where(band_array == 8, 6, 
                                                    np.where(band_array == 9, 7,
                                                       np.where(band_array == 10, 8,
                                                            np.where(band_array == 100, np.nan, np.nan))))))))))))) 
                                # the last line also includes values equal to 255, as they are the only ones left



# create and save reclassified raster as a new file

outDs = gdal.GetDriverByName('GTiff').Create("**PATH**\\reclassified_raster.tif", original.RasterXSize, original.RasterYSize, 1, gdal.GDT_Float32)

outBand = outDs.GetRasterBand(1)
outBand.WriteArray(new_array)

outDs.SetGeoTransform(original.GetGeoTransform())
outDs.SetProjection(original.GetProjection())


# flush cache

outDs.FlushCache()

Dieses Skript spielt mit numpy.where ( https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html ): in allen Schritten außer dem letzten, anstatt einen Wert zurückzugeben, wenn der Bedingung ist nicht erfüllt, es gibt eine andere np.where zurück. Und es geht so lange weiter, bis alle möglichen Werte des Rasters berücksichtigt sind.

Giallo
quelle