Add files via upload

DTC_GEM.py

@@ -0,0 +1,114 @@
+import numpy as np
+import math
+from matplotlib import pyplot as plt
+from scipy.optimize import curve_fit
+from scipy.integrate import quad
+pi=np.pi
+
+# This is the example of GEM-type model (including GEM-eta and GEM-sigma models) writing with Python 3.7
+# Both models have four parameters
+
+def GEM_eta(time,*param):
+    P, T_Ave, eta0, eta1 = param[0], param[1], param[2], param[3]
+    sum = 0
+    h1 = eta0 + eta1 * time_DTC
+    for i in range(1, FourierNumber + 1):
+        M = 1.0 / np.sqrt(i * OmegaCal * P * P + np.sqrt(2 * i * OmegaCal) * P * h1 + h1 * h1)
+        FourierA = 1.0 / pi * quad(FourierKernalA, X1, X2, args=(i))[0]
+        FourierB = 1.0 / pi * quad(FourierKernalB, X1, X2, args=(i))[0]
+        Fi = np.arctan(P * np.sqrt(i * OmegaCal) / (np.sqrt(2.0) * h1 + P * np.sqrt(i * OmegaCal)))
+        Gt = FourierA * np.cos(i * OmegaCal * time - Fi) + FourierB * np.sin(i * OmegaCal * time - Fi)
+
+        sum = sum + Gt * M
+
+    LSTReturn = T_Ave +  sum
+    return LSTReturn
+
+def GEM_sigma(time,*param):
+
+    P, T_Ave, h1, sigma =param[0], param[1], param[2], param[3]
+    sum=0
+    for i in range(1,FourierNumber+1):
+
+        M=1.0/np.sqrt(i*OmegaCal*P*P+np.sqrt(2*i*OmegaCal)*P*h1+h1*h1)
+        FourierA = 1.0 / pi * quad(FourierKernalA, X1, X2, args=(i))[0]
+        FourierB = 1.0 / pi * quad(FourierKernalB, X1, X2, args=(i))[0]
+        Fi=np.arctan(P*np.sqrt(i*OmegaCal)/(np.sqrt(2.0)*h1+P*np.sqrt(i*OmegaCal)))
+        Gt=FourierA*np.cos(i*OmegaCal*time-Fi)+FourierB*np.sin(i*OmegaCal*time-Fi)
+
+        sum=sum+Gt*M
+
+    LSTReturn=T_Ave+sigma*(time-43200)+sum
+    return LSTReturn
+
+def FourierKernalA(x,n):
+    CosSZA=np.cos(delta/180*pi)*np.cos(Latitude/180*pi)*np.cos(x-pi)+np.sin(delta/180*pi)*np.sin(Latitude/180*pi)
+    AtmosphereTransmittance=1-0.2*np.sqrt(1.0/CosSZA)
+    return (1-albedo)*SolarConstant*CosSZA*AtmosphereTransmittance*np.cos(n*x)
+
+def FourierKernalB(x,n):
+    CosSZA=np.cos(delta/180*pi)*np.cos(Latitude/180*pi)*np.cos(x-pi)+np.sin(delta/180*pi)*np.sin(Latitude/180*pi)
+    AtmosphereTransmittance=1-0.2*np.sqrt(1.0/CosSZA)
+    return (1-albedo)*SolarConstant*CosSZA*AtmosphereTransmittance*np.sin(n*x)
+
+
+
+
+# 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)
+
+
+# Basic settings of GEM-type model
+FourierNumber=10
+SolarConstant=1367
+albedo=0.15
+
+# OmegaEarth represents the angular velocity of the earth
+OmegaEarth=2*pi/86400
+DayNumber=1
+OmegaCal=OmegaEarth/DayNumber
+
+
+# delta is the solar declination
+delta = 23.45 * np.sin(2 * pi / 366.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
+
+# Calculating the integrating range
+# TimeChangeOne and TimeChangTwo cloud be regarded as the time range of daytime (represented by seconds)
+TimeChangeOne= 43200 - np.arccos((0.04 - np.sin(delta / 180 * pi) * np.sin(Latitude / 180 * pi)) / (np.cos(delta / 180 * pi) * np.cos(Latitude / 180 * pi))) / OmegaEarth
+TimeChangTwo= 86400 - TimeChangeOne
+#X1 and X2 represent the intergrating range of FourierKernalA and FourierKernalA
+X1 = TimeChangeOne * OmegaEarth
+X2 = TimeChangTwo * OmegaEarth
+
+
+# Setting the initial value
+p0 = [1000 ,np.average(temperature_DTC) ,10, 0]
+
+# Solving the free parameters of GEM-eta model. Bounds and max_nfev are optional
+popt, pcov = curve_fit(GEM_eta, time_DTC*3600, temperature_DTC, p0, bounds=([100, 200, -100, -100], [5000, 350, 100, 100]),max_nfev=10000)
+temperature_modelling = GEM_eta(time_DTC*3600, popt[0], popt[1], popt[2], popt[3])
+print (popt)
+
+# Solving the free parameters of GEM-sigma model. Bounds and max_nfev are optional
+# popt, pcov = curve_fit(GEM_sigma, time_DTC*3600, temperature_DTC, p0, bounds = ([100, 200, -100, -100], [5000, 350, 100, 100]), max_nfev = 10000)
+# temperature_modelling = GEM_sigma(time_DTC*3600, popt[0], popt[1], popt[2], popt[3])
+# print (popt)
+
+print('RMSE of GEM model:',np.sqrt(np.mean(np.square(temperature_DTC-temperature_modelling))))
+
+plt.title('Example of GEM-type 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()
+

DTC_GOT01.py

@@ -0,0 +1,61 @@
+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 GOT01 model writing with Python 3.7
+
+def GOT01(time, *param):
+
+    T0, Ta, tm, ts,deltaT = param[0], param[1], param[2], param[3], param[4]
+    k = omega / (pi * np.tan(pi / omega * (ts - tm)))
+    temperature = np.zeros(np.shape(time),dtype=np.float)
+    Mask=time<ts
+    temperature[Mask] = T0 + Ta * np.cos(pi / omega * (time[Mask] - tm))
+
+    #You can obtain the INA08 model by replacing the exponential function with hyperbolic function
+    temperature[~Mask] = T0 +deltaT + [Ta * np.cos(pi / omega * (ts - tm))-deltaT] * np.exp(-(time[~Mask]- ts) / k)
+
+    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]
+
+# Solving the free parameters of GOT09 model. Bounds and max_nfev are optional
+popt, pcov = curve_fit(GOT01, time_DTC, temperature_DTC, p0, bounds = ([200, 0, 8,12,-5], [400, 60, 17, 23,5]), max_nfev = 10000)
+print(popt)
+
+temperature_modelling = GOT01(time_DTC, popt[0], popt[1], popt[2], popt[3],popt[4])
+
+print('RMSE of GOT01 model:',np.sqrt(np.mean(np.square(temperature_DTC-temperature_modelling))))
+
+plt.title('Example of GOT01 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()
+
+

DTC_GOT09.py

@@ -0,0 +1,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()
+
+

DTC_GOT09_dT_tao.py

@@ -0,0 +1,93 @@
+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-dT-tao model writing with Python 3.7
+
+# GOT09-dT-tao model is the recommended four-parameter DTC model
+# GOT09-dT-tao model is modified from GOT09 model by fixing deltaT as 0 and tao as 0.01
+def GOT09_dT_tao(time, *param):
+    T0, Ta, tm, ts = param[0], param[1], param[2], param[3]
+
+    theta = pi / 12 * (time - tm)
+    theta_s = pi / 12 * (ts - tm)
+    mask_LT = time < ts
+    mask_GE = time >= ts
+
+    #fixing tao (represent the total atmospheric optical thickness) as 0.01
+    tao = 0.01
+
+    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
+    k = k1 * cos_sza_s / (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 + Ta * cos_sza_s * temp1 * 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]
+
+# Solving the free parameters of GOT09 model. Bounds and max_nfev are optional
+popt, pcov = curve_fit(GOT09_dT_tao, time_DTC, temperature_DTC, p0, bounds = ([200, 0, 8,12], [400, 60, 17, 23]), max_nfev = 10000)
+print(popt)
+
+temperature_modelling = GOT09_dT_tao(time_DTC, popt[0], popt[1], popt[2], popt[3])
+
+print('RMSE of GOT09-dT-tao model:',np.sqrt(np.mean(np.square(temperature_DTC-temperature_modelling))))
+
+plt.title('Example of GOT09-dT-tao 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()
+
+