Извлечение растровых значений в точках с использованием Open Source GIS?

21

Как извлечь значения из растра по точкам?

Я предпочитаю не в Arcgis.

Я предпочитаю в Qgis или Mapwindow или других открытых источниках ГИС.

Василис
источник
1
Таким образом, у вас есть точки, и вам нужно извлечь значения из растра под этими точками, или вам нужно преобразовать растровые ячейки в точки. Просто проверяю, прежде чем пытаться выработать ответ.
Натан W
Во-первых, у меня есть точки, и мне нужно извлечь значения из растра, под этими точками. Thnx !!
Василис

Ответы:

37

QGIS "Point Sampling Tool" должен быть тем плагином, который вы ищете.

Вот подробное описание того, как его использовать: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Обновление на основе комментария Паоло:

плагин - не единственное решение и не всегда самое простое решение. Альтернативным решением является функция Saga «Добавить растровые значения в точку» на панели инструментов обработки. Подробнее см. Http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

Подземье
источник
5
Люди все еще находят вышеупомянутый пост через этот Q & A. Однако плагин - не единственное решение и не всегда самое простое решение. Альтернативным решением является функция Saga «Добавить значения сетки в точку» на панели инструментов обработки. Подробности смотрите в этом посте .
Экодив
Внимание. Я только что запустил инструмент выбора точек. 60000 точек и 13 растров. Результаты не прошли мой тест на 30 случайных выборок за каждый год. Этот инструмент имеет проблемы с большими наборами данных. Я бы не использовал это.
Если вы не знаете - просто ГИС
Несмотря на упомянутые проблемы с большими наборами данных, это очень полезно для извлечения всех многополосных значений за один раз. Все другие решения, связанные с QGIS, поддерживают только извлечение одной полосы (например, GRASS и т.д.) или запрещают использование многоканального растра (например, Saga - значения растра в точках).
EikeMike
7

В PostGIS 2.0 вы можете сделать:

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

Убедитесь, что ваш растр очень маленький, когда вы загружаете его (-t 10x10 с загрузчиком).

Пьер Расин
источник
7

У меня были проблемы с инструментами графического интерфейса QGIS и SAGA, упомянутыми в этой теме ( Raster values to pointsпо какой-то причине произошел сбой и возникли бесполезные ошибки, а GRASS v.sampleсоздала совершенно новый слой, который не помог). После некоторого сбоя инструментов GUI я попытался сделать это в Калькуляторе поля. Это сработало довольно хорошо, и я смог контролировать процесс немного лучше, чем позволяют графические интерфейсы, и по пути делать некоторые другие вычисления.

Скажем, у вас есть слой с именем ptsи другим именем rast, оба в одной системе координат. Вы хотите сделать выборку rastдля каждой пары X, Y, представленной в pts.

Если вы раньше не пользовались полевым калькулятором, это довольно просто. Вы введете свой расчет в поле «Выражение», а Q предоставит вам ряд переменных и операций в среднем столбце, с текстом справки в правом столбце. Я разделю этот процесс на четыре этапа:

  1. Откройте таблицу атрибутов ptsслоя, с которым вы хотите выполнить выборку.

  2. В диалоговом окне «Калькулятор полей» выберите, хотите ли вы создать новое поле или изменить существующее поле в своем ptsслое.

  3. Затем создайте выражение, чтобы заполнить новый или существующий ptsстолбец атрибута. Вы могли бы начать с изменения кода выражения, который работал для меня:

raster_value('rast', 1, make_point($x, $y))
  1. Вы должны предоставить raster_value()имя растрового слоя 'rast', номер полосы 1и геометрию точки в make_point(). $xи $yявляются геометрическими переменными, зависящими от местоположения точки в каждой строке таблицы атрибутов.

Этот метод также позволяет выполнять арифметические операции , такие как вычитая значение другого растрового слоя под названием other_rastот rast, который спас мне кучу времени над инструментами графического интерфейса. Пример ниже:

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

Отметим еще раз , что в три слоя pts, rastи other_rastдолжны находиться в одной и той же системе координат , для этого метода к работе.

Ян
источник
1
это лучший ответ на этот вопрос
до н. э.
6

Попробуйте использовать QGIS 3.2.2 и SAGA (установленные по умолчанию в QGIS): функция «Растровые значения в точки» сделает все за вас: она берет файл изображения и преобразует его в форму вектора точек, получая информацию из растрового изображения.

Фернандо М.М.
источник
3

В GRASS GIS вы можете запросить карту в графическом интерфейсе или использовать http://grass.osgeo.org/gdp/html_grass64/r.what.html

markusN
источник
2
Для травы я нашел это: pvanb.wordpress.com/2010/05/05/…
Vassilis
Другим модулем GRASS, способным производить выборку растра, является v.sample, который также доступен через набор инструментов обработки QGIS.
user55937
2

Вот функция, которую я написал, используя python и gdal. Функция берет список растров и кадр данных pandas, содержащий координаты точки, и возвращает кадр данных pandas с координатами точки, центроидами для соответствующих растровых ячеек и соответствующими значениями ячеек. Эта функция является частью пакета chorospy, находящегося в стадии разработки (находится здесь ).

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

Пример того, как его запустить, учитывая, что растры находятся в вашем текущем рабочем каталоге:

rasDf = getValuesAtPoint('.', ['raster1', 'raster2'], inPoints, 'x', 'y')
spyrostheodoridis
источник
1

Если у вас есть доступ к FME, вы можете использовать один из двух трансформаторов в FME Workbench.

RasterCellCoercer ( «Разлагается все входные числовыми особенности растровых точек на отдельные или многоугольники. Один вектор признаков выводятся для каждой ячейки в растре.»)

PointOnRasterValueExtractor ( «принимает в точечных функций и одного опорного растра. Вывод состоит из полосы и палитры значения (ов) в месте расположения каждой точки.»)

Марк Ирландия
источник
Нет у меня нет или я использую FME, это отдельное приложение или плагин?
Василис
0

Быстрая мысль:

  1. gdal_polygonize.py - полигонизировать растровые объекты
  2. Вставьте свои точечные объекты и полигоны в базу данных PostGIS
  3. Используйте функцию st_intersects, чтобы получить все значения высот, где пересекаются объекты

источник
Интересный подход, потому что вчера начали изучать, как использовать Postgis.
Василис
Спасибо, это довольно упрощенно, но работает. Вот что я смог сделать с помощью этого подхода: i.imgur.com/h8CGJ.png