Suchen von Positionen mit den höchsten Werten in Raster mithilfe von ArcGIS Desktop?

12

In ArcGIS 10 habe ich ein Raster, in dem ich das Pixel mit dem Maximalwert im Raster suchen und dessen Position (Pixelmitte) in Dezimalgrad zurückgeben möchte. Ich möchte diesen Prozess durchlaufen und den Ort des zweithöchsten Werts des Rasters, dann den dritten usw. zurückgeben, damit ich am Ende eine Liste von N Orten habe, die die höchsten Werte im Raster in der richtigen Reihenfolge haben.

Ich stelle mir vor, dass dies am einfachsten mit einem Python-Skript möglich ist, bin aber offen für andere Ideen, wenn es einen besseren Weg gibt.

mga
quelle
Haben Sie versucht, das Raster in Punkte umzuwandeln, dann X-, Y-Felder hinzuzufügen und zu sortieren?
Jakub Sisak GeoGraphics
Sind die Rasterwerte Gleitkommazahlen oder Ganzzahlen?
Whuber
@ Jakub - Nein, habe ich nicht. Ich werde mich wahrscheinlich nur für die obersten 1% der Punkte interessieren, daher weiß ich nicht, ob es sich lohnt, x, y-Felder für alle Punkte hinzuzufügen und dann zu sortieren. Vielleicht, wenn es keine effizientere Option gibt?
mga
@whuber - Die Rasterwerte sind Floats.
mga
@mga, es ist einen Versuch wert. Die Konvertierung ist sehr schnell und das Hinzufügen von XY ist ebenfalls ein Standardwerkzeug. Das Löschen unerwünschter Datensätze ist unkompliziert und kann in einem einzigen Modell zusammengefasst werden. Nur eine Idee.
Jakub Sisak GeoGraphics

Antworten:

5

Wenn Sie gerne R verwenden , gibt es ein Paket namens raster . Sie können ein Raster mit dem folgenden Befehl einlesen:

install.packages('raster')
library(raster)
test <- raster('F:/myraster')

Wenn Sie es dann ansehen (durch Eingabe test), sehen Sie die folgenden Informationen:

class       : RasterLayer 
dimensions  : 494, 427, 210938  (nrow, ncol, ncell)
resolution  : 200, 200  (x, y)
extent      : 1022155, 1107555, 1220237, 1319037  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
values      : F:/myraster 
min value   : 0 
max value   : 1 

Möglicherweise gibt es bessere Möglichkeiten, das Raster zu bearbeiten. Eine Möglichkeit, die gewünschten Informationen zu finden, besteht darin, den höchsten Wert zu ermitteln und die Matrixposition abzurufen und diese dann zu den niedrigeren Werten hinzuzufügen.

djq
quelle
1
+1 Für diese Referenz. Nachdem Sie das Raster eingelesen haben R, können Sie Standardfunktionen Roder die getValuesMethode verwenden, um auf Zellenwerte zuzugreifen. Von dort ist es einfach, die höchsten Werte und ihre Positionen zu identifizieren.
Whuber
1
Dank Ihrer Empfehlung habe ich das getan. Die Verwendung des Raster-Pakets in R war im Vergleich zum Testen in ArcGIS ein Kinderspiel. Ich habe auch andere räumliche Analysen in R verwendet und war mit den Ergebnissen sehr zufrieden. Guter Rat!
mga
8

Die Antwort erhalten Sie von indem ein Indikatorraster der oberen 1% der Werte mit Rastern für Breite und Länge kombiniert wird . Der Trick besteht darin, dieses Indikatorraster zu erstellen, da ArcGIS (immer noch! Nach 40 Jahren!) Kein Verfahren zum Rangieren von Rasterdaten hat.

Eine Lösung für Gleitkomma-Raster ist iterativ, aber barmherzig schnell . Sei n die Anzahl der Datenzellen. Die empirische kumulative Werteverteilung besteht aus allen Paaren (z, n (z)), wobei z ein Wert im Raster ist und n (z) die Anzahl der Zellen im Raster mit Werten kleiner oder gleich z ist . Wir erhalten eine Kurve, die (-infinity, 0) mit (+ infinity, n) verbindet, aus der Folge dieser Scheitelpunkte, die nach z geordnet sind . Es definiert dabei eine Funktion f , wobei (z, f (z)) immer auf der Kurve liegt. Sie möchten einen Punkt (z0, 0,99 * n) auf dieser Kurve finden.

