Mengekstraksi nilai raster di titik menggunakan Open Source GIS?

Bagaimana saya bisa mengekstrak nilai dari raster berdasarkan poin?

Saya lebih suka tidak di Arcgis.

Saya lebih suka di Qgis atau Mapwindow atau open source gis lainnya.

QGIS "Alat Pengambilan Sampel Point" harus menjadi plugin yang Anda cari.

Berikut adalah uraian terperinci tentang cara menggunakannya: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Pembaruan berdasarkan komentar Paolo:

plugin bukan satu-satunya solusi, dan tidak selalu solusi termudah lagi. Solusi alternatif adalah fungsi Saga 'Tambahkan nilai raster ke titik' di kotak alat pemrosesan. Lihat untuk detail http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

Di PostGIS 2.0 Anda dapat melakukan:

SELECT ST_Value(rast, geom) val
FROM yourrastertabe, yourpointtable
WHERE ST_Intersects(rast, geom)

Pastikan raster Anda berukuran sangat kecil saat Anda memuatnya (-t 10x10 dengan loader).

Saya mengalami masalah dengan alat QGIS dan SAGA GUI yang disebutkan di utas ini ( Raster values to pointsgagal karena beberapa alasan dan melempar kesalahan yang tidak membantu dan GRASS v.samplemembuat layer baru yang sama sekali tidak membantu). Setelah gagal dengan alat GUI untuk sementara waktu, saya mencoba melakukan ini di Field Calculator. Ini bekerja cukup baik dan saya bisa mengendalikan prosesnya sedikit lebih baik daripada yang diizinkan GUI, dan membuat beberapa perhitungan lain di sepanjang jalan.

Katakanlah Anda memiliki layer bernama ptsdan nama lain rast, keduanya dalam sistem koordinat yang sama. Anda ingin mencicipi rastpada setiap X, pasangan Y yang diwakili dalam pts.

Jika Anda belum pernah menggunakan Field Calculator sebelumnya, ini cukup sederhana. Anda akan memasukkan perhitungan Anda di kotak "Ekspresi", dan Q memberi Anda sejumlah variabel dan operasi di kolom tengah, dengan teks bantuan menyala di kolom kanan. Saya akan memecah proses ini menjadi empat langkah:

1. Buka tabel atribut dari ptslapisan yang ingin Anda sampel.

2. Setelah Anda berada di dialog Kalkulator Bidang, pilih apakah Anda ingin Membuat bidang baru atau Memodifikasi bidang yang ada di ptslapisan Anda .

3. Selanjutnya, buat ekspresi untuk mengisi ptskolom atribut baru atau yang sudah ada . Anda mungkin mulai dengan mengubah kode ekspresi yang berfungsi untuk saya:

raster_value('rast', 1, make_point($x, $y))

4. Anda harus memberi raster_value()nama layer raster 'rast', nomor band 1, dan geometri titik di make_point(). $xdan $yapakah variabel geometri tergantung pada lokasi titik di setiap baris tabel atribut.

Metode ini juga memungkinkan operasi aritmatika seperti mengurangkan nilai lapisan raster lain yang disebut other_rastdari rast, yang menyelamatkan saya banyak waktu selama alat GUI. Contoh di bawah ini:

raster_value('rast', 1, make_point($x, $y)) - raster_value('other_rast', 1, make_point($x, $y))

Perhatikan lagi bahwa ketiga layer pts,, rastdan other_rastharus berada dalam sistem koordinat yang sama agar metode ini dapat berfungsi.

Coba gunakan QGIS 3.2.2 dan SAGA (terinstal secara default di QGIS): Fungsi "Nilai Raster ke Poin" akan melakukan segalanya untuk Anda: Dibutuhkan file gambar dan mengubahnya menjadi bentuk Point-vektor yang mengambil informasi dari gambar raster.

Alat GME Hawthorne Beyer melakukan ini dengan baik melalui baris perintah, dan memungkinkan batching mudah dengan loop 'for'.

isectpntrst(in="path/to/shapefile", raster="path/to/raster", field="fieldname")

Referensi perintah GME isectpntrst

Di GRASS GIS, Anda dapat menanyakan peta di GUI atau menggunakan http://grass.osgeo.org/gdp/html_grass64/r.what.html

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html memiliki panduan langkah demi langkah untuk menggunakan paket R Raster, ekstrak nilai raster dari poin.

Anda dapat menggunakan ini: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

Ada dalam SAGA Toolbox dari Qgis! Itu melakukan segalanya dalam satu langkah :)

Inilah fungsi yang saya tulis menggunakan python dan gdal. Fungsi mengambil daftar raster dan bingkai data panda yang berisi titik koordinat dan mengembalikan bingkai data panda dengan titik koordinat, centroid untuk masing-masing sel raster dan nilai sel masing-masing. Fungsi ini merupakan bagian dari paket chorospy paket yang sedang dikembangkan (ditemukan di sini ).

import pandas
import numpy
from osgeo import gdal

def getValuesAtPoint(indir, rasterfileList, pos, lon, lat):
    #gt(2) and gt(4) coefficients are zero, and the gt(1) is pixel width, and gt(5) is pixel height.
    #The (gt(0),gt(3)) position is the top left corner of the top left pixel of the raster.
    for i, rs in enumerate(rasterfileList):

        presValues = []
        gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
        gt = gdata.GetGeoTransform()
        band = gdata.GetRasterBand(1)
        nodata = band.GetNoDataValue()

        x0, y0 , w , h = gt[0], gt[3], gt[1], gt[5]

        data = band.ReadAsArray().astype(numpy.float)
        #free memory
        del gdata

        if i == 0:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                Xc = x0 + x*w + w/2 #the cell center x
                y = int((p[1][lat] - y0)/h)
                Yc = y0 + y*h + h/2 #the cell center y
                try:
                    if data[y,x] != nodata:
                        presVAL = [p[1][lon],p[1][lat], '{:.6f}'.format(Xc), '{:.6f}'.format(Yc), data[y,x]]
                        presValues.append(presVAL)
                except:
                    pass
            df = pandas.DataFrame(presValues, columns=['x', 'y', 'Xc', 'Yc', rs])
        else:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                y = int((p[1][lat] - y0)/h)
                try:
                    if data[y,x] != nodata:
                        presValues.append(data[y,x])
                except:
                    pass
            df[rs] = pandas.Series(presValues)
    del data, band
    return df

Contoh cara menjalankannya mengingat raster ada di direktori kerja Anda saat ini:

rasDf = getValuesAtPoint('.', ['raster1', 'raster2'], inPoints, 'x', 'y')

Jika Anda memiliki akses ke FME maka Anda dapat menggunakan salah satu dari dua transformer di Workbench FME.

The RasterCellCoercer ( "Terurai semua input numerik fitur raster ke titik-titik individu atau poligon. Salah satu fitur vektor adalah output untuk setiap sel dalam raster.")

The PointOnRasterValueExtractor ("Mengambil fitur titik dan raster referensi tunggal. Output terdiri dari nilai band dan palet di lokasi masing-masing titik.")

Pikiran cepat:

1. gdal_polygonize.py - poligon fitur raster Anda
2. Masukkan fitur poin Anda dan poligon ke dalam database PostGIS
3. Gunakan fungsi st_intersects untuk menarik semua nilai ketinggian tempat fitur saling berpotongan