In [37]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
%matplotlib inline
from scipy.fftpack import ifft, fft, fftfreq, fft2, ifft2, fftn, ifftn, fftshift
In [38]:
from scipy import interpolate
import smopy
In [39]:
running=np.genfromtxt("/Users/juan/Downloads/running_clean_nov.gpx",delimiter=",")
In [40]:
RunLon=running[:,0]
RunLat=running[:,1]
minLat=min(RunLat)
minLon=min(RunLon)
maxLat=max(RunLat)
maxLon=max(RunLon)
wi=max(abs(maxLat-minLat),abs(maxLon-minLon))
map = smopy.Map((4.635,-74.089,4.654,-74.070),z=15)
x, y = map.to_pixels(RunLat,RunLon)
ax = map.show_mpl(figsize=(8,8))
ax.plot(x,y,"r",lw=3)
#plt.savefig("/Users/juan/Documents/Trabajo/Andes/2015-V/MetodosComputacionalesLaboratorio/2015-V/actividades/Talleres/Taller6/trip.png")
plt.show()
In [52]:
ar=np.array([xRun,yRun,zRun])
In [59]:
R=6371 # Radio de la Tierra en Km
latBog=4.5981/180.*np.pi # Cambiar a radianes
lonBog=-74.0758/180.*np.pi # Cambiar a radianes
thetaBog=lonBog
phiBog=np.pi/2.-latBog # Medir respecto al polo norte
basetheta=np.array([-np.sin(thetaBog),np.cos(thetaBog),0.]) # Calcular el vector unitario en dirección oriente
basephi=-np.array([-np.cos(thetaBog)*np.cos(phiBog),np.sin(thetaBog)*np.cos(phiBog),-np.sin(phiBog)]) # Ahora el que va hacia el norte
In [65]:
phiRun=running[:,0]/180.*np.pi # Convertir a radianes
thetaRun=np.pi/2.-running[:,1]/180.*np.pi # Convertir a radianes
xRun=R*np.cos(thetaRun)*np.sin(phiRun) # Rx
yRun=R*np.sin(thetaRun)*np.sin(phiRun) # Ry
zRun=R*np.cos(phiRun) # Rz
# Calcular los desplazamientos
xRun-=xRun[0]
yRun-=yRun[0]
zRun-=zRun[0]
# Ahora proyectar sobre los vectores unitarios N y O
xLocal=np.zeros(len(xRun)) # En dirección oriente
yLocal=np.zeros(len(yRun)) # En dirección norte
for i in range(len(xLocal)):
vecGlobe=np.array([xRun[i],yRun[i],zRun[i]])
xLocal[i]=sum(vecGlobe*basetheta)
yLocal[i]=sum(vecGlobe*basephi)
t=running[:,3]+running[:,4]/60.+running[:,5]/3600.-running[0,3]-running[0,4]/60-running[0,5]/3600.
In [66]:
plt.figure(figsize=(20,10))
plt.plot(xLocal,yLocal,"-")
plt.axis('equal')
plt.grid('on')
plt.xlabel('x/km')
plt.ylabel('y/km')
plt.show()
In [124]:
# calcular la distancia recorrida
difX=np.diff(xLocal)
difY=np.diff(yLocal)
difR=np.sqrt(difX**2+difY**2)
cumR=np.cumsum(difR)
newt=t[:-1:]
plt.figure(figsize=(15,8))
plt.subplot(2,1,1)
plt.plot(newt,cumR,"o-",ms=0.5)
plt.xlim(0,1.6)
plt.xlabel("t/h")
plt.ylabel("d/km")
plt.grid('on')
# ahora calcular rapideces
# calcular la distancia recorrida
deli=20 # salto para suavizar
difX=np.diff(xLocal[::deli])
difY=np.diff(xLocal[::deli])
difR=np.sqrt(difX**2+difY**2)
cumR=np.cumsum(difR)
newt2=t[::deli]
#Calcular rapideces
dift=np.diff(newt2)
print "Δt/s=",np.mean(newt2)*3600
speed=difR/dift
newt3=newt2[:-1:]
# Una interpolación con un cubic spline
cubic_speed = interpolate.interp1d(newt3,speed,kind='cubic')
x_fine=np.linspace(min(newTime2),max(newTime2),200)
# Ahora la gráfica de rapidez
plt.subplot(2,1,2)
plt.plot(newt3,speed,"o--",ms=2.5)
plt.plot(x_fine,cubic_speed(x_fine),"r",lw=2,alpha=0.5)
plt.xlim(0,max(newt2[:-1:]))
plt.ylim(0,15.)
plt.xlabel("t/h")
plt.ylabel("v/kph")
plt.grid('on')
plt.show()
In [130]:
# calcular la distancia recorrida
difX=np.diff(xLocal)
difY=np.diff(yLocal)
difR=np.sqrt(difX**2+difY**2)
cumR=np.cumsum(difR)
newt=t[:-1:]
fun=interpolate.interp1d(cumR,newt)
for i in range(1,10):
print str(i)+" km | " + str(round(fun(i),2))+" h"
In [131]:
# calcular la distancia recorrida
difX=np.diff(xLocal)
difY=np.diff(yLocal)
difR=np.sqrt(difX**2+difY**2)
cumR=np.cumsum(difR)
newt=t[:-1:]
with plt.xkcd():
plt.figure(figsize=(15,8))
plt.subplot(1,1,1)
plt.plot(newt,cumR)
plt.xlim(0,1.6)
plt.xlabel("t/h")
plt.ylabel("d/km")
plt.grid('on')
plt.title("A walk in the park")
#plt.savefig("/Users/juan/Documents/Trabajo/Andes/2015-V/MetodosComputacionalesLaboratorio/2015-V/actividades/talleres/Taller6/tripdist.png")
plt.show()