Perfil de elevación 10 km a cada lado de una línea

14

¿Cómo puedo obtener un perfil de elevación para una banda de terreno?

Se debe tener en cuenta la elevación más alta dentro de los 10 km (a cada lado de la línea definida).

Espero que mi pregunta sea clara. Muchas gracias de antemano.

    
pregunta Kara 06.02.2013 - 09:05

4 respuestas

14

A continuación de los comentarios, aquí hay una versión que funciona con segmentos de línea perpendiculares. Por favor, use con precaución ya que no lo he probado a fondo.

Este método es mucho más torpe que la respuesta de @ whuber, en parte porque no soy un muy buen programador, y en parte porque el procesamiento de vectores es un poco como un error. Espero que al menos lo ayude a comenzar si los segmentos de línea perpendiculares son lo que necesita.

Necesitará tener Shapely , Fiona y Numpy paquetes de Python instalados (junto con sus dependencias) para ejecutar esto.

#-------------------------------------------------------------------------------
# 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()

La imagen de abajo muestra un ejemplo de la salida del script. Se alimenta en un shapefile que representa su línea central, y especifica la longitud de las líneas perpendiculares y su espaciado. La salida es un nuevo shapefile que contiene las líneas rojas en esta imagen, cada una de las cuales tiene un atributo asociado que especifica su distancia desde el inicio del perfil.

[email protected],unavezquellegasaestaetapa,elrestoesbastantefácil.LaimagenacontinuaciónmuestraotroejemploconlasalidaagregadaaArcMap.

Use la herramienta Feature to Raster para convertir las líneas perpendiculares en una trama categórica. Establezca el ráster VALUE para que sea el campo Dist en el archivo de forma de salida. También recuerde configurar la herramienta Environments para que Extent , Cell size y Snap raster sean los mismos que para su DEM subyacente. Deberías terminar con una representación rasterizada de tus líneas, algo como esto:

Finalmente,conviertaesterásteraunacuadrículadeenteros(utilizandolaherramienta Int o la calculadora ráster), y úsela como zonas de entrada para el Estadísticas zonales como tabla herramienta. Deberías terminar con una tabla de salida como esta:

El campo VALUE en esta tabla indica la distancia desde el inicio de la línea de perfil original. Las otras columnas dan varias estadísticas (máximo, media, etc.) para los valores en cada transecto. Puede usar esta tabla para trazar su perfil de resumen.

NB: Un problema obvio con este método es que, si su línea original es muy ondulada, algunas de las líneas del transecto pueden superponerse. Las herramientas de estadísticas zonales en ArcGIS no pueden lidiar con las zonas superpuestas, por lo que cuando esto suceda, una de sus líneas de transecto tendrá prioridad sobre la otra. Esto puede o no ser un problema para lo que estás haciendo.

¡Buena suerte!

    
respondido por el JamesS 13.02.2013 - 23:52
11

La elevación más alta dentro de 10 km es el valor máximo de vecindario calculado con un radio circular de 10 km, así que simplemente extraiga un perfil de esta cuadrícula de vecindario máximo a lo largo de la trayectoria.

Ejemplo

Aquí hay un DEM sombreado con una trayectoria (línea negra que va de abajo hacia arriba):

Estaimagenesdeaproximadamente17por10kilómetros.Elegíunradiodesolo1kmenlugarde10kmparailustrarelmétodo.Subúferde1kmsemuestradelineadoenamarillo.

ElmáximodevecindariodeunDEMsiempreseveráunpocoextraño,porquetenderáaaumentarsuvalorenpuntosdondeunmáximo(unacimadeunacolina,quizás)caigamásde10kmyotromáximoenunaelevacióndiferentevienesolodentrodelos10km.Enparticular,lascolinasquedominansusalrededorescontribuiránconcírculosdevaloresperfectoscentradosenelpuntodeelevaciónmáximalocal:

Más oscuro es más alto en este mapa.

Aquí hay un gráfico de los perfiles del DEM original (azul) y el máximo de vecindario (Rojo):

Se calculó dividiendo la trayectoria en puntos espaciados regularmente a una distancia de 0,1 km (comenzando en el extremo sur), extrayendo las elevaciones en esos puntos y realizando un diagrama de dispersión de los triples resultantes (distancia desde el inicio, la elevación, el máximo elevación). El espaciado de puntos de 0,1 km se eligió para ser sustancialmente más pequeño que el radio del búfer, pero lo suficientemente grande para que el cálculo fuera más rápido (fue instantáneo).

    
respondido por el whuber 06.02.2013 - 20:44
6

Tuve el mismo problema e intenté la solución de James S, pero no conseguí que GDAL trabajara con Fiona.

Luego descubrí el algoritmo SAGA "Perfiles cruzados" en QGIS 2.4, obtuve exactamente el resultado que quería y supongo que también lo está buscando (ver más abajo).

    
respondido por el Don Richardo 02.07.2014 - 11:36
6

Para cualquiera que esté interesado, aquí hay una versión modificada del código JamesS que crea líneas perpendiculares usando solo las bibliotecas numpy y osgeo. Gracias a JamesS, ¡su respuesta me ayudó mucho hoy!

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()
    
respondido por el ket 28.11.2014 - 13:42

Lea otras preguntas en las etiquetas