Wydobywanie wartości rastrowych w punktach za pomocą GIS Open Source?

Jak mogę wyodrębnić wartości z rastra według punktów?

Wolę nie w Arcgis.

Wolę w Qgis, Mapwindow lub innych open source GIS.

QGIS „Point Sampling Tool” powinno być wtyczką, której szukasz.

Oto szczegółowy opis tego, jak z niego korzystać: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Aktualizacja na podstawie komentarza Paolo:

wtyczka nie jest jedynym rozwiązaniem i nie zawsze jest najłatwiejszym rozwiązaniem. Alternatywnym rozwiązaniem jest funkcja Saga „Dodaj wartości rastrowe do punktu” w przyborniku przetwarzania. Zobacz szczegóły http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

W PostGIS 2.0 możesz:

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

Upewnij się, że twój raster jest bardzo mały, gdy go załadujesz (-t 10x10 z modułem ładującym).

Miałem problemy z narzędziami QGIS i SAGA GUI wspomnianymi w tym wątku ( Raster values to pointsz jakiegoś powodu zawodziły i v.samplegenerowały nieprzydatne błędy, a GRASS stworzył zupełnie nową warstwę, która nie była pomocna). Po pewnym czasie awarii narzędzi GUI próbowałem to zrobić w kalkulatorze polowym. Działa całkiem dobrze i byłem w stanie kontrolować proces nieco lepiej, niż pozwalają na to GUI, i po drodze dokonać innych obliczeń.

Załóżmy, że masz nazwaną warstwę, ptsa drugą nazwano rast, obie w tym samym układzie współrzędnych. Chciałbyś próbkować rastna każdej pary X, Y reprezentowanej w pts.

Jeśli nie korzystałeś wcześniej z kalkulatora pola, jest to całkiem proste. Wprowadzisz swoje obliczenia w polu „Wyrażenie”, a Q daje ci szereg zmiennych i operacji w środkowej kolumnie, z tekstem pomocy w prawej kolumnie. Podzielę ten proces na cztery kroki:

1. Otwórz tabelę atrybutów ptswarstwy, z którą chcesz próbkować.

2. Po przejściu do okna dialogowego Kalkulator pola wybierz, czy chcesz utworzyć nowe pole, czy zmodyfikować istniejące pole w ptswarstwie.

3. Następnie zbuduj wyrażenie, aby wypełnić nową lub istniejącą ptskolumnę atrybutu. Możesz zacząć od zmodyfikowania kodu wyrażenia, który dla mnie działał:

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

4. Musisz podać raster_value()nazwę warstwy rastrowej 'rast', numer pasma 1i geometrię punktu w make_point(). $xi $ysą zmiennymi geometrycznymi zależnymi od położenia punktu w każdym rzędzie tabeli atrybutów.

Metoda ta pozwala także operacje arytmetyczne jak odjęcie wartości innej warstwie rastrowej zwanej other_rastod rast, co zaoszczędziło mi sporo czasu nad narzędzi graficznych. Przykład poniżej:

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

Zauważ, że znowu z trzech warstw pts, rasti other_rastmusi być w tym samym układzie współrzędnych dla tej metody do pracy.

Spróbuj użyć QGIS 3.2.2 i SAGA (domyślnie instalowanych w QGIS): Funkcja „Wartości rastrowe do punktów” zrobi wszystko za Ciebie: pobiera plik obrazu i konwertuje go do kształtu wektora punktowego, pobierając informacje z obrazu rastrowego.

Narzędzia GME Hawthorne Beyera ładnie to robią za pomocą wiersza poleceń i pozwalają na łatwe grupowanie w pętle „for”.

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

Odwołanie do komendy GME isectpntrst

W GRASS GIS możesz wykonać zapytanie do mapy w GUI lub użyć http://grass.osgeo.org/gdp/html_grass64/r.what.html

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html ma przewodnik krok po kroku, aby użyć pakietu R Raster wyodrębnić wartości rastrowe z punktów.

Możesz użyć tego: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

Jest w SAGA Toolbox Qgis! Robi wszystko w jednym kroku :)

Oto funkcja, którą napisałem używając Pythona i Gdal. Funkcja przyjmuje listę rastrów i ramkę danych pandy zawierającą współrzędne punktu i zwraca ramkę danych pandy ze współrzędnymi punktu, centroidami dla odpowiednich komórek rastrowych i odpowiednich wartości komórek. Ta funkcja jest częścią pakietu rozwijanego pakietu chorospy (znajdującego się tutaj ).

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

Przykład tego, jak go uruchomić, biorąc pod uwagę, że rastry znajdują się w bieżącym katalogu roboczym:

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

Jeśli masz dostęp do FME , możesz użyć jednego z dwóch transformatorów w FME Workbench.

RasterCellCoercer ( „Rozkłada wejściowych numeryczny wszystkie cechy rastrowych do poszczególnych punktów lub wielokątów. Jedną z cech wektor jest wyjście dla każdej komórki rastra”).

PointOnRasterValueExtractor ( „przyjmuje w funkcji punktów rastrowych i jednym odniesienia. Wynik polega na wartość zakresu i palety (-ów) w miejscu, w każdym punkcie”).

Szybka myśl:

1. gdal_polygonize.py - poligonizuje funkcję rastrową
2. Wstaw swoje punkty i wielokąty do bazy danych PostGIS
3. Użyj funkcji st_intersects, aby pobrać wszystkie wartości rzędnych tam, gdzie przecinają się elementy