Finden Sie die minimale Begrenzungsausdehnung eines bestimmten Pixelwerts innerhalb des Rasters?

9

Ich frage mich, ob es eine Möglichkeit gibt, den Mindestbegrenzungsumfang für ein Raster mit einem bestimmten Wert zu ermitteln. Ich habe ein Raster aus einem globalen Bild ausgeschnitten und die Ausdehnung wird als globale Ausdehnung mit viel NoData-Bereich festgelegt. Ich möchte den NoData-Bereich aus diesem Raster entfernen und nur den Bereich beibehalten, der die Pixel des jeweiligen Werts enthält. Wie kann ich das machen?

Hier ist mein Beispiel: Ich möchte den Wert = 1 (blauer Bereich) extrahieren und die Ausdehnung des blauen Bereichs anstelle der ganzen Welt für die weitere Verarbeitung verwenden.

Beispielbild

Gesehen
quelle
Könnten Sie ein Beispiel posten?
Aaron
"Ich möchte die Nullzeilen und -spalten für dieses Raster löschen." Was genau bedeutet das? Ich verstehe nicht, was das gewünschte Endergebnis ist.
blah238
Suchen Sie nach "minimaler Begrenzungsausdehnung" nach einer rechteckigen Ausdehnung oder einem polygonalen Fußabdruck, der den Bereich des Bildes mit Daten darstellt?
blah238
1
@Tomek, die OP ist auf der Suche zu finden , das Ausmaß, müssen nicht manuell erstellen.
blah238
1
Wenn buchstäblich alles ein faires Spiel ist, verfügt einige Software über integrierte Befehle, um dies zu tun. Siehe zum Beispiel reference.wolfram.com/mathematica/ref/ImageCrop.html .
whuber

Antworten:

6

WENN ich die Frage richtig verstanden habe, klingt es so, als ob Sie den minimalen Begrenzungsrahmen der Werte kennen möchten, die nicht null sind. Möglicherweise können Sie das Raster in Polygone konvertieren, die gewünschten Polygone auswählen und sie dann wieder in ein Raster konvertieren. Sie können sich dann die Eigenschaftswerte ansehen, die Ihnen den Minium-Begrenzungsrahmen geben sollen.

Dango
quelle
1
Alles in allem ist dies angesichts der Grenzen des ArcGIS-Rasterverarbeitungsworkflows wahrscheinlich der beste Ansatz.
blah238
Ich habe das genau gemacht. Gibt es einen automatischen Weg? Ich denke, der Raster-zu-Polygon-Algorithmus hat einen Schritt, um den minimalen Begrenzungsrahmen des Rasters zu extrahieren.
Gesehen am
bist du nach einer Python-Lösung?
Dango
8

Der Trick besteht darin, die Grenzen der Daten mit Werten zu berechnen. Der vielleicht schnellste, natürlichste und allgemeinste Weg, diese zu erhalten, sind zonale Zusammenfassungen: Wenn Sie alle Nicht-NoData-Zellen für die Zone verwenden, liefern die zonalen Min- und Max-Werte der Gitter, die die X- und Y-Koordinaten enthalten, die volle Ausdehnung.

ESRI ändert ständig die Art und Weise, wie diese Berechnungen durchgeführt werden können. Beispielsweise wurden integrierte Ausdrücke für die Koordinatengitter mit ArcGIS 8 gelöscht und scheinen nicht zurückgekehrt zu sein. Nur zum Spaß, hier ist eine Reihe von schnellen, einfachen Berechnungen, die die Arbeit erledigen, egal was passiert.

  1. Konvertieren Sie das Raster in eine einzelne Zone, indem Sie es wie in sich selbst mit sich selbst gleichsetzen

    "Mein Gitter" == "Mein Gitter"

  2. Erstellen Sie ein Spaltenindexgitter, indem Sie ein konstantes Gitter mit dem Wert 1 im Fluss akkumulieren. (Die Indizes beginnen mit 0.) Wenn gewünscht, multiplizieren Sie dies mit der Zellengröße und addieren Sie die x-Koordinate des Ursprungs, um ein x-Koordinatengitter zu erhalten. " X "(unten gezeigt).

  3. Erstellen Sie auf ähnliche Weise ein Zeilenindexgitter ( und dann ein y-Koordinatengitter "Y"), indem Sie ein konstantes Gitter mit dem Wert 64 im Fluss akkumulieren.

  4. Verwenden Sie das Zonenraster aus Schritt (1), um die zonalen Min und Max von "X" und "Y" zu berechnen : Sie haben jetzt die gewünschte Ausdehnung.

