Extraire des valeurs raster à des points à l'aide du SIG Open Source?

Comment extraire des valeurs d'un raster par points?

Je préfère ne pas à Arcgis.

Je préfère dans Qgis ou Mapwindow ou d'autres SIG open source.

QGIS "Point Sampling Tool" devrait être le plugin que vous recherchez.

Voici une description détaillée de son utilisation: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Mise à jour basée sur le commentaire de Paolo:

le plugin n'est plus la seule solution, et n'est plus toujours la solution la plus simple. Une solution alternative est la fonction Saga «Ajouter des valeurs raster au point» dans la boîte à outils de traitement. Voir pour plus de détails http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

Dans PostGIS 2.0, vous pouvez faire:

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

Assurez-vous que votre raster est très petit en mosaïque lorsque vous le chargez (-t 10x10 avec le chargeur).

J'avais des problèmes avec les outils QGIS et SAGA GUI mentionnés dans ce fil ( Raster values to pointséchouait pour une raison quelconque et lançait des erreurs inutiles et GRASS a v.samplecréé une toute nouvelle couche qui n'était pas utile). Après avoir échoué avec les outils de l'interface graphique pendant un certain temps, j'ai essayé de le faire dans la calculatrice de champ. Cela a très bien fonctionné et j'ai pu contrôler le processus un peu mieux que les interfaces graphiques ne le permettent, et faire d'autres calculs en cours de route.

Supposons que vous ayez un calque nommé ptset un autre nommé rast, tous deux dans le même système de coordonnées. Vous souhaitez échantillonner rastà chaque paire X, Y représentée dans pts.

Si vous n'avez jamais utilisé la calculatrice de champ auparavant, c'est assez simple. Vous entrerez votre calcul dans la case "Expression", et Q vous donnera un certain nombre de variables et d'opérations dans la colonne du milieu, avec le texte d'aide dans la colonne de droite. Je vais diviser ce processus en quatre étapes:

1. Ouvrez la table des attributs de la ptscouche avec laquelle vous souhaitez échantillonner.

2. Une fois dans la boîte de dialogue Calculatrice de champ, choisissez si vous souhaitez créer un nouveau champ ou modifier un champ existant dans votre ptscouche.

3. Ensuite, créez une expression pour remplir le nouveau ou l'existant pts colonne d'attribut . Vous pourriez commencer par modifier le code d'expression qui a fonctionné pour moi:

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

4. Vous devez fournir raster_value()un nom de couche raster'rast' , un numéro de bande 1et la géométrie du point à make_point(). $xet $ysont des variables de géométrie qui dépendent de l'emplacement du point dans chaque ligne de la table attributaire.

Cette méthode permet également des opérations arithmétiques comme la soustraction de la valeur d'une autre couche raster appelée other_rast à partir rast, ce qui m'a sauvé un tas de temps sur les outils de l' interface graphique. Exemple ci-dessous:

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

Notez encore que les trois couches pts, rastet other_rastdoit être dans le même système de coordonnées pour que cette méthode fonctionne.

Essayez d'utiliser QGIS 3.2.2 et SAGA (installés par défaut dans QGIS): la fonction "Valeurs raster en points" fera tout pour vous: elle prend un fichier image et le convertit en une forme point-vecteur en prenant les informations de l'image raster.

Les outils GME de Hawthorne Beyer le font bien via la ligne de commande et permettent un batch facile avec des boucles 'for'.

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

Référence de la commande isectpntrst GME

Dans GRASS GIS, vous pouvez soit interroger la carte dans l'interface graphique ou utiliser http://grass.osgeo.org/gdp/html_grass64/r.what.html

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html a un guide étape par étape pour utiliser le package R Raster extraire des valeurs raster à partir de points.

Vous pouvez l'utiliser: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

C'est dans la boîte à outils SAGA de Qgis! Il fait tout en une seule étape :)

Voici une fonction que j'ai écrite en utilisant python et gdal. La fonction prend une liste de rasters et une trame de données pandas contenant les coordonnées du point et renvoie une trame de données pandas avec les coordonnées du point, les centroïdes pour les cellules raster respectives et les valeurs de cellule respectives. La fonction fait partie du paquet chorospy en cours de développement (trouvé ici ).

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

Exemple de la façon de l'exécuter étant donné que les rasters se trouvent dans votre répertoire de travail actuel:

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

Si vous avez accès à FME, vous pouvez utiliser l'un des deux transformateurs de FME Workbench.

Le RasterCellCoercer ( « Decomposes toutes les fonctions de trame numérique d'entrée en des points individuels ou des polygones. Une caractéristique de vecteur est sortie pour chaque cellule de la trame. »)

Le PointOnRasterValueExtractor ( « Prend en fonctions d'aiguillage et une seule trame de référence. La sortie se compose de la valeur de la bande et de la palette (s) à l'emplacement de chaque point ») .

Pensée rapide:

1. gdal_polygonize.py - polygone votre entité raster
2. Insérez vos entités ponctuelles et vos polygones dans la base de données PostGIS
3. Utilisez la fonction st_intersects pour extraire toutes les valeurs d'altitude où les entités se croisent