Extraindo valores rasterizados em pontos usando GIS de código aberto?

Como posso extrair valores de uma varredura por pontos?

Eu prefiro não em Arcgis.

Eu prefiro no Qgis ou Mapwindow ou outros gis de código aberto.

O QGIS "Point Sampling Tool" deve ser o plugin que você está procurando.

Aqui está uma descrição detalhada de como usá-lo: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Atualização baseada no comentário de Paolo:

o plug-in não é a única solução e nem sempre é a solução mais fácil. Uma solução alternativa é a função Saga 'Adicionar valores de varredura ao ponto' na caixa de ferramentas de processamento. Veja para detalhes http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

No PostGIS 2.0, você pode:

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

Certifique-se de que a sua varredura esteja lado a lado muito pequena ao carregá-la (-t 10x10 com a carregadora).

Eu estava tendo problemas com as ferramentas QGIS e SAGA GUI mencionadas neste tópico ( Raster values to pointsestava falhando por algum motivo e gerando erros inúteis e o GRASS v.samplecriou uma nova camada que não foi útil). Depois de falhar com as ferramentas da GUI por um tempo, tentei fazer isso na calculadora de campo. Funcionou muito bem e eu pude controlar o processo um pouco melhor do que as GUIs permitem e fazer alguns outros cálculos ao longo do caminho.

Digamos que você tenha uma camada com nome ptse outra com nome rast, ambos no mesmo sistema de coordenadas. Você gostaria de experimentar rastcada par X, Y representado em pts.

Se você nunca usou a Calculadora de Campo antes, é bem simples. Você digitará seu cálculo na caixa "Expressão" e Q fornecerá várias variáveis ​​e operações na coluna do meio, com o texto de ajuda na coluna da direita. Vou dividir esse processo em quatro etapas:

1. Abra a tabela de atributos da ptscamada com a qual você deseja experimentar.

2. Quando estiver na caixa de diálogo Calculadora de campo, escolha se você deseja criar um novo campo ou modificar um campo existente em sua ptscamada.

3. Em seguida, crie uma expressão para preencher a ptscoluna de atributo novo ou existente . Você pode começar modificando o código de expressão que funcionou para mim:

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

4. Você deve fornecer raster_value()um nome de camada raster 'rast', um número de banda 1e a geometria do ponto em make_point(). $xe $ysão variáveis ​​de geometria dependentes da localização do ponto em cada linha da tabela de atributos.

Este método também permite que as operações aritméticas como subtraindo o valor de uma outra camada raster chamada other_rasta partir rast, que me salvou um monte de tempo sobre as ferramentas GUI. Exemplo abaixo:

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

Observe novamente que as três camadas pts, raste other_rastdeve estar no mesmo sistema de coordenadas para que este método de trabalho.

Tente usar o QGIS 3.2.2 e o SAGA (instalado por padrão no QGIS): A função "Valores de varredura em pontos" fará tudo por você: Pega um arquivo de imagem e o converte em um formato vetorial de ponto, obtendo as informações da imagem rasterizada.

As ferramentas GME da Hawthorne Beyer fazem isso muito bem via linha de comando e permitem lotes fáceis com loops 'for'.

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

Referência de comando isectpntrst do GME

No GRASS GIS, você pode consultar o mapa na GUI ou usar http://grass.osgeo.org/gdp/html_grass64/r.what.html

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html possui um guia passo a passo para usar o pacote R Raster extrair valores raster de pontos.

Você pode usar isso: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

Está na caixa de ferramentas SAGA da Qgis! Faz tudo de uma só vez :)

Aqui está uma função que escrevi usando python e gdal. A função pega uma lista de rasters e um quadro de dados do pandas que contém as coordenadas do ponto e retorna um quadro de dados do pandas com as coordenadas do ponto, os centróides para as respectivas células raster e os respectivos valores das células. A função faz parte do pacote chorospy em desenvolvimento (encontrado aqui ).

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

Exemplo de como executá-lo, uma vez que os rasters estão no seu diretório de trabalho atual:

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

Se você tiver acesso ao FME , poderá usar um dos dois transformadores no FME Workbench.

O RasterCellCoercer ("Decompõe todos os recursos de varredura numérica de entrada em pontos ou polígonos individuais. Um recurso de vetor é gerado para cada célula na varredura.")

O PointOnRasterValueExtractor ("Considera recursos de ponto e uma única varredura de referência. A saída consiste no (s) valor (es) de banda e paleta no local de cada ponto.")

Pensamento rápido:

1. gdal_polygonize.py - poligoniza seu recurso de varredura
2. Inserir seus recursos de ponto e polígonos no banco de dados PostGIS
3. Use a função st_intersects para obter todos os valores de elevação onde os recursos se cruzam