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

21

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.

Vassilis
la source
1
Vous avez donc des points et vous devez extraire les valeurs du raster sous ces points, ou devez-vous convertir les cellules raster en points. Je vérifie juste avant d'essayer de trouver la réponse.
Nathan W
Le premier, j'ai les points et j'ai besoin d'extraire les valeurs du raster, sous ces points. THNX !!
Vassilis

Réponses:

37

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/

obscur
la source
5
Les gens trouvent toujours le poste mentionné ci-dessus grâce à cette Q&R. Cependant, le plugin n'est pas la seule solution, et n'est plus toujours la solution la plus simple. Une solution alternative est la fonction Saga «Ajouter des valeurs de grille au point» dans la boîte à outils de traitement. Voir pour plus de détails cet article .
Ecodiv
Mise en garde. Je viens de lancer l'outil d'échantillonnage de points. 60 000 points et 13 rasters. Les résultats ont échoué mon test de 30 échantillons aléatoires pour chaque année. Cet outil a des problèmes avec les grands ensembles de données. Je ne l'utiliserais pas.
Si vous ne savez pas - juste GIS
Malgré les problèmes mentionnés avec les grands ensembles de données, il est très utile pour extraire toutes les valeurs multibandes en une seule fois. Toutes les autres solutions liées à QGIS ne prennent en charge que l'extraction d'une bande (comme GRASS r.what) ou interdisent l'utilisation de rasters multibandes (comme Saga - valeurs de raster en points).
EikeMike
7

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

Pierre Racine
la source
7

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

Ian
la source
1
c'est la meilleure réponse à cette question
BC B.
6

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.

Fernando MM
la source
4

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

jbaums
la source
2

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

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

Mark Ireland
la source
Non, je n'ai pas ou n'utilise pas FME, est-ce une application autonome ou un plugin?
Vassilis
0

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

la source
approche intéressante, car hier, commencez à étudier comment utiliser Postgis.
Vassilis
Merci, c'est assez simpliste mais ça marche. Voici ce que j'ai pu produire avec cette approche: i.imgur.com/h8CGJ.png