In [1]:
%matplotlib inline
from wmf import wmf
import numpy as np
import pylab as pl
import datetime as dt
Lectura de mapas de direcciones y de elevación:
In [2]:
DEM=wmf.read_map_raster('/media/nicolas/discoGrande/03_SALGAR/raster/dem_salgar.tif',True)
DIR=wmf.read_map_raster('/media/nicolas/discoGrande/03_SALGAR/raster/dir_salgar.tif',True)
wmf.cu.nodata=-9999.0; wmf.cu.dxp=12.7
DIR[DIR<=0]=wmf.cu.nodata.astype(int)
DIR=wmf.cu.dir_reclass(DIR,wmf.cu.ncols,wmf.cu.nrows)
In [3]:
s = wmf.Stream(-76.05,5.99,DEM,DIR)
In [4]:
# Trazado de la cuenca
cuSalgar = wmf.SimuBasin(-75.9808, 5.9647, DEM, DIR, name='Liboriana',
dt = 300, stream = s )
cuSalgar.Save_Basin2Map('/media/nicolas/discoGrande/03_SALGAR/vector/cuenca12m.shp')
cuSalgar.Save_Net2Map('/media/nicolas/discoGrande/03_SALGAR/vector/red12m.shp',umbral=2000)
In [5]:
cuSalgar = wmf.SimuBasin(-75.9808, 5.9647, DEM, DIR, name='Liboriana',
dt = 300, stream = s, modelType = 'hills')
In [6]:
cuSalgar.ncells
Out[6]:
In [7]:
cuSalgar.GetGeo_Cell_Basics()
In [8]:
cuSalgar.set_Geomorphology()
In [9]:
# Establecer Parametros Geomorfo
cuSalgar.GetGeo_Cell_Basics()
cuSalgar.set_Geomorphology(stream_width=cuSalgar.CellLong,umbrales=[200,2000])
In [10]:
cuSalgar.GetQ_Balance(2300)
print cuSalgar.CellQmed[-1]
In [11]:
#Ejemplo turc
cuSalgar.GetQ_Balance(2300, Tipo_ETR=1)
cuSalgar.Plot_basin(cuSalgar.CellETR)
print cuSalgar.CellETR.mean()
In [17]:
cuSalgar.GetQ_Balance(2300, 2)
cuSalgar.Plot_basin(cuSalgar.CellETR)
print cuSalgar.CellETR.mean()
In [15]:
cuSalgar.GetQ_Balance(2300, Tipo_ETR=3)
#cuSalgar.Plot_basin(cuSalgar.CellETR)
print cuSalgar.CellETR.mean()
In [16]:
cuSalgar.CellHeight
Out[16]:
In [8]:
#fi = dt.datetime(2015,5,10)
#ff = dt.datetime(2015,5,19)
#m,p = cuSalgar.rain_radar2basin('radar/','lluvia/Mayo_10_19_12mts',fi,ff,5,
# 'pptint_Salgar_liboriana_','.asc')
Los parámetros físicos del modelo se dividen en velocidades verticales y vel horizontales, las cuales a su ves se dividen en versiones lineales y no lineales de las mismas, en el caso de las velocidades verticales todas las relaciones son lineales, a continuación se describen las variables relacionadas con los flujos que se encuentran en el modelo:
Velocidades Horizontales:
Almacenamientos:
La calibración se compone de 10 parámetros escalares, los cuales son:
In [18]:
# Evaporación en la cuenca estimada por Turc.
Evp=4.658*np.exp(-0.0002*cuSalgar.CellHeight)/288.0
# Infiltración Ks [mm/s]
Ks = 0.00126
# Percolación
Kp = Ks/10.0
# Perdidas se asumen iguales a cero
Kpp = 0.0
# Inclusión de los parámetros en el modelo
Lista=[Evp,Ks,Kp,Kpp]
for pos,var in enumerate(Lista):
cuSalgar.set_PhysicVariables('v_coef',var,pos)
In [19]:
# Velocidades superficial, sub-super, y subte.
v2 = 1.4 * cuSalgar.CellSlope**0.5
v3 = Ks * cuSalgar.CellSlope**0.5
v4 = Kp * cuSalgar.CellSlope**0.5
Lista=[v2,v3,v4]
for pos,var in enumerate(Lista):
cuSalgar.set_PhysicVariables('h_coef',var,pos)
In [20]:
#Calcula variable para el coeficiente horizontal
area = cuSalgar.CellAcum * (12.7**2) #Tamaño de celda al cuadrado
var,w1 = wmf.OCG_param(pend = cuSalgar.CellSlope, area = area)
cuSalgar.set_PhysicVariables('h_coef',var,3)
#Variable del exponente horizontal
cuSalgar.set_PhysicVariables('h_exp',w1,3)
In [21]:
#Genera puntos de control y los transporta al cauce
coordXY = np.array([[-76.0426,5.9858],
[-76.02944,5.9548],
[-76.0159,5.95572]]).T
ids = np.array([200, 201, 202])
cuSalgar.set_Speed_type()
Ids = cuSalgar.set_Control(coordXY,ids,tipo='Q')
order = wmf.models.control[wmf.models.control<>0]
#Puntos de control de humedad
coordHumedad = np.array([[-76.08, 6.00],
[-76.055,5.98],
[-76.03,5.95],
[-76.0014,5.9619]]).T
idsHumedad = np.array([300,301,302,303])
IdsHumedad = cuSalgar.set_Control(coordHumedad,idsHumedad,tipo='H')
wmf.Save_Points2Map(coordXY,ids,'vector/ptos_controlQ.shp')
wmf.Save_Points2Map(coordHumedad,idsHumedad,'vector/ptos_controlH.shp')
A continuación se dan algunas descripciones d elos suelos de salgar, que pueden ser de utilidad para la simulación. También se presentan las ecuaciones usadas para determinar capacidades de almacenamiento, y velocidades de drenaje de acuerdo a la descripción del suelo.
Tipo | Pendiente | Profundidad | Retención | Permeabilidad |
---|---|---|---|---|
Clase III | 12-25 | Media : $0.6m$ | Media | Media |
Clase IV | < 12 | Media : $0.6m$ | Baja | Baja |
Clase VI | 25-30 | Profundo $1.0mts$ | Media | Media |
Clase VII | 30-50 | Superficial $0.3m$ | Muy baja | Lenta |
Clase VIII | > 50 | Superficial $0.2m$ | Muy baja | Lenta |
Estos suelos se describen como suelos de color pardo amarillento a rojizo de textura franco-arcilloza, y cascajos. en general bien drenados y con poca capacidad de retención (poca memoria?). En general también son descritos como suelos profundos (Osorio, 2008) (Unidades de suelo representativas de la Zona Cafetera Colombiana), según la literatura un suelo profundo es más de un metro.
- Capacidad almacenamiento: $0.14 cm/cm$
- $\theta_{cc}$: Humedad a capacidad de campo, se asume igual a: $35\%$
- $\theta_{pmp}$: Humedad en el punto de marchitez, se asume igual a: $21.3\%$
- $\theta_{s}$: Humedad en el punto de saturación, se asume igual a: $47.2\%$
- $K_s$ : Conductividad hidraulica saturada: $4.56 mm/hr$ o $0.00126 mm/seg$
De acuerdo a Saxton, 2006, la conductividad hidraulica del suelo puede variar de acuerdo a las siguientes ecuaciones:
- $K = K_s (\theta / \theta_{s})^{\left[3+(2/\lambda)\right]}$
- $\lambda = 1/B$
- $B = [ln(1500) - ln(33)]/[ln(\theta_{33})-ln(\theta_{1500})]$
- $\theta_{33}$ es la humedad a $33Kpa$, en este caso igual a $37\%$.
- $\theta_{1500}$ es la humedad a $1500Kpa$, en este caso igual a $21.3\%$
Almacenamiento máximo: Se tienen dos almacenamientos máximos, el capilar y el gravitacional, las ecuaciones que los rigen son:
- Capilar: $H_u = Z (\theta_{cc} - \theta_{pmp})$
- Gravitacional : $H_g = Z (\theta_{s} - \theta_{cc})$
La ecuación se describe de la siguiente manera:
$k_{h}=\frac{R*k_{s}S_{0} \Delta x^2}{3H_{g}^2} A ^2$
En este caso $A$ es el área de la sección equivalente en una celda, esta se puede asumir como un porcentaje del área total del flujo que sale por la sección transversal del tanque.
In [22]:
#Mapa de profundidades:
Z = np.zeros(cuSalgar.ncells)
Z[cuSalgar.CellSlope<0.25]=0.6
Z[(cuSalgar.CellSlope>=0.25)&(cuSalgar.CellSlope<0.30)]=1.0
Z[(cuSalgar.CellSlope>=0.30)&(cuSalgar.CellSlope<0.50)]=0.3
Z[cuSalgar.CellSlope>=0.5] = 0.2
#Mapas de almacenamiento maximo capilar y gravitacional
Hu = Z * (0.35 - 0.213) * 1000
Hg = Z * (0.47 - 0.35) * 1000
#cuSalgar.Plot_basin(vec=Z,ruta='imagenes12/ProfEfectivaSuelo.png')
In [32]:
cuSalgar.Plot_basin(Hu,colorbarLabel='a) Capillary maximum storage ($H_u$) [$mm$]', colorTable = pl.get_cmap('viridis'),
ruta='/home/nicolas/Dropbox/Articulos_Prep/Salgar/Figura4_1_Hu.png')
Out[32]:
In [33]:
cuSalgar.set_PhysicVariables('capilar',Hu,0)
cuSalgar.set_PhysicVariables('gravit',Hg,0)
wmf.models.retorno = 1.0
La ecuación se describe de la siguiente manera:
$k_{h}=\frac{R*k_{s}S_{0} \Delta x^2}{3H_{g}^2} A ^2$
En este caso $A$ es el área de la sección equivalente en una celda, esta se puede asumir como un porcentaje del área total del flujo que sale por la sección transversal del tanque.
In [34]:
# Ejemplo de solucionador de onda cinematica por multiples remplazos
def sol(sm,coef,expo,dx,dt):
v = 0.00127
for i in range(4):
Area = sm*0.9 / (dx + v*dt)
vn = coef * Area ** expo
v = (2*vn + v)/3.0
return Area,v
#Calculo de coeficiente por kubota y ensayo de como funciona
ksh=(Ks*cuSalgar.CellSlope*(12.7**2.0))/(3*(Hg*0.9/1000.0)**2)
speed = 1
area,v = sol(80,100,2.0,12.7,300.0)
print v
print area
print 1.0 - 30.0/(v*300.0+30.0)
In [38]:
cuSalgar.Plot_basin(ksh, colorTable=pl.get_cmap('viridis'), colorbarLabel='b) Kubota & Sivapalan coefficient',
ruta = '/home/nicolas/Dropbox/Articulos_Prep/Salgar/Figura4_2_ksh.png')
Out[38]:
In [ ]:
pos = np.random.choice(cuSalgar.ncells,10000)
h,b = np.histogram(ksh,bins=np.linspace(0,200,7))
h = h.astype(float); h = h/h.sum()
b = (b[:-1]+b[1:])/2.0
pl.plot(b,h,lw = 2)
pl.grid(True)
pl.xlabel('$\\frac{R*k_{s}S_{0} \\Delta x^2}{3H_{g}^2}$ $[mm/seg]$',size=15)
pl.ylabel('pdf $[\\%]$',size=15)
pl.savefig('imagenes12/PDF_Coef_khs_12m.png',bbox_inches='tight')
In [ ]:
# Establece el coeficiente y el exponente para la version no lineal de sub-superficial
cuSalgar.set_PhysicVariables('h_coef',ksh,1)
cuSalgar.set_PhysicVariables('h_exp',2.0,1)
wmf.models.speed_type[1] = 2
A continuación se realiza una ejecución ensayo del modelo, ojo esta es una prueba inicial puede fallar facilmente
La calibración se compone de 10 parámetros escalares, los cuales son:
In [17]:
#Actualizar estados
for pos,i in enumerate(Res['Storage']):
cuSalgar.set_Storage(i+0.01,pos)
In [18]:
#Reiniciar estados
#cuSalgar.set_Storage(Hu,0)
#cuSalgar.set_Storage(2,0)
In [19]:
#wmf.models.storage = Res['Storage']
wmf.models.separate_fluxes = 1
#cuSalgar.set_Storage(10,3)
#cuSalgar.set_Storage(0.1,2)
#cuSalgar.set_Storage(0,0)
Npasos = 300
Inicio = 1
Calibracion = [0.1, 2.7, 0.8, 0.0, 1.0, 1.5, 0.01, 0.996, 1.0, 1.0]
ruta_lluvia = 'lluvia/Mayo_18_19_12mts.bin'
Res = cuSalgar.run_shia(Calibracion,ruta_lluvia,Npasos,Inicio)#ruta_storage='storage_17_18.bin')
S = wmf.read_mean_rain('lluvia/Mayo_18_19_12mts.hdr',Npasos,Inicio)
In [21]:
cuSalgar.Save_SimuBasin('/media/nicolas/discoGrande/binCuencas/cuencaSalgar_2_12mts.nc',
ruta_dem = '/media/nicolas/discoGrande/raster/dem_salgar.tif',
ruta_dir = '/media/nicolas/discoGrande/raster/dir_salgar.tif')
In [37]:
wmf.plot_sim_single(Res['Fluxes'][3][:2],figsize = (9,4.5),mrain=S.values,Dates = S.index.to_pydatetime(),
ruta='imagenes12/Separacion_Flujos.png',ids = ['Run-off','Sub-Sup'])
In [29]:
cuSalgar.Save_SimuBasin('/media/nicolas/discoGrande/binCuencas/cuencaSalgar_12mts.nc',
ruta_dem = '/media/nicolas/discoGrande/raster/dem_salgar.tif',
ruta_dir = '/media/nicolas/discoGrande/raster/dir_salgar.tif')
In [30]:
cuSalgar.Plot_basin(wmf.models.retorned,ZeroAsNaN='si')
In [20]:
#Caudal simulado
wmf.plot_sim_single(Res['Qsim'],figsize = (9,4.5),mrain=S.values,Dates = S.index.to_pydatetime(),
ids=[0]+ids.tolist(),)#ruta='imagenes12/Qsim_18_Condiciones0.png')
In [29]:
#Mapa Humedad Simulado
Humedad = (Res['Storage'][2]+Res['Storage'][0])/(Hu+Hg)
cuSalgar.Plot_basin(vec=Humedad,ZeroAsNaN='si',)#ruta='imagenes12/Mapa_Humedad_18_2:10.png')
In [30]:
#Serie de humedad simulada
pos = np.where(wmf.models.control_h<>0)[1]
HumedadRel=[]
for i,j,p in zip(Res['Humedad'],[300,301,302,303],pos):
HumedadRel.append(i/(Hu[p]+Hg[p]))
HumedadRel=np.array(HumedadRel)
wmf.plot_sim_single(HumedadRel,mrain=S.values,Dates=S.index.to_pydatetime(),
ylabel='Soil Moisture $[\\theta]$',ids=[300,301,302,303],figsize=(10,5),)
#ruta = 'imagenes12/HumedadSimulada_18.png')
In [31]:
cuSalgar.set_Storage(Hu*0.8,0)
QsimOrig = Res['Qsim'][0]
In [32]:
Res = cuSalgar.run_shia(Calibracion,ruta_lluvia,Npasos,Inicio)
S = wmf.read_mean_rain('lluvia/Mayo_18_19_12mts.hdr',Npasos,Inicio)
QsimEdit = Res['Qsim'][0]50.0
In [33]:
QsimJoin = np.vstack([QsimOrig,QsimEdit])
wmf.plot_sim_single(QsimJoin,legend=False,Dates=S.index.to_pydatetime(),figsize=(10,5),
ruta = 'imagenes12/Qsim_comparado_Hu',mrain=S.values)
Se observa como varía la humedad en la cuenca en diferentes puntos de control, en el punto de contol 2000 (ubicado en la zona alta de la cuenca) se observa como durante el evento del 18 de Mayo se da una saturación importante 1 hora antes del evento, el suelo se encontraba con un valor de humedad de alrededor del 55%, este valor se mantiene, y cuando sucede el aguacero de las 2:00 se hace más crudo.
In [75]:
#Actualizar estados
for pos,i in enumerate(Res['Storage']):
cuSalgar.set_Storage(i,pos)
cuSalgar.set_Storage(0,0)
Npasos = 225
Inicio = 1
Res = cuSalgar.run_shia(Calibracion,ruta_lluvia,Npasos,Inicio)
S = wmf.read_mean_rain('lluvia/Mayo_18_19_12mts.hdr',Npasos,Inicio)
In [77]:
#Series de humedad en puntos de control
pos = np.where(wmf.models.control_h<>0)[1]
HumedadRel=[]
for i,j,p in zip(Res['Humedad'],[300,301,302,303],pos):
HumedadRel.append(i/(Hu[p]+Hg[p]))
HumedadRel=np.array(HumedadRel)
wmf.plot_sim_single(HumedadRel,mrain=S.values,Dates=S.index.to_pydatetime(),
ylabel='Soil Moisture $[\\theta]$',ids=[300,301,302,303],figsize=(10,5),
ruta = 'imagenes12/HumedadSimulada_Antes18.png')
#Mapa de humedad antes del evento
Humedad = (Res['Storage'][2]+Res['Storage'][0])/(Hu+Hg)
cuSalgar.Plot_basin(vec=Humedad,ZeroAsNaN='si',ruta='imagenes12/Mapa_Humedad_Antes18.png')