0

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.