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

21

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.

Vassilis
fonte
1
Portanto, você tem pontos e precisa extrair os valores da varredura sob esses pontos, ou precisa converter as células da varredura em pontos. Basta verificar antes de tentar descobrir a resposta.
Nathan W
O primeiro, eu tenho os pontos e preciso extrair os valores da varredura, nesses pontos. THNX !!
Vassilis

Respostas:

37

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/

underdark
fonte
5
As pessoas ainda estão encontrando o post acima mencionado através deste Q&A. No entanto, 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 grade ao ponto' na caixa de ferramentas de processamento. Veja para detalhes este post .
Ecodiv
Cuidado. Acabei de executar a Point Sampling Tool. 60.000 pontos e 13 rasters. Os resultados falharam no meu teste de 30 amostras aleatórias para cada ano. Essa ferramenta tem problemas com grandes conjuntos de dados. Eu não usaria.
Se você não sabe- apenas GIS
Apesar dos problemas mencionados com grandes conjuntos de dados, é muito útil para extrair todos os valores de várias bandas de uma só vez. Todas as outras soluções relacionadas ao QGIS suportam apenas a extração de uma banda (como GRASS r.what) ou proíbem o uso de varredura multibanda (como Saga - valores raster em pontos).
EikeMike
7

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).

Pierre Racine
fonte
7

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))
  1. 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.

Ian
fonte
1
esta é a melhor resposta para esta pergunta
BC B.
6

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.

Fernando MM
fonte
4

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

jbaums
fonte
2

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')
spyrostheodoridis
fonte
1

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.")

Mark Ireland
fonte
Não, eu não tenho ou uso o FME, é um aplicativo independente ou um plug-in?
Vassilis
0

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

fonte
abordagem interessante, porque ontem comecei a estudar como usar o Postgis.
Vassilis
Obrigado, é bastante simplista, mas funciona. Aqui está o que eu era capaz de produzir com esta abordagem: i.imgur.com/h8CGJ.png