Endgültiges Bild

(Das Ausmaß, wie in den beiden Tabellen der Zonenstatistik gezeigt, ist in dieser Figur durch einen rechteckigen Umriss dargestellt. Das Gitter "I" ist das in Schritt (1) erhaltene Zonenraster.)

Um weiter zu gehen, müssen Sie diese vier Zahlen aus ihren Ausgabetabellen extrahieren und sie verwenden, um den Analyseumfang zu begrenzen. Das Kopieren des ursprünglichen Rasters mit dem eingeschränkten Analyseumfang schließt die Aufgabe ab.

whuber
quelle
8

Hier ist eine Version der @ whubers-Methode für ArcGIS 10.1+ als Python-Toolbox (.pyt).

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Raster Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [ClipNoData]


class ClipNoData(object):
    def __init__(self):
        """Clip raster extent to the data that have values"""
        self.label = "Clip NoData"
        self.description = "Clip raster extent to the data that have values. "
        self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        params = []

        # First parameter
        params+=[arcpy.Parameter(
            displayName="Input Raster",
            name="in_raster",
            datatype='GPRasterLayer',
            parameterType="Required",
            direction="Input")
        ]

        # Second parameter
        params+=[arcpy.Parameter(
            displayName="Output Raster",
            name="out_raster",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Output")
        ]

        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return arcpy.CheckExtension('spatial')==u'Available'

    def execute(self, parameters, messages):
        """See https://gis.stackexchange.com/a/55150/2856
           ##Code comments paraphrased from @whubers GIS StackExchange answer
        """
        try:
            #Setup
            arcpy.CheckOutExtension('spatial')
            from arcpy.sa import *
            in_raster = parameters[0].valueAsText
            out_raster = parameters[1].valueAsText

            dsc=arcpy.Describe(in_raster)
            xmin=dsc.extent.XMin
            ymin=dsc.extent.YMin
            mx=dsc.meanCellWidth
            my=dsc.meanCellHeight
            arcpy.env.extent=in_raster
            arcpy.env.cellSize=in_raster
            arcpy.AddMessage(out_raster)

            ## 1. Convert the grid into a single zone by equating it with itself
            arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
            zones = Raster(in_raster) == Raster(in_raster)

            ## 2. Create a column index grid by flow-accumulating a constant grid
            ##    with value 1. (The indexes will start with 0.) Multiply this by
            ##    the cellsize and add the x-coordinate of the origin to obtain
            ##    an x-coordinate grid.
            arcpy.AddMessage(r'Create a constant grid...')
            const = CreateConstantRaster(1)

            arcpy.AddMessage(r'2. Create an x-coordinate grid...')
            xmap = (FlowAccumulation(const)) * mx + xmin

            ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
            ##    constant grid with value 64.
            arcpy.AddMessage(r'3. Create a y-coordinate grid...')
            ymap = (FlowAccumulation(const * 64)) * my + ymin

            ## 4. Use the zone grid from step (1) to compute the zonal min and
            ##    max of "X" and "Y"
            arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

            xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORY\xrange", "NODATA", "MIN_MAX")
            xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORY\xrange", ["MIN","MAX"]).next()

            yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORY\yrange", "NODATA", "MIN_MAX")
            ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORY\yrange", ["MIN","MAX"]).next()

            arcpy.Delete_management(r"IN_MEMORY\xrange")
            arcpy.Delete_management(r"IN_MEMORY\yrange")

            # Write out the reduced raster
            arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
            output = Raster(in_raster) * 1
            output.save(out_raster)

        except:raise
        finally:arcpy.CheckInExtension('spatial')
user2856
quelle
Sehr netter Luke. In sich geschlossen, ausführbar, verwendet in_memory und zum Booten gut kommentiert. Ich musste die Hintergrundverarbeitung deaktivieren ( Geoverarbeitung> Optionen> ... ), damit sie funktioniert.
Matt Wilkie
1
Ich habe das Skript aktualisiert und canRunInBackground = False gesetzt. Ich werde bemerken, dass es sich lohnt, die Arbeitsbereichs- / Arbeitsbereichsumgebungen in einen lokalen Ordner (nicht FGDB) zu ändern, wie ich festgestellt habe, als ich sie als Standard (dh <Netzwerkbenutzerprofil> \ Default.gdb) belassen habe. Das Skript dauerte 9 Minuten !!! auf einem 250x250-Zellenraster laufen. Das Wechseln zu einer lokalen FGDB dauerte 9 Sekunden und zu einem lokalen Ordner 1-2 Sekunden ...
user2856
Ein guter Punkt zu lokalen Ordnern und vielen Dank für die schnelle Hintergrundkorrektur (viel besser als Anweisungen für alle zu schreiben, an die ich sie weitergebe). Könnte es wert sein, dies auf bitbucket / github / gcode / etc. Zu werfen.
Matt Wilkie
+1 Danke für diesen Beitrag, Luke! Ich weiß es zu schätzen, dass Sie die (ziemlich große) Lücke in meiner Antwort ausgefüllt haben.
whuber
4

Ich habe eine gdal- und numpy-basierte Lösung entwickelt. Es zerlegt die Rastermatrix in Zeilen und Spalten und löscht alle leeren Zeilen / Spalten. In dieser Implementierung ist "leer" weniger als 1, und nur einzelne Bandraster werden berücksichtigt.

(Während ich schreibe, wird mir klar, dass dieser Scanline-Ansatz nur für Bilder mit Nodata-Halsbändern geeignet ist. Wenn Ihre Daten Inseln in Meeren von Nullen sind, wird auch der Abstand zwischen den Inseln fallen gelassen, wodurch alles zusammengedrückt wird und die Georeferenzierung völlig durcheinander gebracht wird .)

Die Geschäftsteile (müssen ausgearbeitet werden, funktionieren nicht so wie sie sind):

    #read raster into a numpy array
    data = np.array(gdal.Open(src_raster).ReadAsArray())
    #scan for data
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
        # assumes data is any value greater than zero
    crop_box = (min(non_empty_rows), max(non_empty_rows),
        min(non_empty_columns), max(non_empty_columns))

    # retrieve source geo reference info
    georef = raster.GetGeoTransform()
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    # write to disk
    band = out_raster.GetRasterBand(1)
    band.WriteArray(new_data)
    band.FlushCache()
    out_raster = None

In einem vollständigen Skript:

import os
import sys
import numpy as np
from osgeo import gdal

if len(sys.argv) < 2:
    print '\n{} [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
    sys.exit(1)

src_raster = sys.argv[1]
out_raster = sys.argv[2]

def main(src_raster):
    raster = gdal.Open(src_raster)

    # Read georeferencing, oriented from top-left
    # ref:GDAL Tutorial, Getting Dataset Information
    georef = raster.GetGeoTransform()
    print '\nSource raster (geo units):'
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]
    cols, rows = raster.RasterYSize, raster.RasterXSize
    print '  Origin (top left): {:10}, {:10}'.format(xmin, ymax)
    print '  Pixel size (x,-y): {:10}, {:10}'.format(xcell, ycell)
    print '  Columns, rows    : {:10}, {:10}'.format(cols, rows)

    # Transfer to numpy and scan for data
    # oriented from bottom-left
    data = np.array(raster.ReadAsArray())
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
    crop_box = (min(non_empty_rows), max(non_empty_rows),
        min(non_empty_columns), max(non_empty_columns))

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    new_rows, new_cols = new_data.shape # note: inverted relative to geo units
    #print cropped_transform

    print '\nCrop box (pixel units):', crop_box
    print '  Stripped columns : {:10}'.format(cols - new_cols)
    print '  Stripped rows    : {:10}'.format(rows - new_rows)

    print '\nCropped raster (geo units):'
    print '  Origin (top left): {:10}, {:10}'.format(new_xmin, new_ymax)
    print '  Columns, rows    : {:10}, {:10}'.format(new_cols, new_rows)

    raster = None
    return new_data, cropped_transform


