使用开源GIS提取点上的栅格值？

#1 楼

QGIS“点采样工具”应该是您要寻找的插件。

以下是如何使用它的详细说明：http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

根据Paolo的评论进行更新：

插件不是唯一的解决方案，也不总是最简单的解决方案。另一种解决方案是处理工具箱中的Saga函数“将栅格值添加到点”。有关详细信息，请参见http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

#2 楼

我在此线程中提到的QGIS和SAGA GUI工具遇到了问题（Raster values to points由于某种原因失败并抛出了无用的错误，而GRASS v.sample创建了一个全新的无用层）。在使用GUI工具一段时间后，我尝试在Field Calculator中执行此操作。它工作得很好，并且我能够控制该过程比GUI允许的要好一些，并在此过程中进行了其他一些计算。在同一坐标系中。您想对pts中表示的每个X，Y对采样rast。

如果您以前没有使用过字段计算器，这很简单。您将在“表达式”（Expression）框中输入计算，Q在中间列为您提供许多变量和运算，在右列中显示帮助文本。我将这个过程分为四个步骤：


打开要采样的rast图层的属性表。
一旦进入“字段计算器”对话框，选择是要创建新字段还是修改pts图层中的现有字段。
接下来，构建一个表达式以填充新的或现有的pts属性列。您可能会先修改对我有用的表达式代码：

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

您必须为pts提供一个栅格图层名称pts，波段编号raster_value()，点几何在'rast'。 1和make_point()是几何变量，取决于点在属性表的每一行中的位置。我花了很多时间在GUI工具上。以下示例：

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

再次注意，三层$x，$y和other_rast必须在同一坐标系中才能起作用。

#3 楼

在PostGIS 2.0中，您可以执行以下操作：

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

确保在加载栅格时将栅格平铺（使用-t 10x10加载）。

#4 楼

尝试使用QGIS 3.2.2和SAGA（默认情况下已安装在QGIS中）：“栅格值转换为点”功能将为您完成所有工作：它将图像文件转换为点矢量形状，从而从光栅图像中获取信息。

#5 楼

Hawthorne Beyer的GME工具可以通过命令行很好地做到这一点，并允许使用“ for”循环轻松进行批处理。

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

GME isectpntrst命令参考

#6 楼

在GRASS GIS中，您可以在GUI中查询地图，也可以使用http://grass.osgeo.org/gdp/html_grass64/r.what.html

#7 楼

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html
具有逐步指南，可使用R Raster包从点提取栅格值。

#8 楼

您可以使用此工具：
http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

它在Qgis的SAGA工具箱中！它一步就能完成所有事情：)

#9 楼

这是我使用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

#10 楼

如果可以访问FME，则可以在FME Workbench中使用两个转换器之一。

RasterCellCoercer（“将所有输入的数字栅格要素分解为单个点或多边形。为每个像元输出一个矢量要素在栅格中。”）

PointOnRasterValueExtractor（“获取点要素和单个参考栅格。输出由每个点位置处的波段和调色板值组成。”）

#11 楼

快速思考：


gdal_polygonize.py-将栅格要素多边形化
将点要素和面插入PostGIS数据库中
使用st_intersects函数将所有高程值提取到要素相交