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.
Antworten:
Wenn Sie gerne R verwenden , gibt es ein Paket namens raster . Sie können ein Raster mit dem folgenden Befehl einlesen:
Wenn Sie es dann ansehen (durch Eingabe
test
), sehen Sie die folgenden Informationen: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.
quelle
R
, können Sie StandardfunktionenR
oder diegetValues
Methode verwenden, um auf Zellenwerte zuzugreifen. Von dort ist es einfach, die höchsten Werte und ihre Positionen zu identifizieren.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.
quelle
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
argsort
andere. 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).Zeigt Folgendes für meine Test-Raster-Datei an:
quelle
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 bestehtargsort
darin, die NODATA-Werte aus der Analyse zu entfernen oder sie einfach in der for-Schleife zu überspringen, wenn sie NaN / maskiert sind.