############################################################ # # # Virtual Laboratory of Geomathematics in Python # # # # Study of the topographic profile of a hill(26.09.2017) # # # # Complutense University of Madrid, Spain # # # # THIS SCRIPT IS PROVIDED BY THE AUTHORS "AS IS" AND # # CAN BE USED BY ANYONE FOR THE PURPOSES OF EDUCATION # # AND RESEARCH. # # # ############################################################ #!/usr/bin/env python import sys from pylab import * import xlrd import numpy as np from scipy.signal import savgol_filter # Ventana informacion print("In this script:\n 1) We selected a file with a height profile\n 2) The height profile is plotted\n 3) First derivative (slope) and second derivative are calculated\n 4) Both of them are plotted") # Eleccion de fichero fichero=raw_input("Please choose an excel file (including path)\n") print(fichero) # Funcion de apertura y lectura de la columna n de un fichero excel def read_excel(fichero,n): book = xlrd.open_workbook(fichero) sh = book.sheet_by_index(0) # Metemos la columna n en una lista coln=[] for rx in range(1,sh.nrows): coln.append(float(sh.cell_value(colx=n-1,rowx=rx))) return coln # Si no se ha seleccionado el fichero se cierra el programa if fichero == None: print("No file selected") sys.exit() # Lee las posiciones en x numx=int(input('Column number for the position (1-100):')) if numx == None: print("No column selected") sys.exit() posx=read_excel(fichero,numx) # Lee las posiciones en y numy = int(input('Column number for the height (1-100):')) if numy == None: print("No column selected") sys.exit() altura=read_excel(fichero,numy) # Suavizado de la altura con un filtro Savitzky-Golay smaltura=savgol_filter(altura, 11, 4) # Calculamos la primera derivada pos1=[] der1=[] # Primera derivada for i in range(0,len(posx)-2): pos1.append((posx[i+1]+posx[i])*0.5) der1.append((altura[i+1]-altura[i])/(posx[i+1]-posx[i])) # Suavizado de la primera derivada con un filtro Savitzky-Golay smder1=savgol_filter(der1, 21, 4, mode='nearest') # Calculamos la segunda derivada pos2=[] der2=[] # Segunda derivada der2s=[] # Segunda derivada a partir de la primera suaviazada for i in range(0,len(pos1)-2): pos2.append((pos1[i+1]+pos1[i])*0.5) der2.append((der1[i+1]-der1[i])/(pos1[i+1]-pos1[i])) der2s.append((smder1[i+1]-smder1[i])/(pos1[i+1]-pos1[i])) # Suavizado de la segunda derivada con un filtro Savitzky-Golay smder2=savgol_filter(der2, 21, 4, mode='nearest') # Suavizado de la segunda derivada obtenida a partir de la primera suavizada con un filtro Savitzky-Golay smder2s=savgol_filter(der2s, 21, 4, mode='nearest') # Limites en el eje x de la grafica etlim=['X min', 'X max'] minx=int(min(posx)) maxx=int(max(posx)) ans=raw_input('Limits in position:\nMinimum = '+str(minx) +'? (y/n)\n') if ans != 'y': print('New value for the minimum position?') minx=float(input()) print('Minimum in position: '+str(minx)) ansm=raw_input('Limits in position:\nMaximum = '+str(maxx) +'? (y/n)\n') if ansm != 'y': print('New value for the maximum position?') maxx=float(input()) print('Maximum in position: '+str(maxx)) xlimi=(minx,maxx) # GRAFICAS fig=figure() fig.subplots_adjust(hspace=0.0) # Grafica de la altura ax1=fig.add_subplot(311) # Posicion de la grafica de altura ax1.tick_params(direction='in') # Ticks en la parte interior de la grafica ax1.yaxis.set_ticks_position('both') # Ticks en ambos ejes verticales ax1.xaxis.set_ticks_position('both') # Ticks en ambos ejes horizontales ax1.plot(posx,altura,color="blue") # Dibuja la altura en azul ax1.set_title(fichero) # Nombre del fichero como titulo ax1.set_ylabel("Height (m)") # Etiqueta vertical ax1.set_xticklabels([]) # No pone numeros en el eje horizontal ax1.set_xlim(float(xlimi[0]),float(xlimi[1])) # Limites de la grafica # Grafica de la pendiente ax2=fig.add_subplot(312) # Posicion de la grafica de pendientes ax2.set_xlim(float(xlimi[0]),float(xlimi[1])) # Limites de la grafica ax2.tick_params(direction='in') # Ticks en la parte interior de la grafica ax2.yaxis.set_ticks_position('both') # Ticks en ambos ejes verticales ax2.xaxis.set_ticks_position('both') # Ticks en ambos ejes horizontales ax2.plot(pos1,der1,color="blue") # Dibuja la pendiente en azul ax2.plot(pos1,smder1,'r--',lw=0.6) # Dibuja la pendiente suavizada en rojo con linea de trazos fina #ax2.plot(posx,derv,color="red") # Dibuja la derivada leida del fichero ax2.set_ylabel("Slope\n (First derivative)") # Etiqueta vertical ax2.set_xticklabels([]) # No pone numeros en el eje horizontal # Grafica de la segunda derivada ax3=fig.add_subplot(313) # Posicion de la grafica de segundas derivadas ax3.set_xlim(float(xlimi[0]),float(xlimi[1])) # Limites de la grafica ax3.tick_params(direction='in') # Ticks en la parte interior de la grafica ax3.yaxis.set_ticks_position('both') # Ticks en ambos ejes verticales ax3.xaxis.set_ticks_position('both') # Ticks en ambos ejes horizontales ax3.plot(pos2,der2,color="blue") # Dibuja la segunda derivada en azul ax3.plot(pos2,smder2s,'r--',lw=0.6) # Dibuja la segunda derivada suavizada a partir de datos suavizados en rojo ax3.set_xlabel("Position (m)") # Etiqueta del eje horizontal ax3.set_ylabel("Curvature\n (Second derivative)") # Etiqueta vertical plt.show() sys.exit()