Mit anderen Worten, die Aufgabe besteht darin, eine Null von f (z) - (1-0,01) * n zu finden . Tun Sie dies mit einer beliebigen Nullfindungsroutine (die beliebige Funktionen handhaben kann: Diese ist nicht differenzierbar). Das einfachste und häufig effizienteste ist das Raten und Prüfen: Zunächst wissen Sie, dass z0 zwischen dem Mindestwert zMin und dem Höchstwert zMax liegt. Errate einen vernünftigen Wert genau zwischen diesen beiden. Wenn die Schätzung zu niedrig ist, setzen Sie zMin = z0; ansonsten setze zMax = z0. Wiederholen Sie jetzt. Sie werden sich schnell der Lösung annähern. Sie sind nah genug, wenn zMax und zMin nah genug sind. Um konservativ zu sein, wählen Sie den endgültigen Wert von zMin als Lösung: Es werden möglicherweise einige zusätzliche Punkte gesammelt, die Sie später verwerfen können. Weitere Informationen finden Sie in Kapitel 9 der Numerischen Rezepte (der link geht zu einer älteren kostenlosen version).

Rückblickend auf diesen Algorithmus zeigt sich, dass Sie nur zwei Arten von Rasteroperationen ausführen müssen : (1) Wählen Sie alle Zellen aus, die kleiner oder gleich einem Zielwert sind, und (2) zählen Sie ausgewählte Zellen. Diese gehören zu den einfachsten und schnellsten Operationen, die es gibt. (2) kann als Zonenzählung oder durch Lesen eines Datensatzes aus der Attributtabelle des Auswahlrasters erhalten werden.

whuber
quelle
7

Ich habe dies vor einiger Zeit getan, obwohl meine Lösung GDAL verwendet (dies gilt also nicht nur für ArcGIS). Ich denke, Sie können ein NumPy-Array aus einem Raster in ArcGIS 10 abrufen, aber ich weiß es nicht genau. NumPy bietet eine einfache und leistungsstarke Array-Indizierung, wie auch argsortandere. In diesem Beispiel werden keine NODATA- oder Transformationskoordinaten von projiziert nach lat / long verarbeitet (dies ist jedoch mit osgeo.osr, das mit GDAL bereitgestellt wird, nicht schwierig).

import numpy as np
from osgeo import gdal

# Open raster file, and get GeoTransform
rast_src = gdal.Open(rast_fname)
rast_gt = rast_src.GetGeoTransform()

def get_xy(r, c):
    '''Get (x, y) raster centre coordinate at row, column'''
    x0, dx, rx, y0, ry, dy = rast_gt
    return(x0 + r*dx + dx/2.0, y0 + c*dy + dy/2.0)

# Get first raster band
rast_band = rast_src.GetRasterBand(1)

# Retrieve as NumPy array to do the serious work
rast = rast_band.ReadAsArray()

# Sort raster pixels from highest to lowest
sorted_ind = rast.argsort(axis=None)[::-1]

# Show highest top 10 values
for ind in sorted_ind[:10]:
    # Get row, column for index
    r, c = np.unravel_index(ind, rast.shape)
    # Get [projected] X and Y coordinates
    x, y = get_xy(r, c)
    print('[%3i, %3i] (%.3f, %.3f) = %.3f'%
          (r, c, x, y, rast[r, c]))

Zeigt Folgendes für meine Test-Raster-Datei an:

[467, 169] (2813700.000, 6353100.000) = 844.538
[467, 168] (2813700.000, 6353200.000) = 841.067
[469, 168] (2813900.000, 6353200.000) = 840.705
[468, 168] (2813800.000, 6353200.000) = 840.192
[470, 167] (2814000.000, 6353300.000) = 837.063
[468, 169] (2813800.000, 6353100.000) = 837.063
[482, 166] (2815200.000, 6353400.000) = 833.038
[469, 167] (2813900.000, 6353300.000) = 832.825
[451, 181] (2812100.000, 6351900.000) = 828.064
[469, 169] (2813900.000, 6353100.000) = 827.514
Mike T
quelle
+1 Danke, dass du das geteilt hast. Ich sehe die Unfähigkeit, mit NoData umzugehen, nicht als Einschränkung: Konvertieren Sie einfach alle NoData-Werte in extrem negative Werte, bevor Sie fortfahren. Beachten Sie auch, dass sich die Antworten bei einer erneuten Projektion des Rasters wahrscheinlich aufgrund einer erneuten Abtastung des Rasters ändern. Daher möchte man normalerweise nicht, dass die erneute Projektion während einer solchen Berechnung automatisch erfolgt. Stattdessen können die gemeldeten Koordinaten anschließend neu projiziert werden. Somit ist Ihre Lösung vollkommen allgemein.
whuber
Die NODATA-Behandlung kann implementiert werden, indem zuerst der Wert aus dem Raster abgerufen wird NODATA = rast_band.GetNoDataValue()und dann entweder ein NaN-Wert ( rast[rast == NODATA] = np.nan) oder ein maskiertes Array ( rast = np.ma.array(rast, mask=(rast == NODATA))) verwendet wird. Der kompliziertere Trick besteht argsortdarin, die NODATA-Werte aus der Analyse zu entfernen oder sie einfach in der for-Schleife zu überspringen, wenn sie NaN / maskiert sind.
Mike T