Profil d'élévation 10 km de chaque côté d'une ligne

15

Comment puis-je obtenir un profil d'altitude pour une bande de terrain?

L'élévation la plus élevée dans un rayon de 10 km (de chaque côté de la ligne définie) doit être prise en compte.

J'espère que ma question est claire. Merci beaucoup d'avance.

Kara
la source
La ligne définissant votre profil est-elle une simple ligne droite ou est-elle composée de plusieurs segments avec des coins?
Jake
La ligne se compose de plusieurs segments. Mais tous les segments sont en ligne droite. :) Je veux dire qu'il n'y a pas de courbes.
Kara
Juste ... comme on dit ... spitballing ... mais pourriez-vous tamponner la ligne avec un tampon de 10 km. puis sélectionnez toutes les fonctionnalités dans le tampon ... puis sélectionnez les valeurs les plus élevées?
Ger
1
Pourriez-vous fournir une image de ce que vous voulez réaliser?
Alexandre Neto
@ Alex: le résultat que je veux est un graphique d'élévation régulier. Mais avec un tampon de 10 km afin que la valeur la plus élevée de 10 km de chaque côté du chemin sélectionné soit affichée sur le graphique.
Kara

Réponses:

14

À la suite des commentaires, voici une version qui fonctionne avec des segments de ligne perpendiculaires. Veuillez utiliser avec prudence car je ne l'ai pas testé à fond!

Cette méthode est beaucoup plus maladroite que la réponse de @ whuber - en partie parce que je ne suis pas un très bon programmeur et en partie parce que le traitement vectoriel est un peu un problème. J'espère que cela vous permettra au moins de commencer si vous avez besoin de segments de ligne perpendiculaires.

Vous aurez besoin d'avoir les packages Shapely , Fiona et Numpy Python installés (ainsi que leurs dépendances) pour exécuter ceci.

#-------------------------------------------------------------------------------
# Name:        perp_lines.py
# Purpose:     Generates multiple profile lines perpendicular to an input line
#
# Author:      JamesS
#
# Created:     13/02/2013
#-------------------------------------------------------------------------------
""" Takes a shapefile containing a single line as input. Generates lines
    perpendicular to the original with the specified length and spacing and
    writes them to a new shapefile.

    The data should be in a projected co-ordinate system.
"""

import numpy as np
from fiona import collection
from shapely.geometry import LineString, MultiLineString

# ##############################################################################
# User input

# Input shapefile. Must be a single, simple line, in projected co-ordinates
in_shp = r'D:\Perp_Lines\Centre_Line.shp'

# The shapefile to which the perpendicular lines will be written
out_shp = r'D:\Perp_Lines\Output.shp'

# Profile spacing. The distance at which to space the perpendicular profiles
# In the same units as the original shapefile (e.g. metres)
spc = 100

# Length of cross-sections to calculate either side of central line
# i.e. the total length will be twice the value entered here.
# In the same co-ordinates as the original shapefile
sect_len = 1000
# ##############################################################################

# Open the shapefile and get the data
source = collection(in_shp, "r")
data = source.next()['geometry']
line = LineString(data['coordinates'])

# Define a schema for the output features. Add a new field called 'Dist'
# to uniquely identify each profile
schema = source.schema.copy()
schema['properties']['Dist'] = 'float'

# Open a new sink for the output features, using the same format driver
# and coordinate reference system as the source.
sink = collection(out_shp, "w", driver=source.driver, schema=schema,
                  crs=source.crs)

# Calculate the number of profiles to generate
n_prof = int(line.length/spc)

