Skip to content

Instantly share code, notes, and snippets.

@mbstacy

mbstacy/fread.py

Created June 8, 2012 14:47
Show Gist options
  • Save mbstacy/2895998 to your computer and use it in GitHub Desktop.
Save mbstacy/2895998 to your computer and use it in GitHub Desktop.
Download ZIP
DDA code for TECO
Raw
fread.py
import sys
from numpy import *
def readobs(filename,headlen=1):
"Read the observational file for the TECO model"
obs=loadtxt(filename,skiprows=headlen)
return obs
def readinit(filename,headlen=0):
"Read initial conditions for TECO model"
initcon=loadtxt(filename,skiprows=headlen)
return initcon
def fout(filename,data,head):
"Output array of data with header from TECO model"
# data=array_str(data)
# f=open(filename,'w')
# f.write(head + '\n')
# f.close
# savetxt(filename,data,fmt='%9.6f', delimiter=' ',header=head)
# savetxt(filename,data,fmt='%1.14e',delimiter=' ')
savetxt(filename,data,fmt='%1.14g',delimiter=' ')
# savetxt(filename,data,fmt='%9.6f', delimiter=' ')
# savetxt(filename,data,fmt='%9.6f')
def readpara(filename,headlen=1):
"Read information from site parameter file"
f=open(filename)
for i in range(headlen):
f.readline()
params=f.readline()
params=params.split()
f.close()
return params
Raw
teco_control_work.py
import fread
import sys
import numpy
import math
import os
from tkfilter import kalman
from teco import *
def blcpg(tair,tavg72,omega,glmax,grmax,gsmax,lai,laimax,laimin,sla,tau_l,bmleaf,bmroot,bmstem,bmplant,rootmax,stemmax,saps,storage,gdd5,stor_use,onset,accumulation,gddonset,nsc,fnsc,nscmin,nscmax,add,tcold,gamma_wmax,gamma_tmax):
"Brians re-write of the plant growth model"
bml=bmleaf
bmr=bmroot
bms=bmstem
bmp=bmplant
bw=2.0
bt=2.0
twarm=30.0
rs0=1.0
rs=bmr/bml
if(bml<nsc/.333*.5):
bml=nsc/.333*.5
if(bmr<nsc/.333*.5):
bmr=nsc/.333*.5
if(bms<nsc/.334*.5):
bms=nsc/.334*.5
stemsap=min(stemmax,saps*bms)
rootsap=bmr
if(stemsap<.001):
stemsap=.001
if(rootsap<.001):
rootsap=.001
if(gdd5>gddonset and onset==0 and storage>stor_use):
onset=1
if(onset==1 and storage>stor_use):
if(lai<laimax):
add=stor_use
storage=storage-add
else:
add=0.0
onset=0.0
if(accumulation<(nscmax+.005*rootsap)):
store=.005*nsc
else:
store=0.0
accumulation=accumulation+store
sps=1-fnsc
sps=max(.001,sps)
st=1/(1+19*math.e**(-.2*tair))
sw=max(.333,.333+omega)
w=min(1,50*omega)
ss=min(1,2*fnsc)
sl_rs=rs/(rs+rs0*(1.5-omega))
sr_rs=(rs0*(1.5-omega))/(rs+rs0*(1.5-omega))
slai=min(1,2.333*(laimax-lai)/(laimax-laimin))
growthl=glmax*bml*st*sw*fnsc*sl_rs*slai*0.45
growthr=grmax*rootsap*st*sw*fnsc*sr_rs*0.5 *0.45
growths=gsmax*stemsap*st*sw*fnsc*0.5*0.5 *0.45
growthl=max(growthl,0)
growthr=max(growthr,0)
growths=max(growths,0)
growthp=growthl+growthr+growths
if(growthp>nsc*.5):
growthl=.5*nsc*growthl/growthp
growthr=.5*nsc*growthr/growthp
growths=.5*nsc*growths/growthp
npp=add+growthl+growthr+growths
if(npp==0):
alpha_l=.333
alpha_w=.333
alpha_r=.333
else:
alpha_l=(growthl+add)/npp
alpha_w=growths/npp
alpha_r=growthr/npp
if(tair>tcold+10):
beta_t=1
else:
if(tair>tcold):
beta_t=(tair-tcold)/10
if(tair<=tcold):
beta_t=0
if(tau_l<8760):
gamma_w=(1.-w)**bw*gamma_wmax
gamma_t=(1-beta_t)**bt*gamma_tmax
gamma_n=0
else:
gamma_w=0
gamma_t=0
gamma_n=1/(tau_l*sw)
if(lai<laimin):
gamma_w=0
gamma_t=0
gamma_n=0
l_fall=bmleaf*.45*min((gamma_w+gamma_t+gamma_n),.99)
results=(stemsap,rootsap,storage,onset,accumulation,sps,store,add,l_fall,npp,alpha_l,alpha_w,alpha_r)
return results
add=0
parafile=sys.argv[1]
outfile=sys.argv[2]
tfilter=int(sys.argv[3])
aint=int(sys.argv[4])
output=0
parameters=fread.readpara(parafile)
(site,vegtype,climfile,NEE_file,out_d,lat,longi,wsmax,wsmin,gddonset,LAIMAX,LAIMIN,rdepth,Rootmax,Stemmax,SapR,SapS,SLA,GLmax,GRmax,Gsmax,stom_n,a1,Ds0,Vcmx0,extkU,xfang,alpha,co2ca,tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass)=parameters
lat=float(lat)
longi=float(longi)
wsmax=float(wsmax)
wsmin=float(wsmin)
gddonset=float(gddonset)
LAIMAX=float(LAIMAX)
LAIMIN=float(LAIMIN)
rdepth=float(rdepth)
Rootmax=float(Rootmax)
Stemmax=float(Stemmax)
SapR=float(SapR)
SapS=float(SapS)
SLA=float(SLA)
GLmax=float(GLmax)
GRmax=float(GRmax)
Gsmax=float(Gsmax)
stom_n=float(stom_n)
a1=float(a1)
Ds0=float(Ds0)
Vcmx0=float(Vcmx0)
extkU=float(extkU)
xfang=float(xfang)
alpha=float(alpha)
co2ca=float(co2ca)
tau_L=float(tau_L)
tau_W=float(tau_W)
tau_R=float(tau_R)
tau_F=float(tau_F)
tau_C=float(tau_C)
tau_Micr=float(tau_Micr)
tau_Slow=float(tau_Slow)
tau_Pass=float(tau_Pass)
climate=fread.readobs(climfile,2)
neeObs_days=list()
neeSim_days=list()
tau_L =tau_L *8760
tau_W =tau_W *8760
tau_R =tau_R *8760
tau_F =tau_F *8760
tau_C =tau_C *8760
tau_Micr=tau_Micr*8760
tau_Slow=tau_Slow*8760
tau_Pass=tau_Pass*8760
GLmax=float(GLmax)/24
GRmax=float(GRmax)/24
Gsmax=float(Gsmax)/24
dayCount = 1
nonValidCount=0
sumObs=0.0
sumSim=0.0
sumObsObs=0.0
sumSimObs=0.0
co2=co2ca
accumulation=0
Q_root=0
year_data=climate[:,0]
day_data=climate[:,1]
hour_data=climate[:,2]
input_data=climate[:,3:]
hour_data=hour_data+1
yr_data=year_data[-1]-year_data[0]+1
inputstep=hour_data[1]-hour_data[0]
lines=len(year_data)-1
tpara=45
#apara=19
#apara=17
apara=11
#apara=9
state=numpy.zeros((lines+1,tpara))
astate=numpy.zeros((lines+1,apara))
updated=numpy.zeros((lines+1,tpara))
cov_mat=0
NEE_obs=numpy.zeros(lines+1)
oerror=[[0],[0],[0]]
#oerror=[0]
#oerror=numpy.zeros(lines+1)
have_obs=False
if(os.path.exists(NEE_file)):
have_obs=True
obs=fread.readobs(NEE_file)
year_NEE=obs[:,0]
day_NEE=obs[:,1]
hour_NEE=obs[:,2]
reader_NEE=obs[:,3:]
for i in range(len(reader_NEE[:,12])):
if reader_NEE[i,10]==-999.:
reader_NEE[i,10]=0
NEE_obs=reader_NEE[:,10]*12./1000./1000000.
yr_NEE=year_NEE[-1]-year_NEE[0]+1
step_NEE=hour_NEE[1]-hour_NEE[0]
Cumol2gram=3600.*step_NEE*1000
year_obs=int(year_NEE[0])
NEE_annual=0.0
thksl=[10.0,10.0,10.0,10.0,10.0,20.0,20.0,20.0,20.0,20.0 ]
FRLEN=[0.30,0.20,0.15,0.15,0.1,0.05,0.05,0.0,0.0,0.0 ]
#thksl=numpy.array(thksl)
#FRLEN=numpy.array(FRLEN)
(pi,tauL,rhoL,rhoS,emleaf,emsoil,Rconst,sigma,cpair,Patm,Trefk,H2OLv0,airMa,H2OMw,chi,Dheat,wleaf,gsw0,Vcmx0,eJmx0,theta,conKc0,conKo0,Ekc,Eko,o2ci,Eavm,Edvm,Eajm,Edjm,Entrpy,gam0,gam1,gam2)=consts()
header_out=' year doy hour carbon_Leaf carbon_Wood carbon_Root carbon_Flitter carbon_Clitter SOM_Miro SOM_SLOW SOM_Pass carbon_canopy carbon_biomass carbon_soil gpp npp nee_simulate nee_observe resp_auto resp_hetero resp_tot ET Transpiration Runoff LatentHeat LAI RootMoist SoilWater soilwater_1 soilwater_2 soilwater_3 soilwater_4 soilwater_5 soilwater_6 soilwater_7 soilwater_8 soilwater_9 soilwater_10 satfrac_1 satfrac_2 satfrac_3 satfrac_4 satfrac_5 satfrac_6 satfrac_7 satfrac_8 satfrac_9 satfrac_10 CO2concentratio AirSpecificHumu rain_div_3600 scale_sw TairPlus273_15'
init_opt=fread.readinit('initial_opt.txt')
wsmax=init_opt[0]
wsmin=init_opt[1]
gddonset=init_opt[2]
LAIMAX=init_opt[3]
LAIMIN=init_opt[4]
rdepth=init_opt[5]
Rootmax=init_opt[6]
Stemmax=init_opt[7]
SapR=init_opt[8]
SapS=init_opt[9]
SLA=init_opt[10]
GLmax=init_opt[11]/24.
GRmax=init_opt[12]/24.
Gsmax=init_opt[13]/24.
a1=init_opt[14]
Ds0=init_opt[15]
Vcmx0=init_opt[16]
alpha=init_opt[17]
tau_L=init_opt[18]*8760.
tau_W=init_opt[19]*8760.
tau_R=init_opt[20]*8760.
tau_F=init_opt[21]*8760.
tau_C=init_opt[22]*8760.
tau_Micr=init_opt[23]*8760.
tau_Slow=init_opt[24]*8760.
tau_Pass=init_opt[25]*8760.
TminV=init_opt[26]
TmaxV=init_opt[27]
ToptV=init_opt[28]
Tcold=init_opt[29]
Gamma_Wmax=init_opt[30]
Gamma_Tmax=init_opt[31]
NSC=init_opt[32]
Q_leaf=init_opt[33]
Q_wood=init_opt[34]
Q_root1=init_opt[35]
Q_root2=init_opt[36]
Q_root3=init_opt[37]
Q_coarse=init_opt[38]
Q_fine=init_opt[39]
Q_micr=init_opt[40]
Q_slow=init_opt[41]
Q_pass=init_opt[42]
S_w_min=init_opt[43]
Q10_h=init_opt[44]
WILTPT=wsmin/100.0
FILDCP=wsmax/100.0
wscontent=WILTPT
fwsoil=1.0
topfws=1.0
omega=1.0
wcl=[FILDCP,FILDCP,FILDCP,FILDCP,FILDCP,FILDCP,FILDCP,FILDCP,FILDCP,FILDCP]
Storage=32.09
times_storage_use=600
stor_use=Storage/times_storage_use
onset=0
duration=0
offset=0
dormancy=1
LAI=LAIMIN
bmstem=Q_wood/0.45
bmroot=(Q_root1+Q_root2+Q_root3)/0.45
bmleaf=Q_leaf/0.45
bmplant=bmstem+bmroot+bmleaf
mult=1
writer=yr_data*mult
num=0
Tavg72=5.0
neeDif_all=0.0
eJmx0=Vcmx0*2.7
m=0
n=0
ta=0
Sps=0.17
cov_run=True
for yr in range(1,int(writer+yr_data)+2):
if((year_data[m%lines]% 4)==0):
if((year_data[m%lines]%100)==0):
if((year_data[m%lines]%400)==0):
idays=366
else:
idays=365
else:
idays=366
else:
idays=365
GDD5=0.0
onset=0
gpp_yr=0.0
NPP_yr=0.0
Rh_yr =0.0
NEE_yr=0.0
neeDif_year=0.0
for days in range(1,idays+1):
StemSap=min(Stemmax,SapS*bmstem)
RootSap=min(Rootmax,SapR*bmroot)
NSCmin=5.
NSCmax=0.05*(StemSap+RootSap+Q_leaf)
#changed by blc, fortran bug
if(days==1 and yr==1):
# ta=(input_data[0,0]-273.15)
ta=0.0
if(ta>0):
GDD5=GDD5+ta
ta=0.0
#canopy and soil model
gpp_t =0.0
transp_t=0.0
Hcanop_d=0.0
evap_t =0.0
Ts=0.0
rain_t=0.0
Trunoff=0.0
RaL=0.0
RaS=0.0
RaR=0.0
Rauto=0.0
neeDif_day=0.0
neeSim_daily=0.0
neeObs_daily=0.0
dtimes=24
for i in range(1,dtimes+1):
ti=i
if(m>lines):
output=1
if(m>lines):
m=0
if(n>lines):
n=0
year=year_data[m]
day=day_data[m]
hour =hour_data[m]
Tair=input_data[m,0]-273.15
rain=input_data[m,6]*3600.
radsol=input_data[m,10]
Qair=input_data[m,2]
wind=math.fabs(input_data[m,4])
co2ca=input_data[m,14] *1.0E-6
Pa_air=101325.0
LWdown=input_data[m,12]
Tsoil=Tair*0.8
windU0=min(1.0,wind)
if(radsol==0):
radsol=.01
eairP=Qair/(Qair+0.62198)*Pa_air
Dair=max(50.01,(esat(Tair)-eairP))
RH=eairP/(eairP+Dair)*100.0
wethr=1
Rnet=0.8*radsol
if(radsol>10):
G=-25
else:
G=0.5
Esoil=.05*radsol
if(radsol<=10):
Esoil=.5*G
Hcrop=0.1
Ecstot=0.1
Anet=0.1
DepH2O=0.2
ta= ta + Tair/24.0
Ts=Ts+Tsoil/24.0
if(NSC<=NSCmin):
fnsc=0.0
if(NSC>=NSCmax):
fnsc=1.0
if(NSC<NSCmax):
if(NSC>NSCmin):
fnsc=(NSC-NSCmin)/(NSCmax-NSCmin)
[gpp,evap,transp,Acanop,Hcanop,Sps]=canopy(fwsoil,topfws,wscontent,LAI,Sps,day,hour,radsol,Tair,Dair,eairP,windU0,rain,wethr,Rnet,G,Esoil,Hcrop,Ecstot,Anet,Tsoil,DepH2O,wsmax,wsmin,lat,co2ca,a1,Ds0,Vcmx0,extkU,xfang,alpha,stom_n,pi,tauL,rhoL,rhoS,emleaf,emsoil,Rconst,sigma,cpair,Patm,Trefk,H2OLv0,airMa,H2OMw,chi,Dheat,wleaf,gsw0,eJmx0,theta,conKc0,conKo0,Ekc,Eko,o2ci,Eavm,Edvm,Eajm,Edjm,Entrpy,gam0,gam1,gam2,TminV,TmaxV,ToptV)
gpp=max(0,gpp)
# call resperation
[RaLeaf,RaStem,RaRoot,Rauto]=respiration(LAIMIN,gpp,Tair,Tsoil,DepH2O,LAI,SLA,bmstem,bmroot,bmleaf,StemSap,RootSap,NSC,fnsc)
# call soilwater
[transp,wcl,evap,runoff,wscontent,fwsoil,topfws,omega,omega_S,WaterR,waters,satfracl]=soilwater(wsmax,wsmin,rdepth,FRLEN,rain,Tair,transp,wcl,Tsoil,RH,thksl,LAI,omega)
omega=omega_S
ET=evap+transp
NA_WC_soil=wscontent
NA_WC_root=WaterR
# call plantgrowth
#bmroot odd after first
[StemSap,RootSap,Storage,onset,accumulation,Sps,store,add,L_fall,npp,alpha_L,alpha_W,alpha_R]=blcpg(Tavg72,Tavg72,omega,GLmax,GRmax,Gsmax,LAI,LAIMAX,LAIMIN,SLA,tau_L,bmleaf,bmroot,bmstem,bmplant,Rootmax,Stemmax,SapS,Storage,GDD5,stor_use,onset,accumulation,gddonset,NSC,fnsc,NSCmin,NSCmax,add,Tcold,Gamma_Wmax,Gamma_Tmax)
NSC=NSC+gpp-Rauto-(npp-add)-store
# call TCS
[Q_leaf,Q_wood,Q_root,Q_fine,Q_coarse,Q_micr,Q_slow,Q_pass,Rh_f,Rh_c,Rh_Micr,Rh_Slow,Rh_Pass]=tcs(Tair,Tsoil,omega,npp,alpha_L,alpha_W,alpha_R,L_fall,tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,Q_leaf,Q_wood,Q_root,Q_fine,Q_coarse,Q_micr,Q_slow,Q_pass,S_w_min,Q10_h)
Rhetero=Rh_f + Rh_c + Rh_Micr + Rh_Slow + Rh_Pass
NEE=Rauto+Rhetero - gpp
Q_soil=Q_micr + Q_slow + Q_pass
bmroot=Q_root/0.45
bmleaf=Q_leaf/0.45
bmstem=Q_wood/0.45
bmplant=bmleaf+bmroot+bmstem
LAI=bmleaf*SLA
# if(writer< yr<(writer+yr_data+1)):
if(output==1):
NA_L=Q_leaf/1000.
NA_W=Q_wood/1000.
NA_R=Q_root/1000.
NA_F=Q_fine/1000.
NA_C=Q_coarse/1000.
NA_Micr=Q_micr/1000.
NA_slow=Q_slow/1000.
NA_pass=Q_pass/1000.
NA_can=NA_L+NA_W*0.7
NA_BM=NA_L+NA_W+NA_R
NA_S=NA_Micr+NA_slow+NA_pass
NA_GPP=gpp/(1000.*3600.)
NA_NPP=npp/(1000.*3600.)
NA_NEE=NEE/(1000.*3600.)
NA_Ra=Rauto/(1000.*3600.)
NA_Rh=Rhetero/(1000.*3600.)
NA_Rt=NA_Ra+NA_Rh
NA_ET=ET
NA_TR=transp
NA_Run=runoff
if(NEE_obs[n]>-500*12/1000/1000000):
NEE_diff = (NA_NEE-NEE_obs[n])*100000000/12
else:
NEE_diff = 0.0
Qle=ET*((2.501-0.00236*Tair)*1000000.0)/3600.
# astate[m]=[Q_leaf,Q_wood,Q_root,Q_fine,Q_coarse,Q_micr,Q_slow,Q_pass,tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,days,i,NEE/(1000*3600)]
# astate[m]=[tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,NEE/(1000*3600)]
# astate[m]=[tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,days,i,NEE/(1000*3600)]
astate[m]=[tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,days,i,NEE]
# astate[m]=[tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,NEE]
if(output==0 and m>=200 and tfilter==1 and have_obs==True and (m%aint==0) and (NEE_obs[m]!=0) and cov_run==False):
[updated,cov_mat]=kalman(astate,(days,ti,NEE_obs[m]*1000*3600),oerror,cov_mat)
# [updated,cov_mat]=kalman(astate,NEE_obs[m],oerror,cov_mat)
# (Q_leaf,Q_wood,Q_root,Q_fine,Q_coarse,Q_micr,Q_slow,Q_pass,tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,tdays,ti,NEE)=updated
# (tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,NEE)=updated
(tau_L,tau_W,tau_R,tau_F,tau_C,tau_Micr,tau_Slow,tau_Pass,tdays,ti,NEE)=updated
# NEE=NEE*1000*3600
state[m]=[Q_leaf,Q_wood,Q_root,Q_fine,Q_coarse,Q_micr,Q_slow,Q_pass,gpp,npp,Rauto,Rhetero,ET,transp,runoff,LAI,WaterR,wscontent,waters[0],waters[1],waters[2],waters[3],waters[4],waters[5],waters[6],waters[7],waters[8],waters[9],satfracl[0],satfracl[1],satfracl[2],satfracl[3],satfracl[4],satfracl[5],satfracl[6],satfracl[7],satfracl[8],satfracl[9],co2ca,Qair,rain,radsol,Tair,Qle,NEE]
m=m+1
# sum day
gpp_t=gpp_t + gpp*(24./dtimes)
transp_t=transp_t + transp*(24./dtimes)
evap_t=evap_t + evap*(24./dtimes)
Hcanop_d=Hcanop_d+Hcanop/(24./dtimes)
Trunoff=Trunoff+runoff
# sum year
gpp_yr=gpp_yr+gpp
NPP_yr=NPP_yr+npp
Rh_yr =Rh_yr +Rhetero
NEE_yr=NEE_yr+NEE
if(output==1):
neeDif_day= neeDif_day + NEE_diff*NEE_diff
neeSim_daily=neeSim_daily+NA_NEE*1000000000./12
neeObs_daily=neeObs_daily+NEE_obs[n]*1000000000./12
n=n+1
if(output==1):
neeObs_days.append(neeObs_daily)
neeSim_days.append(neeSim_daily)
sumObs=sumObs+neeObs_daily
sumSim=sumSim+neeSim_daily
sumObsObs=sumObsObs+neeObs_daily*neeObs_daily
sumSimObs=sumSimObs+neeSim_daily*neeObs_daily
dayCount = dayCount + 1
neeDif_year=neeDif_year + neeDif_day
Tavg72=ta
storage=accumulation
if(cov_run==True and have_obs==True and tfilter==True and m==lines+1):
cov_mat=numpy.cov(astate,rowvar=0)
cov_run=False
stor_use=Storage/times_storage_use
accumulation=0.0
onset=0
neeDif_all=neeDif_all + neeDif_year
outvars=53
output=numpy.array(numpy.zeros((lines+1,outvars)))
output[:,0]=year_data
output[:,1]=day_data
output[:,2]=hour_data
output[:,3]=state[:,0]/1000
output[:,4]=state[:,1]/1000
output[:,5]=state[:,2]/1000
output[:,6]=state[:,3]/1000
output[:,7]=state[:,4]/1000
output[:,8]=state[:,5]/1000
output[:,9]=state[:,6]/1000
output[:,10]=state[:,7]/1000
output[:,11]=output[:,3]+output[:,4]*.07
output[:,12]=output[:,3]+output[:,4]+output[:,5]
output[:,13]=output[:,8]+output[:,9]+output[:,10]
output[:,14]=state[:,8]/(1000*3600)
output[:,15]=state[:,9]/(1000*3600)
output[:,16]=state[:,-1]/(1000*3600)
output[:,17]=NEE_obs
output[:,18]=state[:,10]/(1000*3600)
output[:,19]=state[:,11]/(1000*3600)
output[:,20]=output[:,18]+output[:,19]
output[:,21]=state[:,12]
output[:,22]=state[:,13]
output[:,23]=state[:,14]
output[:,24]=state[:,-2]
output[:,25]=state[:,15]
output[:,26]=state[:,16]
output[:,27]=state[:,17]
output[:,28:38]=state[:,18:28]
output[:,38:48]=state[:,28:38]
output[:,48]=state[:,38]
output[:,49]=state[:,39]
output[:,50]=state[:,40]/3600
output[:,51]=state[:,41]
output[:,52]=state[:,42]+273.15
fread.fout(outfile,output,header_out)
Raw
tkfilter_working.py
import sys
from numpy import *
def kalman(state, obs, oberror):
"updates the forcast given the history and state matrix, observations, and observational error"
oberror=matrix(oberror)
state=matrix(state)
# X=state[:,-1]
X=state[-1,:].T
Z=matrix(obs).T
Pf=cov(state,rowvar=0)
tmp=size(obs)
H=matrix(zeros((tmp,size(X))))
# H=matrix(zeros((size(X),tmp)))
startdim=(shape(H)[1])-tmp
for i in range(0,shape(H)[0]):
# H[i+startdim,i]=1
H[i,i+startdim]=1
R=oberror.T*oberror
PHT=Pf*H.T
D=(H*PHT)+R
# Dn=linalg.pinv(D)
# Dn=linalg.inv(D)
# K=(PHT)*Dn
# K=linalg.lstsq(D.T,PHT.T)[0]
K=linalg.solve(D.T,PHT.T)
K=K.T
tmp=K*H
Pu=(eye(shape(tmp)[0])-tmp)*Pf
Xu=X+(K*(Z-(H*X)))
Xu=list(Xu)
for i in range(len(Xu)):
Xu[i]=float(Xu[i])
return Xu
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment