Erhalten Sie mit Python gdal Koordinaten und entsprechende Pixelwerte von GeoTiff und speichern Sie sie als Numpy-Array

10

Wie kann ich projizierte Koordinaten sowie die tatsächlichen Pixelwerte an diesen Koordinaten aus einer GeoTiff-Datei abrufen und sie dann in einem Numpy-Array speichern? Ich habe die Datei arsenci020l.tif und ihre Koordinaten sind in Metern. Unten ist die gekürzte Ausgabe von gdalinfo, die ich darauf ausgeführt habe.

~$ gdalinfo arsenci020l.tif 
Driver: GTiff/GeoTIFF
Files: arsenci020l.tif
       arsenci020l.tfw
Size is 10366, 7273
Coordinate System is:
PROJCS["Lambert Azimuthal Equal Area projection with arbitrary plane grid; projection center 100.0 degrees W, 45.0 degrees N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-6086629.000000000000000,4488761.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
...

Es gab hier eine ähnliche Frage zum Abrufen von Lat / Long-Koordinaten von tiff (Erhalten von Latitude und Longitude aus einer GeoTIFF-Datei), und die Antwort zeigte, wie nur die oberen linken x- und y-Pixelkoordinaten erhalten werden. Ich muss ALLE projizierten Pixelkoordinaten erhalten sowie die Pixelwerte erhalten und sie in einem numpy-Array speichern. Wie kann ich es tun?

irakhman
quelle
Sie möchten 10366 × 7273 = über 75 Millionen Punkte?
Mike T
@MikeT Ich denke schon, ich weiß nicht wirklich, wie ich das Problem, das ich zu lösen versuche, besser lösen kann: Ich muss jedem Schwerpunkt des US-Blocks die nächstgelegene Pixelkoordinate aus diesem Datensatz finden und dann die zuweisen Entsprechender Pixelwert für diesen Block. Beim Durchsuchen wurde mir klar, dass die cKDTree-Abfrage mir bei der Suche nach dem nächsten Nachbarn helfen wird. Die Python-Funktion für den Algorithmus fordert einen "Baum" an, der als numpy-Array abgefragt werden soll. Um einen "Baum" zu erstellen. Von allen Pixelkoordinaten aus diesem Datensatz muss ich sie alle irgendwie speichern. Wenn Sie eine bessere Lösung haben, lassen Sie es mich bitte wissen!
irakhman

Antworten:

7

würde als Kommentar hinzufügen, aber ein bisschen lang - falls Sie gdal / ogr in Python verwenden möchten - könnte so etwas funktionieren (zusammen gehackt von einem anderen Code, den ich hatte - nicht getestet!). Dies setzt auch voraus, dass es nicht den nächsten findet Rasterpixel auf einen Polygonschwerpunkt, Sie fragen das Raster einfach am xy des Schwerpunkts ab. Ich habe keine Ahnung, wie hoch der Geschwindigkeitskompromiss sein könnte ...

from osgeo import gdal,ogr

fc='PathtoYourVector'
rast='pathToYourRaster'

def GetCentroidValue(fc,rast):
    #open vector layer
    drv=ogr.GetDriverByName('ESRI Shapefile') #assuming shapefile?
    ds=drv.Open(fc,True) #open for editing
    lyr=ds.GetLayer(0)

    #open raster layer
    src_ds=gdal.Open(rast) 
    gt=src_ds.GetGeoTransform()
    rb=src_ds.GetRasterBand(1)
    gdal.UseExceptions() #so it doesn't print to screen everytime point is outside grid

    for feat in lyr:
        geom=feat.GetGeometryRef()
        mx=geom.Centroid().GetX()
        my=geom.Centroid().GetY()

        px = int((mx - gt[0]) / gt[1]) #x pixel
        py = int((my - gt[3]) / gt[5]) #y pixel
        try: #in case raster isnt full extent
            structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_Float32) #Assumes 32 bit int- 'float'
            intval = struct.unpack('f' , structval) #assume float
            val=intval[0]
        except:
            val=-9999 #or some value to indicate a fail

       feat.SetField('YOURFIELD',val)
       lyr.SetFeature(feat)

    src_ds=None
    ds=None

GetCentroidValue(fc,rast)
flüssige Bewegung
quelle
14

Das sollte dich zum Laufen bringen. Die Rasterwerte werden mit Rasterio gelesen , und die Pixelmittenkoordinaten werden mit affinen Werten in Ost- / Nordwerte konvertiert , die dann mit pyproj in Breiten- / Längengrade konvertiert werden . Die meisten Arrays haben dieselbe Form wie das Eingabe-Raster.

import rasterio
import numpy as np
from affine import Affine
from pyproj import Proj, transform

fname = '/path/to/your/raster.tif'

# Read raster
with rasterio.open(fname) as r:
    T0 = r.transform  # upper-left pixel corner affine transform
    p1 = Proj(r.crs)
    A = r.read()  # pixel values

# All rows and columns
cols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))

# Get affine transform for pixel centres
T1 = T0 * Affine.translation(0.5, 0.5)
# Function to convert pixel row/column index (from 0) to easting/northing at centre
rc2en = lambda r, c: (c, r) * T1

# All eastings and northings (there is probably a faster way to do this)
eastings, northings = np.vectorize(rc2en, otypes=[np.float, np.float])(rows, cols)

# Project all longitudes, latitudes
p2 = Proj(proj='latlong',datum='WGS84')
longs, lats = transform(p1, p2, eastings, northings)
Mike T.
quelle
1
Wenn ich dies verwende, erhalte ich die Meldung "AttributeError: 'DatasetReader' Objekt hat kein Attribut 'affine'" für die Zeile "T0 = r.affine"
Mitchus
@mitchus Anscheinend affineist es nur ein Alias ​​für transformund der Alias ​​wurde aus der neuesten Version von Rasterio entfernt. Ich habe die Antwort bearbeitet, aber es sieht so aus, als müsste sie einer Peer-Review unterzogen werden, da ich neu hier bin. :)
Autumnsault
1
Es sieht auch so aus, als wären die Indizes falsch A.shape, was nur zwei Dimensionen hat.
Autumnsault