# Start iterating along the line
for prof in range(1, n_prof+1):
    # Get the start, mid and end points for this segment
    seg_st = line.interpolate((prof-1)*spc)
    seg_mid = line.interpolate((prof-0.5)*spc)
    seg_end = line.interpolate(prof*spc)

    # Get a displacement vector for this segment
    vec = np.array([[seg_end.x - seg_st.x,], [seg_end.y - seg_st.y,]])

    # Rotate the vector 90 deg clockwise and 90 deg counter clockwise
    rot_anti = np.array([[0, -1], [1, 0]])
    rot_clock = np.array([[0, 1], [-1, 0]])
    vec_anti = np.dot(rot_anti, vec)
    vec_clock = np.dot(rot_clock, vec)

    # Normalise the perpendicular vectors
    len_anti = ((vec_anti**2).sum())**0.5
    vec_anti = vec_anti/len_anti
    len_clock = ((vec_clock**2).sum())**0.5
    vec_clock = vec_clock/len_clock

    # Scale them up to the profile length
    vec_anti = vec_anti*sect_len
    vec_clock = vec_clock*sect_len

    # Calculate displacements from midpoint
    prof_st = (seg_mid.x + float(vec_anti[0]), seg_mid.y + float(vec_anti[1]))
    prof_end = (seg_mid.x + float(vec_clock[0]), seg_mid.y + float(vec_clock[1]))

    # Write to output
    rec = {'geometry':{'type':'LineString', 'coordinates':(prof_st, prof_end)},
           'properties':{'Id':0, 'Dist':(prof-0.5)*spc}}
    sink.write(rec)

# Tidy up
source.close()
sink.close()

L'image ci-dessous montre un exemple de la sortie du script. Vous introduisez un fichier de formes représentant votre ligne médiane et spécifiez la longueur des lignes perpendiculaires et leur espacement. La sortie est un nouveau fichier de formes contenant les lignes rouges de cette image, chacune ayant un attribut associé spécifiant sa distance depuis le début du profil.

Exemple de sortie de script

Comme @whuber l'a dit dans les commentaires, une fois arrivé à ce stade, le reste est assez facile. L'image ci-dessous montre un autre exemple avec la sortie ajoutée à ArcMap.

entrez la description de l'image ici

Utilisez l' outil Entité vers raster pour convertir les lignes perpendiculaires en un raster catégoriel. Définissez le raster VALUEcomme Distchamp dans le fichier de formes en sortie. Rappelez - vous aussi de mettre l'outil de Environmentstelle sorte que Extent, Cell sizeet Snap rastersont les mêmes que pour votre DEM sous - jacente. Vous devriez vous retrouver avec une représentation raster de vos lignes, quelque chose comme ceci:

entrez la description de l'image ici

Enfin, convertissez ce raster en une grille entière (à l'aide de l' outil Int ou de la calculatrice raster) et utilisez-le comme zones d'entrée pour l' outil Statistiques zonales en tant que tableau . Vous devriez vous retrouver avec une table de sortie comme celle-ci:

entrez la description de l'image ici

Le VALUEchamp de ce tableau indique la distance depuis le début de la ligne de profil d'origine. Les autres colonnes donnent diverses statistiques (maximum, moyenne, etc.) pour les valeurs de chaque transect. Vous pouvez utiliser ce tableau pour tracer votre profil de résumé.

NB: Un problème évident avec cette méthode est que, si votre ligne d'origine est très ondulée, certaines des lignes de transect peuvent se chevaucher. Les outils de statistiques zonales dans ArcGIS ne peuvent pas traiter les zones qui se chevauchent, donc lorsque cela se produit, l'une de vos lignes de transect aura la priorité sur l'autre. Cela peut ou non être un problème pour ce que vous faites.

Bonne chance!

JamesS
la source
3
+1 C'est un bon début pour une grande contribution! Si vous regardez attentivement la deuxième figure, vous remarquerez des transects plus courts: ce sont ceux qui traversent près des virages. En effet, votre algorithme de calcul des transects suppose de manière incorrecte que le déplacement de chaque segment sera égal spc, mais les virages raccourcissent les déplacements. Au lieu de cela, vous devez normaliser le vecteur de direction transversale (diviser ses composants par la longueur du vecteur), puis multiplier cela par le rayon souhaité du transect.
whuber
Vous avez tout à fait raison - merci pour les commentaires @whuber! Avec un peu de chance maintenant ...
JamesS
Cher James, je vais essayer ça merci beaucoup. Cette solution convient parfaitement.
Kara
11

L'élévation la plus élevée à moins de 10 km est la valeur maximale du quartier calculée avec un rayon circulaire de 10 km, il suffit donc d'extraire un profil de cette grille maximale du quartier le long de la trajectoire.

Exemple

Voici un DEM ombragé avec une trajectoire (ligne noire allant de bas en haut):

DEM

Cette image mesure environ 17 kilomètres sur 10. J'ai choisi un rayon de seulement 1 km au lieu de 10 km pour illustrer la méthode. Son tampon de 1 km est illustré en jaune.

