¿Extracción de valores ráster en puntos utilizando Open Source GIS?

¿Cómo puedo extraer valores de un ráster por puntos?

Prefiero no en Arcgis.

Prefiero en Qgis o Mapwindow u otros gis de código abierto.

QGIS "Point Sampling Tool" debería ser el complemento que estás buscando.

Aquí hay una descripción detallada de cómo usarlo: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Actualización basada en el comentario de Paolo:

el complemento no es la única solución, y ya no es siempre la solución más fácil. Una solución alternativa es la función Saga 'Agregar valores ráster al punto' en la caja de herramientas de procesamiento. Ver para más detalles http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

En PostGIS 2.0 puedes hacer:

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

Asegúrese de que su ráster esté en mosaico muy pequeño cuando lo cargue (-t 10x10 con el cargador).

Estaba teniendo problemas con las herramientas GUI de QGIS y SAGA mencionadas en este hilo ( Raster values to pointsestaba fallando por alguna razón y arrojando errores inútiles y GRASS v.samplecreó una capa completamente nueva que no fue útil). Después de fallar con las herramientas GUI por un tiempo, intenté hacer esto en la Calculadora de campo. Funcionó bastante bien y pude controlar el proceso un poco mejor de lo que permiten las GUI, y hacer algunos otros cálculos en el camino.

Digamos que tiene una capa con nombre ptsy otra con nombre rast, ambas en el mismo sistema de coordenadas. Te gustaría probar rasten cada par X, Y representado en pts.

Si no ha usado la Calculadora de campo antes, es bastante simple. Ingresará su cálculo en el cuadro "Expresión", y Q le dará una serie de variables y operaciones en la columna central, con el texto de ayuda en la columna derecha. Romperé este proceso en cuatro pasos:

1. Abra la tabla de atributos de la ptscapa con la que desea muestrear.

2. Una vez que esté en el cuadro de diálogo Calculadora de campo, elija si desea Crear un nuevo campo o Modificar un campo existente en su ptscapa.

3. A continuación, cree una expresión para llenar la ptscolumna de atributos nuevos o existentes . Puede comenzar modificando el código de expresión que funcionó para mí:

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

4. Debe proporcionar raster_value()un nombre de capa ráster 'rast', un número de banda 1y la geometría del punto en make_point(). $xy $yson variables de geometría que dependen de la ubicación del punto en cada fila de la tabla de atributos.

Este método también permite operaciones aritméticas como restando el valor de otra capa de trama llamada other_rastdesde rastque me salvó un montón de tiempo sobre las herramientas GUI. Ejemplo a continuación:

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

Observe nuevamente que las tres capas pts, rasty other_rastdeben estar en el mismo sistema de coordenadas para que este método funcione.

Intente usar QGIS 3.2.2 y SAGA (instalado de manera predeterminada en QGIS): la función "Valores ráster a puntos" hará todo por usted: toma un archivo de imagen y lo convierte en una forma de vector de punto tomando la información de la imagen ráster.

Las herramientas GME de Hawthorne Beyer lo hacen muy bien a través de la línea de comandos y permiten un procesamiento por lotes fácil con bucles 'for'.

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

Referencia de comando GME isectpntrst

En GRASS GIS, puede consultar el mapa en la GUI o usar http://grass.osgeo.org/gdp/html_grass64/r.what.html

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html tiene una guía paso a paso para usar el paquete R Raster para extraer valores ráster de puntos.

Puede usar esto: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

¡Está en la caja de herramientas SAGA de Qgis! Hace todo en un solo paso :)

Aquí hay una función que escribí usando python y gdal. La función toma una lista de rásteres y un marco de datos de pandas que contiene las coordenadas del punto y devuelve un marco de datos de pandas con las coordenadas del punto, los centroides para las celdas de ráster respectivas y los valores de celda respectivos. La función es parte del paquete chorospy del paquete en desarrollo (que se encuentra aquí ).

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

Ejemplo de cómo ejecutarlo dado que los rásteres están en su directorio de trabajo actual:

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

Si tiene acceso a FME , puede usar uno de los dos transformadores en FME Workbench.

El RasterCellCoercer ( "descompone todas las características de mapa de bits de entrada numérico en puntos individuales o polígonos. Un vector de características es de salida para cada celda en el ráster.")

El PointOnRasterValueExtractor ( "Toma en entidades de punto y una sola trama de referencia. La salida consiste en la banda y la paleta de valor (s) en la ubicación de cada punto.")

Pensamiento rápido:

1. gdal_polygonize.py - poligonaliza tu raster raster
2. Inserte sus entidades de puntos y polígonos en la base de datos PostGIS
3. Use la función st_intersects para extraer todos los valores de elevación donde las entidades se cruzan