forked from faluhong/ATC-and-DTC-Code
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDTC_GOT09.py
89 lines (62 loc) · 3.5 KB
/
DTC_GOT09.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
import numpy as np
import math
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
pi=np.pi
# This is the example of GOT09 model writing with Python 3.7
def GOT09(time, *param):
T0, Ta, tm, ts, DeltaT, tao = param[0], param[1], param[2], param[3],param[4],param[5]
theta = pi / 12 * (time - tm)
theta_s = pi / 12 * (ts - tm)
mask_LT = time < ts
mask_GE = time >= ts
Re, H = 6371, 8.43
cos_sza = math.sin(delta / 180 * pi) * math.sin(Latitude / 180 * pi) + math.cos(delta / 180 * pi) * math.cos(
Latitude / 180 * pi) * np.cos(theta)
cos_sza_s = math.sin(delta / 180 * pi) * math.sin(Latitude / 180 * pi) + math.cos(delta / 180 * pi) * math.cos(
Latitude / 180 * pi) * math.cos(theta_s)
sin_sza_s = math.sqrt(1 - cos_sza_s * cos_sza_s)
cos_sza_min = math.sin(delta / 180 * pi) * math.sin(Latitude / 180 * pi) + math.cos(delta / 180 * pi) * math.cos(
Latitude / 180 * pi)
m_val = -Re / H * cos_sza + np.sqrt(pow((Re / H * cos_sza), 2) + 2 * Re / H + 1)
m_sza_s = -Re / H * cos_sza_s + math.sqrt(pow((Re / H * cos_sza_s), 2) + 2 * Re / H + 1)
m_min = -Re / H * cos_sza_min + math.sqrt(pow((Re / H * cos_sza_min), 2) + 2 * Re / H + 1)
sza_derive_s = math.cos(delta / 180 * pi) * math.cos(Latitude / 180 * pi) * math.sin(theta_s) / math.sqrt(
1 - cos_sza_s * cos_sza_s)
m_derive_s = Re / H * sin_sza_s - pow(Re / H, 2) * cos_sza_s * sin_sza_s / math.sqrt(
pow(Re / H * cos_sza_s, 2) + 2 * Re / H + 1)
k1 = 12 / pi / sza_derive_s
k2 = tao * cos_sza_s * m_derive_s
k3=DeltaT/Ta*(cos_sza_min)/np.exp((m_min-m_sza_s)*tao)
k = k1 * (cos_sza_s-k3) / (sin_sza_s + k2)
temperature1 = T0 + Ta * cos_sza[mask_LT] * np.exp(tao * (m_min - m_val[mask_LT])) / cos_sza_min
temp1 = math.exp(tao * (m_min - m_sza_s)) / cos_sza_min
temp2 = np.exp(-12 / pi / k * (theta[mask_GE] - theta_s))
temperature2 = (T0+DeltaT) + (Ta * cos_sza_s * temp1-DeltaT) * temp2
temperature = np.concatenate((temperature1, temperature2))
return temperature
# Example data
Latitude=40.3581
Longitude=-3.9803
DOY=139
# DTC model starts from the sunrise to the next sunrise, the values greater than 24 refer to the time in the next day
time_DTC=np.arange(4.75,28.5,1)
temperature_DTC=np.array([283,286,290,295,301,306,309,313,314,314,312,309,305,301,297,294,293,291,290,288,287,286,285,284],dtype=np.float)
# delta is the solar declination
delta = 23.45 * np.sin(2 * pi / 365.0 * (284 + DOY))
# omega is the duration of daytime
omega = 2.0 / 15 * math.acos(-math.tan(Latitude / 180.0 * pi) * math.tan(delta / 180.0 * pi)) * 180.0 / pi
sunrisetime = 12 - omega / 2
sunsettime=12+omega/2
# Setting the initial value
p0 = [np.average(temperature_DTC), np.max(temperature_DTC) - np.min(temperature_DTC), 13.0, sunsettime-1,0,0.01]
# Solving the free parameters of GOT09 model. Bounds and max_nfev are optional
popt, pcov = curve_fit(GOT09, time_DTC, temperature_DTC, p0, bounds = ([200, 0, 8,12,-5,0], [400, 60, 17, 23,5,1]), max_nfev = 10000)
print(popt)
temperature_modelling = GOT09(time_DTC, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
print('RMSE of GOT09 model:',np.sqrt(np.mean(np.square(temperature_DTC-temperature_modelling))))
plt.title('Example of GOT09 model')
plt.plot(time_DTC,temperature_DTC,'.g',label='LST observations')
plt.plot(time_DTC, temperature_modelling, 'r', label='DTC modelling results')
plt.legend(loc='best')
plt.show()