I'm working on a scientific programming project. I have a DEM and I'm supposed to calculate the slope all the long of trajectory. I had a problem with the interrpolation when I've plotted the graph of slope(distance) it shows discontinuity points.
Here's the code:
# -*- coding: utf-8 -*-
>
Created on Wed Jan 22 14:23:30 2020
@author: kingD
import math
import numpy as np
from PIL import Image
import scipy as sp
from time import sleep
import json
from Calcul_carto import *
import matplotlib.pyplot as plt
def symetrie_image(matrice):
"""Fonction qui fait lq symetrie par rapport a l'axe verticale
"""
d = np.shape(matrice)
M = np.zeros_like(matrice)
for i in range(d[0]):
for j in range(d[1]):
M[i,j] = matrice[i, d[1]-1-j]
return M
#Definition des besoins
path1 = "route.geojson"
path2 = "trace.geojson"
path3 = "La_Blanche.tif"
coord_origine = np.genfromtxt("La_Blanche.tfw")
x0, y0 = coord_origine[4], coord_origine[5]
altitudes = Geotiff_Array(path3)
altitudes = symetrie_image(altitudes)
altitudes = np.rot90(altitudes, 1)
dx, dy = 5, -5
pentes = Carte_Slope_Aspect(altitudes,dx,dy)[0]
#charger les donnees
with open (path1, "r") as read_file :
data1 = json.load( read_file )
with open (path2, "r") as read_file :
data2 = json.load( read_file )
chemin = data1["features"][0]["geometry"]["coordinates"]
tracer = data2["features"][0]["geometry"]["coordinates"]
def get_indices(x, y):
"""
Fonction qui recupere les indices correspondants au point de coordonnees (x, y)
:paramchemin: liste de coordonnées
:type chemin: liste(liste(float))
:return: indices des cellules correspondantes sur le MNT
:rtype: array((n, 2),int))
"""
global x0, y0
return (int(abs((x-x0)/5)), int(abs((y-y0)/5)))
def position_centre(i, j):
"""Fonction qui calcule les coordonnees du centre du pxel (i,j)
:param i,j: ligne, colonne
:type i,j: int
:return: coordonnees du centre
:rtype: tuple(int)
"""
global x0, y0
xc, yc = x0+5*i, y0+5*j
return (xc, yc)
def calculerDelta(a,b,c):
"""Fonction qui calcule delta
"""
return b*b-4*a*c
def resoudreEquationSecondDegre(a,b,c):
delta = calculerDelta(a,b,c)
racineDeDelta=math.sqrt(delta)
x1, x2 = (-b-racineDeDelta)/(2*a), (-b+racineDeDelta)/(2*a)
return (x1, x2)
def droite(A, B):
"""Fontio qui determine l'equation d'une droite
:param A, B: points
:type A, B: tuples
:return: caracteristiques de la droite(a,b) y = ax + b
:rtype: tuple
"""
XA ,XB, YA, YB = A[0],B[0], A[1], B[1]
a = (YB - YA)/(XB - XA)
b = YA - a*XA
return (a,b)
def coordonnees_pole(numero, A, B):
"""Fonction qui calcule les coordonnees d'un pole entre A et B
:param numero: numero du pole
:type numero: int
:param A, B: segment a interpoler
:type A, B: tuples
:return: coordonnes du pole portant le numero
:rtype: tuple
"""
a, b = droite(A, B)
xA, yA = A[0], A[1]
n, d = intercale_points(A, B)
x1, x2 = resoudreEquationSecondDegre((a**2)+1, (2*a*(b-yA)-2*xA), ((b-yA)**2)+(xA**2)-(numero*d)**2)
y1, y2 = a*x1 + b, a*x2 + b
if (droite(A,(x1, y1))[0])*a < 0: #si le signe de la pente entre A et le pole(x1,y1) est different de a
X, Y = x2, y2
else:
X, Y = x1, y1
return (X,Y)
def distance(A, B):
"""Fonction qui determine la distance entre A et B
:param A,B: points
:type A,B: floats
:return: distance entre A et B
:rtype: float
"""
xA, yA = A[0], A[1]
xB, yB = B[0], B[1]
return np.sqrt((xB-xA)**2 + (yB-yA)**2)
def intercale_points(A, B):
"""Fonction qui determine le nombre de points a intercaler sur le segment [AB]
:param A, B: coordonnees des points du segment
:type A, B: tuple
:return: nopmbre de points a intercaler etre A et B et la distance entre deux poles
:rtype: tuple(int, float)
"""
d = distance(A,B)
n = int(d/(50**0.50))
return (n, d/n)
def extract_decimal(x):
"""Fonction qui extrait la partie decimale d'un flottant
"""
return x - int(x)
def interpolation(emprise, P):
"""Fonction qui interpolele points P(x,y)
:param emprise: matrice qui represente les valeur des pixels
:type emprise: array
:pram P: point a determiner la valeur d'interpolation
:type P: tuple
:return: valeur d'interpolation
:rtype: float
"""
global x0, y0
xP, yP = P
ligne, colonne = get_indices(xP, yP)
xc, yc = position_centre(ligne, colonne)
if xP <= xc and yP >= yc:
i, j = ligne-1, colonne-1
elif xP >= xc and yP <= yc:
i, j = ligne, colonne
elif xP >= xc and yP >= yc:
i, j = ligne-1, colonne
elif xP <= xc and yP <= yc:
i, j = ligne, colonne-1
x, y = (xP-x0)/5, (yP-y0)/5 # coordonee de P dans le repere qui est la matrice
t, u = extract_decimal(x-i), extract_decimal(y-j)
VA = (1.0-u)*emprise[i][j] + u*emprise[i][j+1]
VB = (1.0-u)*emprise[i+1][j] + u*emprise[i+1][j+1]
val = (1.0-t)*VA + t*VB
return val
def segment(emprise,A , B, Dist):
"""Fonction qui retourne les differents poles entre A et B
:param emprise: matrice de valeurs
:type emprise: array
:param A,B: segment
:type: tuples
:param Dist: distance parcourue sur l'itineraire jusqu'a A
:type Dist: float
:return: tableau de valeur des poles et abscisse
:rtype: array
"""
xA, yA = A
xB, yB = B
iA, jA = get_indices(xA, yA)
iB, jB = get_indices(xB, yB)
T = []
n = intercale_points(A, B)[0]
for numero in range(n):
coor = coordonnees_pole(numero+1, A, B)
d = Dist + distance((coor[0], coor[1]),A)
z = interpolation(emprise, (coor[0], coor[1]))
T.append((d, z, (coor[0], coor[1])))
Dist += distance(A,B)
T.append((Dist, emprise[iB][jB], (xB, yB)))
return T, Dist
def pente_terrain(emprise, chem):
"""Fonction qui trace la rouite le long du ficheuir route.json
:param emprise: matrice de pixels
:type emprise: array
:param chem: liste de points correspondant de la route
:type chem: array
:return: liste de couple(abscisse, altitude)
:rtype: liste(tuple)
"""
Dep = chem[0]
Dist = 0
iDep, jDep = get_indices(Dep[0], Dep[1])
T = [(Dist, emprise[iDep][jDep],(Dep[0], Dep[1]))] # le point de depart de l'itineraire
for i in range(len(chem)-1):
A = (chem[i][0], chem[i][1])
B = (chem[i+1][0], chem[i+1][1])
S, Dist = segment(emprise, A , B, Dist)
T = T+S
return Transforme_en_vecteurs(T)
def pente_itineraire(emprise, chem):
"""Fonction qui trace la rouite le long du ficheuir route.json
:param emprise: matrice de pixels
:type emprise: array
:param chem: liste de points correspondant de la route
:type chem: array
:return: liste de couple(abscisse, altitude)
:rtype: liste(tuple)
"""
Dep = chem[0]
Dist = 0
iDep, jDep = get_indices(Dep[0], Dep[1])
T = [(Dist, emprise[iDep][jDep], (Dep[0], Dep[1]))] # le point de depart de l'itineraire
for i in range(len(chem)-1):
A = (chem[i][0], chem[i][1])
B = (chem[i+1][0], chem[i+1][1])
S, Dist = segment(emprise, A , B, Dist)
T = T+S
# calcul des pentes sur l'itineraire
D, Pente = [], []
for i in range(len(T)-1):
x, y = T[i][2]
z = T[i][1]
x1, y1 = T[i+1][2]
z1 = T[i+1][1]
p = np.arctan((z-z1)/distance((x, y),(x1, y1)))
Pente.append(p*180/np.pi)
D.append(T[i][0])
Pente.append(0) #on met la derniere pente a 0
D.append(T[i+1][0])
return D, Pente
def Transforme_en_vecteurs(T):
"""Fonction qui transforme une liste de tuple en patrice
:param T: liste de tuples
:type T: liste(tuple)
:return: matrice (2,n)
:rtype: array
"""
D, Z = [], []
for i in T:
D.append(i[0])
Z.append(i[1])
return (D, Z)
if __name__=="__main__":
# avec les altitudes
D, P1 = pente_itineraire(altitudes, chemin)
plt.plot(D, P1)
plt.xlabel('distances a l'' origine (m)')
#plt.ylabel('altitudes (m)')
#plt.title('altitudes sur l'' itineraire')
# avec les pentes
P2 = pente_terrain(pentes, chemin)[1]
plt.plot(D, P2)
plt.ylabel('pentes (m)')
plt.title('pentes sur l'' itineraire')
plt.show()
