############################################################ # # # 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. # # # ############################################################ import sys from pylab import * import xlrd import easygui as eg import numpy as np from scipy.signal import savgol_filter # Ventana informacion eg.msgbox("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",title='Problem 1', ok_button='Continue') # Eleccion de fichero seleccion = ["*.xls","*.xlsx","*.ods"] fichero=eg.fileopenbox(msg="Please choose a file",title="Data file",default="*",filetypes=seleccion) 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 = eg.integerbox(msg='Column number for the position (1-100):', title='Position', default=4, lowerbound=1, upperbound=100, image=None) if numx == None: print("No column selected") sys.exit() posx=read_excel(fichero,numx) # Lee las posiciones en y numy = eg.integerbox(msg='Column number for the height (1-100):', title='Height', default=3, lowerbound=1, upperbound=100, image=None) if numy == None: print("No column selected") sys.exit() altura=read_excel(fichero,numy) ## Lee las derivadas #numz = eg.integerbox(msg='Column number for the slope (1-100):', # title='Slope', # default='', # lowerbound=1, # upperbound=100, # image=None) #if numz == None: # print("No column selected") # sys.exit() #derv=read_excel(fichero,numz) # 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)) xlimi=eg.multenterbox(msg='Limits in position', title='X limits',fields=etlim,values=(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()