¿Medir la distancia entre dos puntos considerando el raster de elevación?

Question

¿Es posible calcular la distancia entre 2 coordenadas (x1,y1,z1) - (x2,y2,z2) considerando un raster de elevación (celdas que contienen elevación)? En la imagen, hay 2 coordenadas (cruces rojas). Me gustaría obtener la distancia en verde. Sé cómo obtener la distancia directa teniendo en cuenta la elevación (línea discontinua roja en la imagen), pero esto no tiene en cuenta que puede haber una montaña o un valle entre los 2 puntos. Estoy usando R, pero si hay una buena manera de hacerlo en QGIS que sería genial también.

Answer 1

En R esto es bastante sencillo por dos puntos. Hay funciones para devolver la matriz de distancia y para extraer la elevación de los puntos basados en un DEM.

En primer lugar, añada los paquetes necesarios y, aquí, cree una trama ficticia a modo de ejemplo. Puede utilizar la función "raster" para leer su raster de elevación desde el disco.

library(raster)
library(sp)

e <- as(extent(-1942755,-1847241,2524840,2624901), "SpatialPolygons") 
  proj4string(e) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

r <- raster(e, res=1000)
  r[] <- runif(ncell(r),500,2500)

Ahora, utilizando raster::extract podemos extraer los valores de elevación de los dos puntos de su trayectoria y devolver la distancia en línea recta no corregida entre los puntos utilizando spDists .

xy <- list( coordinates( sampleRandom(r, 2, sp=TRUE) ) )
  elev <- extract(r, SpatialPoints(xy))
  d <- spDists( SpatialPoints(xy) )[2] 

A partir de los datos generados anteriormente, es sencillo derivar la distancia corregida.

elev.delta <- max(elev) - min(elev)
sqrt(d^2 + elev.delta^2)

La cosa se complica si se quiere corregir la distancia a lo largo de todo el trayecto, que en este caso vendría definida por el tamaño de la celda del raster subyacente a la línea recta o, incluso, a un trayecto irregular. Tendría que establecer una estructura de bucle para calcular las distancias corregidas para cada segmento de línea.

Un punto de partida es coaccionar sus puntos a un objeto de línea sp, calcular la longitud de la línea usando sp::SpatialLinesLengths (mismo resultado que sp::spDists ) y extraer todos los valores de elevación subyacentes con raster::extract .

En este punto, se configuraría para realizar un bucle a través de cada valor de elevación para devolver la distancia corregida para todos los deltas de elevación para luego añadirla a la distancia no corregida [d]. Esto sería efectivamente la distancia corregida de la pendiente.

l <- SpatialLines(list(Lines(list(Line(xy)),ID="1"))) 
( d <- SpatialLinesLengths(l) )
( elev <- unlist(extract(r, l)) )

Answer 2

Asumiendo que se puede acceder fácilmente a las coordenadas de los puntos y a sus valores de elevación, la distancia se puede calcular con un simple PyQGIS script. Para probar mi enfoque, tengo un plugin de QGIS que puede obtener coordenadas y valores de puntos para un raster. Para 2 puntos arbitrarios sobre mi raster (ver siguiente imagen):

El código python es el siguiente:

from math import sqrt

pt1 =  QgsPointV2(QgsWKBTypes.PointZ, 418119.497699, 4456713.69199, 1751)

pt2 = QgsPointV2(QgsWKBTypes.PointZ, 420665.705483, 4456708.34282, 2257)

distance = sqrt((pt1.x() - pt2.x())**2 + (pt1.y() - pt2.y())**2 + (pt1.z() - pt2.z())**2)

print "distance: {:.4f} m".format(distance)

Para mi ejemplo, después de ejecutar el código anterior en la consola de Python de QGIS obtuve:

distance: 2596.0044 m

Desgraciadamente, QgsPointV2 no tiene todavía una función de distancia como QgsPoint clase y QgsGeometry (también tiene una función de distancia) no puede ser instanciada con un QgsPointV2 objeto.

Answer 3

Posible solución: (tarda mucho en calcularse)

library(raster)
library(sp)

# create elevation raster
e <- as(extent(-1942755,-1847241,2524840,2624901), "SpatialPolygons") 
proj4string(e) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

r <- raster(e, res=1000)
r[] <- runif(ncell(r),500,2500)

# plot raster
x11()
plot(r)

# creat coordinates
xy <- data.frame(x= c(-1859255, -1876255), y=c(2586401, 2552401))

# add coordinates to plot
points(-1859255, 2586401, col="red", cex=3, pch= 19)
points(-1876255, 2552401, col="red", cex=3, pch= 19)

#creat spatial line
l <- SpatialLines(list(Lines(list(Line(xy)),ID="1"))) 

#add spatial line to plot
plot(l, add=TRUE, lwd=3, col="Blue")

#create points along line
library(rgeos)
library(sp)
numOfPoints  <-  gLength(l)*10 #segments of 10cm if l in meters
points<-spsample(l, n = numOfPoints, type = "regular")

#add points to plot
plot(points, add=TRUE)

#get elevation for each points
points.df<-as.data.frame(points)
points.df$elevation<-raster::extract(r, points.df)

#comput distance between each points

head(points.df)
output<-data.frame(length=1:length(points.df$x))
i.2<-(length(points.df$x)-1) #last point has no further point for a distance
for(i in 1:i.2){

x1<-points.df[i,"x"]
x2<-points.df[i+1,"x"]
y1<-points.df[i,"y"]
y2<-points.df[i+1,"y"]
z1<-points.df[i,"elevation"]
z2<-points.df[i+1,"elevation"]

Distance<-((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)^0.5

output[i,]<-Distance
}
output[length(points.df$x),]<-0 #last point has no further point for a distance

sum(output) # Length between A and B

Este método permite calcular la distancia entre dos coordenadas (A,B) considerando un raster de elevación. La "línea" se fragmenta en segmentos delimitados por subcoordenadas (puntos) trazados a igual distancia en la línea. Así, la línea une los puntos A y B siguiendo la superficie de la capa ráster. Cuando dos subcoordenadas están en la misma celda, la distancia es la horizontal (azul). Cuando las dos coordenadas están en celdas diferentes, la distancia es la longitud de la hipotenusa (rojo).

Answer 4

Su pregunta puede dividirse en dos partes:

• ¿Cómo crear una polilínea 3D sobre un dem? En este paso, tendrías que dibujar una polilínea 2D para la geometría 2D, luego convertirla en una línea 3D con cada vértice en un píxel de tu dem (o al menos lo suficientemente denso como para obtener una buena aproximación de la longitud) y cada valor Z es igual al dem debajo del vértice.

• ¿Cómo calcular la longitud de una polilínea 3d? En este paso, calcularás la distancia de cada segmento individual, teniendo en cuenta el valor del paso que asignes entre vértices + la diferencia de Z. Es una operación de álgebra geométrica que podrías hacer en Excel, supongo.

@xunilk ha proporcionado esta solución con GRASS (al que deberías poder acceder a través de QGIS) para una pregunta similar : ¿Cómo calcular la longitud 3D de un segmento de línea en QGIS?