def write_raster(template, array, transform, filename):
    '''Create a new raster from an array.

        template = raster dataset to copy projection info from
        array = numpy array of a raster
        transform = geo referencing (x,y origin and pixel dimensions)
        filename = path to output image (will be overwritten)
    '''
    template = gdal.Open(template)
    driver = template.GetDriver()
    rows,cols = array.shape
    out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
    out_raster.SetGeoTransform(transform)
    out_raster.SetProjection(template.GetProjection())
    band = out_raster.GetRasterBand(1)
    band.WriteArray(array)
    band.FlushCache()
    out_raster = None
    template = None

if __name__ == '__main__':
    cropped_raster, cropped_transform = main(src_raster)
    write_raster(src_raster, cropped_raster, cropped_transform, out_raster)

Das Skript befindet sich in meinem Code-Stash auf Github, wenn der Link ein bisschen herumjagt; Diese Ordner sind reif für eine Neuorganisation.

matt wilkie
quelle
1
Dies funktioniert nicht für wirklich große Datensätze. Ich bekomme einenMemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError
user32882
2

Bei aller analytischen Leistungsfähigkeit fehlen in ArcGIS grundlegende Rastermanipulationen, die Sie mit herkömmlichen Desktop-Bildeditoren wie GIMP finden können . Er erwartet , dass Sie die gleiche Analyse Umfang für Ihre Ausgabe - Raster als Eingabe - Raster verwendet werden soll, es sei denn , Sie manuell die außer Kraft setzen Ausgabeausdehnung Umgebungseinstellung. Da dies genau das ist, was Sie suchen und nicht festlegen möchten, steht die ArcGIS-Vorgehensweise im Weg.

