Mehr

Wie extrahiert man Werte aus einem Raster basierend auf Lat und Lon?

Wie extrahiert man Werte aus einem Raster basierend auf Lat und Lon?


Ich habe mehrere netcdf-Dateien, die ich mit der netCDF4-Bibliothek in Python lesen kann. Die Daten sind in 4320 Zeilen und 2160 Spalten gerastert.

Ich habe auch eine Liste von Lat / Lon (in Dezimalgrad) und versuche, den Wert der Gitterzelle zu finden, in der die Koordinaten liegen.

Derzeit konvertiere ich die Koordinaten durch eine einfache lineare Transformation in Gitterzellengrenzen und frage die netcdf-Datei mit diesen neuen Punkten ab. Ist dieser Ansatz gültig?

# prelim from netCDF4 import Dataset import numpy as np # Koordinatendaten lons = [40.822651, 40.69685, 40.750038, 40.798931, 40.74612] lats = [19.83832, 20.007561, 19.97426, 19.86334, 19.843889] # Grenzen der Gitterzellen entsprechend lat/ lon finden = np.ceil(12*(lons + 180)) lats_grid = np.ceil(12*(lats + 90)) # Wert in Gitterzelle für gegebenen lat/lon finden fnc = Dataset(filename.nc, 'r') lat = fnc.variables['latitude'][:] lon = fnc.variables['longitude'][:] variable = fnc.variables['variable'][:] print variable[lons_grid, lats_grid]

Wenn nicht, gibt es eine Möglichkeit, Punktdaten mit einem Raster in Python zu schneiden, ohne arcpy zu verwenden? (Meine netcdf-Dateien sind sehr groß und können mit arcpy nicht gelesen werden.)


Wenn Sie die Pixel-/Zellengröße und den minimalen lon/Lat des Rasters kennen, ist dies eine einfache Berechnung.

px = (lon - minlon) / xcellsize py = (lat - minlat) / ycellsize

Sie können auch min/max lon/lat und die Anzahl der Spalten/Zeilen verwenden:

xcellsize = (maxlon - minlon) / ncols ycellsize = (maxlat - minlat) / nrows #dann verwenden Sie die vorherige Berechnung, um die Pixelkoordinaten zu erhalten

Um die entsprechenden Zellenkoordinaten aus Ihrem Breiten- und Längengrad zu erhalten, müssen Sie die Abdeckung der netCDF-Datei kennen (z. B. die Koordinate der oberen linken Zelle und die Größe der Zellen in x- und y-Richtung). @Lukes Lösung unten funktioniert für ein gerades x/y-Raster ohne Zellendrehung. Wenn Sie das nicht haben, können Sie sich die affin Bibliothek, um Breiten- und Längengrad in Zellkoordinaten umzuwandeln. Gegeben GDAL GeoTransform in der Form:

  1. Zelle oben linksxKoordinate
  2. xZellengröße in projizierten Einheiten
  3. xRichtungsdrehung (normalerweise 0)
  4. Zelle oben linksjaKoordinate
  5. jaRichtungsdrehung (normalerweise 0)
  6. NegativjaZellengröße in projizierten Einheiten

Also zum Beispiel die-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0würde als obere linke Zelle des Rasters mit Koordinaten übersetzen-237481.5, 237536.4, und ein425.0Einheitsquadrat ohne Drehung.

Verwendung deraffinBibliothek können Sie dies wie folgt in ein Transformationsobjekt umwandeln:

aus affinem Import Affine aff = Affine.from_gdal(-237481.5, 425,0, 0.0, 237536.4, 0.0, -425.0)

Welches kann Zellenkoordinaten in projizierte Koordinaten transformieren:

xs = np.array([0, 1, 2]) ys = np.array([3, 4, 5]) x_coords, y_coords = aff * (xs, ys)

Oder (in Ihrem Fall) von Koordinaten zurück zu Zellenkoordinaten:

xs, ys = ~aff * (np.array(lons), np.array(lats))

Diese Werte sind Gleitkommazahlen, daher müssen Sie sie in ganze Zahlen umwandeln, um Zellkoordinaten zu erhalten, die Sie verwenden können.

xs = np.round(xs).astype(np.int) ys = np.round(ys).astype(np.int)

Sie können diese dann als Indizes für Ihr netCDF4-Array verwenden (verwenden Sie die neueste Version von netCDF4 - 1.2.1, da dies bedeutet, dass Sie Duplikate nicht sortieren oder entfernen müssen).

variable = fnc.variables['variable'][xs, ys] # Achten Sie auf Ihre Indexreihenfolge

Dies gibt aufgrund der leicht unterschiedlichen netCDF-Indizierung ein quadratisches Array zurück, aber Sie können die tatsächlichen Werte, die Sie suchen, als Diagnose erhalten:

Werte = variabel.diagonal()

Schau das Video: mosaic raster dataset Landsat ArcGis