Least Squares Estimation in Python - Versionsgeschichte

93.104.14.225 am 8. Februar 2021 um 10:20 Uhr

2021-02-08T10:20:06Z

93.104.14.225 am 8. Februar 2021 um 10:18 Uhr

2021-02-08T10:18:30Z

93.104.14.225 am 8. Februar 2021 um 10:17 Uhr

2021-02-08T10:17:54Z

93.104.14.225: Die Seite wurde neu angelegt: „print 'calculate trend' db1 = Database("DAHITI_Waterlevels") db1.connect() parameter = ["date","height","std"] tablename = "tid_%06d" % (int(dah…“

2021-02-08T10:16:51Z

Die Seite wurde neu angelegt: „print 'calculate trend' db1 = Database("DAHITI_Waterlevels") db1.connect() parameter = ["date","height","std"] tablename = "tid_%06d" % (int(dah…“

Neue Seite

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]
# print (d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10])
jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0)
# print jd
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()

@@ Zeile 1: / Zeile 1: @@
   print  'calculate trend'
   db1 = Database("DAHITI_Waterlevels")
-			db1.connect()
+ db1.connect()
-			parameter = ["date","height","std"]
+ parameter = ["date","height","std"]
-			tablename = "tid_%06d" % (int(dahiti_id))
+ tablename = "tid_%06d" % (int(dahiti_id))
-			data = db1.select(tablename,parameter,orderby=["date ASC"])
+ data = db1.select(tablename,parameter,orderby=["date ASC"])
-			print 'Points',len(data)
+ print 'Points',len(data)
-			A = np.ones((len(data),2))
+ A = np.ones((len(data),2))
-			b = np.ones((len(data),1))
+ b = np.ones((len(data),1))
-			for i in range(0,len(data)):
+ for i in range(0,len(data)):
-				d = data[i]
+     d = data[i]
-		#		print (d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10])
+     jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0)
-				jd = julianDay2000((d['date'][0:4]),(d['date'][5:7]),(d['date'][8:10]),12,0,0)
+     b[i] = float(d['height'])
-		#		print jd
+     A[i][1] = jd
-				b[i] = float(d['height'])
+ db1.disconnect()
-				A[i][1] = jd
-			db1.disconnect()
+ A = np.matrix(A)
+ b = np.matrix(b)
-			A = np.matrix(A)
+ x = ((A.T * A).I * A.T) * b
-			b = np.matrix(b)
+ v = (A * x - b)
-			x = ((A.T * A).I * A.T) * b
+ sigma2 = float(((v.T * v) / (b.size - 2)))
-			v = (A * x - b)
+ Qxx = (A.T*A).I
-			sigma2 = float(((v.T * v) / (b.size - 2)))
+ Sigma_xx = sigma2 * Qxx
-			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
-			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)
-			start_ts = julianDay2000((data[0]['date'][0:4]),(data[0]['date'][5:7]),(data[0]['date'][8:10]),12,0,0)
+ years = (end_ts - start_ts) / 365.25
-			end_ts = julianDay2000((data[-1]['date'][0:4]),(data[-1]['date'][5:7]),(data[-1]['date'][8:10]),12,0,0)
-</syntaxhighlight>
+ print "Sea Level Trend: %.3f +- %.3f mm/year in %.1f years" % (trend,error,years)