Trotz meiner Bemühungen und ohne auf die Programmierung zurückzugreifen, konnte ich keine Möglichkeit finden, das Ausmaß Ihrer gewünschten Teilmenge des Bildes zu ermitteln (ohne die rechnerisch verschwenderische Konvertierung von Raster in Vektor).

Stattdessen habe ich GIMP verwendet, um den blauen Bereich mit dem Farbauswahlwerkzeug auszuwählen, und dann die Auswahl invertiert, auf Löschen geklickt, um die restlichen Pixel zu entfernen, die Auswahl erneut invertiert, das Bild auf Auswahl zugeschnitten und es schließlich wieder exportiert PNG. GIMP hat es als 1-Bit-Tiefenbild gespeichert. Das Ergebnis ist unten:

Ausgabe

Da Ihrem Beispielbild eine räumliche Referenzkomponente fehlte und GIMP räumlich nicht bekannt ist, ist das Ausgabebild natürlich ungefähr so ​​nützlich wie Ihre Beispieleingabe. Sie müssen es georeferenzieren, damit es für eine räumliche Analyse von Nutzen ist.

blah238
quelle
Eigentlich Diese Operation verwendete einfach von Spatial Analyst in früheren Versionen werden: der zonale max und min der beiden Koordinatenraster (X und Y), mit der Funktion als die Zone, gibt das Ausmaß genau. (Nun, vielleicht möchten Sie es um eine halbe Zellengröße in alle vier Richtungen erweitern.) In ArcGIS 10 müssen Sie kreativ sein (oder Python verwenden), um ein Koordinatengitter zu erstellen. Egal, das Ganze kann vollständig innerhalb SA erfolgt nur Netzbetrieb mit und ohne manuelle Eingriffe.
whuber
@whuber Ich habe Ihre Lösung woanders gesehen, bin mir aber immer noch nicht sicher, wie ich Ihre Methode implementieren kann. Können Sie mir etwas mehr Details dazu zeigen?
Gesehen am
@Seen Bei einer schnellen Suche auf dieser Website wird unter gis.stackexchange.com/a/13467 ein Konto für die Methode gefunden .
whuber
1

Hier ist eine Möglichkeit, SAGA GIS zu verwenden: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

In SAGA GIS gibt es ein "Crop to Data" -Modul (in der Grid Tools-Modulbibliothek), das die Aufgabe ausführt.

Dazu müssten Sie jedoch Ihr Geotif mit dem GDAL-Importmodul importieren, in SAGA verarbeiten und schließlich mit dem GDAL-Exportmodul erneut als Geotif exportieren.

Eine weitere Möglichkeit nur ArcGIS GP - Tools wäre einen TIN von Ihrem Raster mit bauen Raster zu TIN berechne seine Grenze mit TIN - Domäne , und Clip durch die Grenze (oder dessen Umschlag mit Raster - Funktion Umschlag zu Polygon ).

blah238
quelle