Dapatkan koordinat dan nilai piksel yang sesuai dari GeoTiff menggunakan python gdal dan simpan sebagai array numpy

10

Bagaimana saya bisa mendapatkan koordinat yang diproyeksikan serta nilai piksel aktual pada koordinat tersebut dari file GeoTiff dan kemudian menyimpannya ke dalam array numpy? Saya punya file arsenci020l.tif, dan koordinatnya dalam meter. Di bawah ini adalah keluaran singkat dari gdalinfo yang saya jalankan di sana.

~$ 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)
...

Ada pertanyaan serupa di sini tentang mendapatkan koordinat lat / long dari tiff (Dapatkan Latitude dan Longitude dari File GeoTIFF) dan jawabannya menunjukkan bagaimana mendapatkan hanya koordinat x dan y piksel kiri atas. Saya perlu mendapatkan SEMUA koordinat piksel yang diproyeksikan serta mendapatkan nilai-nilai piksel dan menyimpannya dalam array numpy. Bagaimana saya bisa melakukannya?

irakhman
sumber
Anda ingin 10366 × 7273 = lebih dari 75 juta poin?
Mike T
@ MikeT Saya pikir begitu, saya tidak benar-benar tahu solusi yang lebih baik tentang bagaimana mendekati masalah yang saya coba pecahkan: Saya perlu mencari koordinat piksel terdekat dari dataset ini untuk setiap centroid blok AS dan kemudian menetapkan nilai pixel yang sesuai untuk blok itu. Dari pencarian sekitar saya menyadari bahwa permintaan cKDTree akan membantu saya dengan pencarian tetangga terdekat. Fungsi phyton untuk algoritma meminta "pohon" untuk query sebagai array numpy. Untuk membuat "pohon" dari semua koordinat piksel dari dataset ini, saya perlu menyimpan semuanya entah bagaimana. Jika Anda memiliki solusi yang lebih baik, tolong beri tahu saya!
irakhman

Jawaban:

7

akan menambahkan sebagai komentar, tetapi agak lama - jika Anda ingin menggunakan gdal / ogr dalam python - sesuatu seperti ini mungkin bekerja (diretas bersama dari beberapa kode lain yang saya miliki - tidak diuji!) Ini juga mengasumsikan bahwa daripada mencari yang terdekat pixel raster ke centroid poligon, Anda cukup menanyakan raster pada xy centroid. saya tidak tahu apa pengorbanan kecepatan mungkin ...

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)
pemindahan cairan
sumber
14

Ini seharusnya membuat Anda pergi. Nilai raster dibaca menggunakan rasterio , dan koordinat pusat piksel dikonversi ke Eastings / Northings menggunakan affine , yang kemudian dikonversi ke Latitude / Longitude menggunakan pyproj . Kebanyakan array memiliki bentuk yang sama dengan input 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
sumber
1
Saat menggunakan ini, saya mendapatkan pesan "AttributeError: objek 'DatasetReader' tidak memiliki atribut 'affine'" untuk baris "T0 = r.affine"
mitchus
@mitchus Rupanya affinehanya alias untuk transform, dan alias telah dihapus dari versi rasterio terbaru. Saya mengedit jawabannya tetapi sepertinya harus ditinjau oleh rekan sejak saya baru di sini. :)
Autumnsault
1
Itu juga terlihat seperti indeks yang salah A.shape, yang hanya memiliki dua dimensi.
Autumnsault