Le maximum de voisinage d'un DEM sera toujours un peu étrange, car il aura tendance à augmenter de valeur aux points où un maximum (un sommet de colline, peut-être) tombe juste au-delà de 10 km et un autre maximum à une altitude différente vient juste à moins de 10 km . En particulier, les sommets qui dominent leur environnement contribueront à des cercles de valeurs parfaits centrés au point d'élévation maximale locale:

Quartier max

Plus sombre est plus haut sur cette carte.

Voici un tracé des profils du DEM d'origine (bleu) et du maximum de voisinage (rouge):

Profils

Il a été calculé en divisant la trajectoire en points régulièrement espacés à 0,1 km de distance (en commençant à la pointe sud), en extrayant les élévations à ces points et en faisant un nuage de points joint des triplets résultants (distance depuis le début, l'élévation, l'élévation maximale). L'espacement des points de 0,1 km a été choisi pour être sensiblement plus petit que le rayon du tampon mais suffisamment grand pour que le calcul se fasse rapidement (il était instantané).

whuber
la source
Ce n'est pas tout à fait exact cependant, non? Au lieu d'un tampon circulaire autour de chaque point, une ligne orthogonale de 20 km ne devrait-elle pas être utilisée pour échantillonner la trame sous-jacente? C'est du moins ainsi que j'interpréterais l'exigence de Kara selon laquelle "la valeur la plus élevée dans un rayon de 10 km de chaque côté de la ligne" serait prise en compte.
Jake
4
@jake Je ne dirais pas "inexact": vous proposez simplement une interprétation alternative. "De chaque côté de" est un terme vague qui pourrait utiliser une meilleure qualification. Je peux proposer des solutions pour des interprétations comme la vôtre; une méthode utilise un maximum zonal. Cependant, son exécution est plus compliquée et beaucoup plus lente. Pourquoi ne voyons-nous pas d'abord ce que le PO pense de cette solution simple?
whuber
Mauvais choix de mots, je n'aurais pas dû utiliser "précis" - désolé
Jake
1
Parce que vous savez comment utiliser l'outil Profil, vous avez presque terminé. QGIS a des interfaces avec GRASS qui incluent des opérations de voisinage. Il suffit d'appliquer l'opération maximale de voisinage à l'aide de r.neighbors et de profiler son résultat.
whuber
1
@JamesS Vous ne voulez pas faire un décalage parallèle, vous voulez rendre chaque profil transversal perpendiculaire à la ligne centrale. (L'approche de décalage parallèle peut être implémentée exactement comme je le décris ici en utilisant un voisinage long et maigre approprié pour le calcul du voisinage max.) Je suis presque sûr que vous pouvez trouver du code sur ce site pour construire des ensembles de segments de ligne perpendiculaires également espacés le long d'une polyligne; c'est la partie difficile. Tout le reste consiste simplement à extraire les valeurs DEM le long de ces segments et à les résumer.
whuber
6

J'ai eu le même problème et j'ai essayé la solution de James S, mais je n'ai pas réussi à faire fonctionner le GDAL avec Fiona.

J'ai ensuite découvert l'algorithme SAGA "Cross Profiles" dans QGIS 2.4, et j'ai obtenu exactement le résultat que je voulais et que je suppose que vous recherchez aussi (voir ci-dessous).

entrez la description de l'image ici

Don Richardo
la source
Salut, je viens de tomber sur ce post il y a quelques années. Je suis confronté au même problème que le démarreur de fil, et, étant très nouveau dans (Q) GIS, je suis en quelque sorte heureux d'avoir atteint la photo ci-dessus. Comment puis-je utiliser les données? La couche Profils croisés montre l'élévation pour chaque point échantillonné, mais je voudrais demander de l'aide 1) trouver l'élévation maximale pour chaque ligne transversale 2) trouver les coordonnées de l'intersection avec le chemin d'origine 3) associer les élévations maximales à partir de 1 avec les coordonnées de 2. Quelqu'un peut-il aider, s'il vous plaît? Merci d'avance! Mal
Cpt Reynolds
6

Pour toute personne intéressée, voici une version modifiée du code JamesS créant des lignes perpendiculaires en utilisant uniquement les bibliothèques numpy et osgeo. Grâce à JamesS, sa réponse m'a beaucoup aidé aujourd'hui!

