Least Squares Estimation in Python: Unterschied zwischen den Versionen
Aus Wiki
| Zeile 1: | Zeile 1: | ||
print 'calculate trend' | print 'calculate trend' | ||
db1 = Database("DAHITI_Waterlevels") | db1 = Database("DAHITI_Waterlevels") | ||
| − | + | db1.connect() | |
| − | + | parameter = ["date","height","std"] | |
| − | + | tablename = "tid_%06d" % (int(dahiti_id)) | |
| − | + | data = db1.select(tablename,parameter,orderby=["date ASC"]) | |
| − | + | print 'Points',len(data) | |
| − | + | A = np.ones((len(data),2)) | |
| − | + | b = np.ones((len(data),1)) | |
| − | + | for i in range(0,len(data)): | |
| − | + | d = data[i] | |
| − | + | jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0) | |
| − | + | b[i] = float(d['height']) | |
| − | + | A[i][1] = jd | |
| − | + | db1.disconnect() | |
| − | + | ||
| − | + | A = np.matrix(A) | |
| − | + | b = np.matrix(b) | |
| − | + | x = ((A.T * A).I * A.T) * b | |
| − | + | v = (A * x - b) | |
| − | + | sigma2 = float(((v.T * v) / (b.size - 2))) | |
| − | + | Qxx = (A.T*A).I | |
| − | + | Sigma_xx = sigma2 * Qxx | |
| − | + | ||
| − | + | error = np.sqrt(float(Sigma_xx[1,1])) * 365.25 * 1000.0 | |
| − | + | trend = float(x[1])*365.25*1000.0 | |
| − | + | ||
| − | + | start_ts = julianDay2000((data[0]['date'][0:4]),(data[0]['date'][5:7]),(data[0]['date'][8:10]),12,0,0) | |
| − | + | end_ts = julianDay2000((data[-1]['date'][0:4]),(data[-1]['date'][5:7]),(data[-1]['date'][8:10]),12,0,0) | |
| − | + | years = (end_ts - start_ts) / 365.25 | |
| − | + | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | print "Sea Level Trend: %.3f +- %.3f mm/year in %.1f years" % (trend,error,years) | |
| + | #~ print " Min. SLA: %.3f mm (%s)" % (start_sla,data[0]['date']) | ||
| + | #~ print " Max. SLA: %.3f mm (%s)" % (end_sla,data[-1]['date']) | ||
| + | |||
| + | print path+"/www_timeseries_trend.dat" | ||
| + | output = open(path+"/www_timeseries_trend.dat","w") | ||
| + | output.write(julianDayDate(start_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * start_ts + x[0])))) | ||
| + | output.write(julianDayDate(end_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * end_ts + x[0])))) | ||
| + | output.close() | ||
| + | |||
| + | trend = "Trend: %.1f \261 %.1f mm/year (%.1f years)" % (trend,error,years) | ||
| + | output = open('/tmp/www_timeseries_trend.legend','w') | ||
| + | output.write('G -0.025i\n') | ||
| + | output.write('N 1\n') | ||
| + | output.write('S 0.12c c 0.1c darkgreen 0.25p 0.3c '+trend+'\n') | ||
| + | output.close() | ||
Aktuelle Version vom 8. Februar 2021, 10:20 Uhr
print 'calculate trend'
db1 = Database("DAHITI_Waterlevels")
db1.connect()
parameter = ["date","height","std"]
tablename = "tid_%06d" % (int(dahiti_id))
data = db1.select(tablename,parameter,orderby=["date ASC"])
print 'Points',len(data)
A = np.ones((len(data),2))
b = np.ones((len(data),1))
for i in range(0,len(data)):
d = data[i]
jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0)
b[i] = float(d['height'])
A[i][1] = jd
db1.disconnect()
A = np.matrix(A)
b = np.matrix(b)
x = ((A.T * A).I * A.T) * b
v = (A * x - b)
sigma2 = float(((v.T * v) / (b.size - 2)))
Qxx = (A.T*A).I
Sigma_xx = sigma2 * Qxx
error = np.sqrt(float(Sigma_xx[1,1])) * 365.25 * 1000.0
trend = float(x[1])*365.25*1000.0
start_ts = julianDay2000((data[0]['date'][0:4]),(data[0]['date'][5:7]),(data[0]['date'][8:10]),12,0,0)
end_ts = julianDay2000((data[-1]['date'][0:4]),(data[-1]['date'][5:7]),(data[-1]['date'][8:10]),12,0,0)
years = (end_ts - start_ts) / 365.25
print "Sea Level Trend: %.3f +- %.3f mm/year in %.1f years" % (trend,error,years)
#~ print " Min. SLA: %.3f mm (%s)" % (start_sla,data[0]['date'])
#~ print " Max. SLA: %.3f mm (%s)" % (end_sla,data[-1]['date'])
print path+"/www_timeseries_trend.dat"
output = open(path+"/www_timeseries_trend.dat","w")
output.write(julianDayDate(start_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * start_ts + x[0]))))
output.write(julianDayDate(end_ts)['date'][0:10]+' %.3f\n' % (float((x[1] * end_ts + x[0]))))
output.close()
trend = "Trend: %.1f \261 %.1f mm/year (%.1f years)" % (trend,error,years)
output = open('/tmp/www_timeseries_trend.legend','w')
output.write('G -0.025i\n')
output.write('N 1\n')
output.write('S 0.12c c 0.1c darkgreen 0.25p 0.3c '+trend+'\n')
output.close()