import osgeo
from osgeo import ogr
import numpy as np

# ##############################################################################
# User input

# Input shapefile. Must be a single, simple line, in projected co-ordinates
in_shp = r'S:\line_utm_new.shp'

# The shapefile to which the perpendicular lines will be written
out_shp = r'S:\line_utm_neu_perp.shp'

# Profile spacing. The distance at which to space the perpendicular profiles
# In the same units as the original shapefile (e.g. metres)
spc = 100

# Length of cross-sections to calculate either side of central line
# i.e. the total length will be twice the value entered here.
# In the same co-ordinates as the original shapefile
sect_len = 1000
# ##############################################################################

# Open the shapefile and get the data
driverShp = ogr.GetDriverByName('ESRI Shapefile')
sourceShp = driverShp.Open(in_shp, 0)
layerIn = sourceShp.GetLayer()
layerRef = layerIn.GetSpatialRef()

# Go to first (and only) feature
layerIn.ResetReading()
featureIn = layerIn.GetNextFeature()
geomIn = featureIn.GetGeometryRef()

# Define a shp for the output features. Add a new field called 'M100' where the z-value 
# of the line is stored to uniquely identify each profile
outShp = driverShp.CreateDataSource(out_shp)
layerOut = outShp.CreateLayer('line_utm_neu_perp', layerRef, osgeo.ogr.wkbLineString)
layerDefn = layerOut.GetLayerDefn() # gets parameters of the current shapefile
layerOut.CreateField(ogr.FieldDefn('M100', ogr.OFTReal))

# Calculate the number of profiles/perpendicular lines to generate
n_prof = int(geomIn.Length()/spc)

# Define rotation vectors
rot_anti = np.array([[0, -1], [1, 0]])
rot_clock = np.array([[0, 1], [-1, 0]])

# Start iterating along the line
for prof in range(1, n_prof):
    # Get the start, mid and end points for this segment
    seg_st = geomIn.GetPoint(prof-1) # (x, y, z)
    seg_mid = geomIn.GetPoint(prof)
    seg_end = geomIn.GetPoint(prof+1)

    # Get a displacement vector for this segment
    vec = np.array([[seg_end[0] - seg_st[0],], [seg_end[1] - seg_st[1],]])    

    # Rotate the vector 90 deg clockwise and 90 deg counter clockwise
    vec_anti = np.dot(rot_anti, vec)
    vec_clock = np.dot(rot_clock, vec)

    # Normalise the perpendicular vectors
    len_anti = ((vec_anti**2).sum())**0.5
    vec_anti = vec_anti/len_anti
    len_clock = ((vec_clock**2).sum())**0.5
    vec_clock = vec_clock/len_clock

    # Scale them up to the profile length
    vec_anti = vec_anti*sect_len
    vec_clock = vec_clock*sect_len

    # Calculate displacements from midpoint
    prof_st = (seg_mid[0] + float(vec_anti[0]), seg_mid[1] + float(vec_anti[1]))
    prof_end = (seg_mid[0] + float(vec_clock[0]), seg_mid[1] + float(vec_clock[1]))

    # Write to output
    geomLine = ogr.Geometry(ogr.wkbLineString)
    geomLine.AddPoint(prof_st[0],prof_st[1])
    geomLine.AddPoint(prof_end[0],prof_end[1])
    featureLine = ogr.Feature(layerDefn)
    featureLine.SetGeometry(geomLine)
    featureLine.SetFID(prof)
    featureLine.SetField('M100',round(seg_mid[2],1))
    layerOut.CreateFeature(featureLine)

# Tidy up
outShp.Destroy()
sourceShp.Destroy()
ket
la source
Merci ket - j'ai essayé cela mais cela ne fonctionne pas complètement pour moi malheureusement. J'ai donné au script un fichier de formes avec une seule entité polyligne mais ma sortie est juste la table d'attributs avec beaucoup de zéros pour la valeur "M100" - aucune entité n'apparaissant sur la carte. Des idées?
davehughes87
Peu importe - réalisé maintenant que votre script semble calculer des lignes perpendiculaires aux extrémités de chaque segment de la polyligne, pas tous les mètres "spc". Cela signifie que je manquais de polyligne pour travailler avec avant que n_prof ne soit atteint dans la boucle et que des valeurs "nan" soient produites.
davehughes87