-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsolarcalc.py
360 lines (293 loc) · 10.9 KB
/
solarcalc.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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# Copyright © SolarCalc Project Contributors
# https://github.com/cgq-qgc/solarcalc
#
# This code is licensed under the terms of the MIT License as published
# by the Open Source Initiative. For more details, see the MIT license at
# https://opensource.org/licenses/MIT/.
#
# This file is a derivative work of codes taken from the SolarCalc 1.0 program,
# licensed under the terms of the Creative Commons CCZero License as published
# by the Creative Commons nonprofit organization. For more details, see the
# CC0 1.0 License at https://creativecommons.org/publicdomain/zero/1.0/.
#
# Source:
# USDA Agricultural Research Service (2019). SolarCalc 1.0. United States
# Department of Agriculture. https://data.nal.usda.gov/dataset/solarcalc-10/.
# Accessed 2021-11-19.
# -----------------------------------------------------------------------------
"""
Global Solar Radiation Estimator
Predicts hourly global solar radiation on a horizontal surface from
daily average temperature extremes and total precipitation records.
"""
from __future__ import annotations
import os.path as osp
import numpy as np
import pandas as pd
__version__ = "1.1"
__project_url__ = "https://github.com/cgq-qgc/solarcalc"
# To avoid "RuntimeWarning: overflow encountered in power" warnings when
# calculating tau**m, which become very large just outside of the
# sunrise and sunset times. This is not relevant here since Sd values
# before sunrise and after sunset are forced to 0 anyway.
np.seterr(over='ignore')
def load_demo_climatedata() -> pd.DataFrame:
"""
Return a pandas dataframe containing climate data to run
a demo of solarcalc.
"""
return pd.read_csv(
osp.join(osp.dirname(__file__), 'demo_climatedata.csv'),
parse_dates=['datetime'],
index_col='datetime',
dtype={'tamin_degC': 'float',
'tamax_degC': 'float',
'ptot_mm': 'float'}
)
def calc_eqn_of_time(dayofyear: int) -> float:
"""
Calculate the Equation of Time correction (typically a 15-20 minutes
correction depending on calendar day).
Source
------
Equation 11.4 in Campbell, G.S. and J.M. Norman (1998). An Introduction to
Environmental Biophysics.
Parameters
----------
dayofyear : int
Day of year.
Returns
-------
float
Calculated Equation of Time value in hours.
"""
f = (279.575 + 0.98565 * dayofyear) * np.pi / 180
return (
-104.7 * np.sin(f) +
596.2 * np.sin(2 * f) +
4.3 * np.sin(3 * f) +
-12.7 * np.sin(4 * f) +
-429.3 * np.cos(f) +
-2.0 * np.cos(2 * f) +
19.3 * np.cos(3 * f)
) / 3600
def calc_long_corr(lon_dd: float) -> float:
"""
Calculate the longitudinal correction.
Parameters
----------
lon_dd : float
Longitude of fieldsite in decimal degrees. Negative value if West
of meridian and positive value if East of meridian.
Returns
-------
float
Longitudinal correction in hours.
"""
# Calculate the local standard time meridian.
lstm = lon_dd - lon_dd % (np.sign(lon_dd) * 15)
# Calculate the longitudinal correction in hours.
long_corr = (lon_dd - lstm) * (24 / 360)
return long_corr
def calc_solar_declination(dayofyear: int) -> float:
"""
Calculate the solar declination angle.
Source
------
Equation 11.2 in Campbell, G.S. and J.M. Norman (1998). An Introduction to
Environmental Biophysics.
Parameters
----------
dayofyear : int
Day of year.
Returns
-------
temp2 : float
Solar declination angle in radians.
"""
temp2 = (
278.97 + 0.9856 * dayofyear +
1.9165 * np.sin((356.6 + 0.9856 * dayofyear) * np.pi / 180)
)
temp2 = np.sin(temp2 * np.pi / 180)
temp2 = np.arcsin(0.39785 * temp2)
return temp2
def calc_zenith_angle(lat_rad: float, solar_dec: float, time: float,
solar_noon: float) -> float:
"""
Zenith angle approximation.
Source
------
Equation 11.1 in Campbell, G.S. and J.M. Norman (1998). An Introduction to
Environmental Biophysics.
Parameters
----------
lat_rad : float
Latitude in radians.
solar_dec : float
Solar declination angle in radians.
time : float
Time of the day in hours.
solar_noon : float
Solar noon value in hours.
Returns
-------
zenith_angle : float
Zenith angle approximation in radians.
"""
hour_angle = 15 * (solar_noon - time) * np.pi / 180
zenith_angle = np.arccos(
np.sin(lat_rad) * np.sin(solar_dec) +
np.cos(lat_rad) * np.cos(solar_dec) * np.cos(hour_angle)
)
return zenith_angle
def calc_halfdaylength(solar_dec: float, lat_rad: float,
twilight: bool = False) -> float:
"""
Calculation of 1/2 solar day length.
Source
------
Equation 11.6 in Campbell, G.S. and J.M. Norman (1998). An Introduction to
Environmental Biophysics.
Parameters
----------
solar_dec : float
Solar Declination in radians.
lat_rad : float
Latitude of site (field) in radians.
twilight : bool
Whether to include twilight in the half day lenght calculation.
Returns
-------
float
1/2 day length in hours.
"""
psi = 96 if twilight else 90
A = np.cos(psi * np.pi / 180) - np.sin(lat_rad) * np.sin(solar_dec)
B = np.cos(lat_rad) * np.cos(solar_dec)
return np.arccos(A/B) * (180 / np.pi) / 15
def calc_solar_rad(lon_dd: float, lat_dd: float, alt: float,
climate_data: pd.DataFrame) -> pd.DataFrame:
"""
Predict hourly global solar radiation on a horizontal surface.
Parameters
----------
lon_dd : float
Site longitude in degree decimal.
lat_dd : float
Site latitude in degree decimal.
alt : float
Altitude of the site in m.
climate_data : pd.DataFrame
A pandas dataframe containing the consecutive (no missing value) daily
meterological data that are required for the calculatations in the
following columns :
* tamax_degC : daily maximum temperature in Celcius.
* tamin_degC : daily minimum temperature in Celcius.
* ptot_mm : daily total precipitations in mm.
The index of the dataframe must contain the date for each daily
reading in a pandas DatetimeIndex.
Returns
-------
A pandas dataframe that contain the hourly values estimated by the model
in the following columns:
* solar_rad_W/m2 : hourly global solar radiation on a horizontal
surface in W/m2.
* deltat_degC : daily temperature total variation.
* tau : atmospheric transmittance.
The index of the dataframe contains the date and time for each hourly
reading in a pandas DatetimeIndex.
"""
# Convert rain[day of year] = 1 if rain and rain = 0 if no rain.
rain = (climate_data['ptot_mm'] > 1).values.astype(int)
deltaT = (climate_data['tamax_degC'] -
climate_data['tamin_degC']).abs().round(5).values
# Estimate the atmospheric transmittance (tau) based on the
# input climate data.
n = len(climate_data)
tau = np.zeros(n)
for i in range(n):
tau[i] = 0.70 # Clear sky value as given in Gates (1980)
# If it is raining then assume it is overcast (cloud cover).
if rain[i] == 1:
tau[i] = 0.40 # from Liu and Jordan (1960)
# If it has been raining for two days then even
# darker (denser cloud cover).
if i > 0 and rain[i] == 1 and rain[i - 1] == 1:
tau[i] = 0.30
# Assign pre-rain days to 80% of tau value ?
if i < (n - 1) and rain[i] == 0 and rain[i + 1] == 1:
tau[i] = 0.60
# If airtemperature rise is less than 10 --> lower tau value
# unless by poles
if np.abs(lat_dd) < 60:
if deltaT[i] <= 10 and deltaT[i] != 0:
tau[i] = tau[i] / (11 - deltaT[i])
# Expand rain, deltaT and tau from daily to hourly time frame.
rain = np.repeat(rain, 24)
deltaT = np.repeat(deltaT, 24)
tau = np.repeat(tau, 24)
# Prepare the dataframe that will be used to store the results and that
# will be returned as the output of this function.
dataframe = pd.DataFrame(
[],
index=pd.date_range(
start=climate_data.index[0],
end=climate_data.index[-1] + pd.Timedelta('23H'),
freq='H'))
dataframe.index.name = 'datetime'
dayofyear = dataframe.index.dayofyear.values
time = dataframe.index.hour.values
# Calculate the longitudinal correction to solar noon.
LC = calc_long_corr(lon_dd)
# Calculate the correction for Equation of Time.
ET = calc_eqn_of_time(dayofyear)
# Calculate solar noon value.
solarnoon = 12 - LC - ET
# Calculate the solar declination angle
solarD = calc_solar_declination(dayofyear)
# Calculate the length of solar day (excluding twilight time).
halfdaylength = calc_halfdaylength(solarD, np.radians(lat_dd))
# Calculate the sunrise and sunset time.
sunrise = solarnoon - halfdaylength
sunset = solarnoon + halfdaylength
# Calculate the Zenith Angle.
zenith_angle = calc_zenith_angle(
np.radians(lat_dd), solarD, time, solarnoon)
# Calculate the atmospheric pressure at the observation site using
# Equation 3.7 in Campbell and Norman (1998).
Pa = 101.3 * np.exp(-1 * alt / 8200)
# Calculate the optical air mass number using Equation 11.12 in
# Campbell and Norman (1998).
m = Pa / 101.3 / np.cos(zenith_angle)
Spo = 1360 # Solar constant in W/m2
# Calculate the hourly beam irradiance on a horizontal surface (Sb)
# using Equations 11.8 and 11.11 in Campbell and Norman (1998).
Sb = Spo * np.power(tau, m) * np.cos(zenith_angle)
Sb[(time < sunrise) | (time > sunset)] = 0
# Calculate the diffuse sky irradiance on horizontal plane (Sd)
# using Equation 11.13 in Campbell and Norman (1998).
Sd = 0.3 * (1 - np.power(tau, m)) * Spo * np.cos(zenith_angle)
Sd[(time < sunrise) | (time > sunset)] = 0
# The global solar radiation on a horizontal surface is the sum of the
# horizontal direct beam (Sb) and diffuse radiation (Sd).
St = Sb + Sd
# Check for never daylight northern latitudes.
St[St < 0] = 0
# Fill NaN values with zeros.
St[pd.isnull(St)] = 0
# Fill the output dataframe with the results.
dataframe['solar_rad_W/m2'] = np.round(St, 2)
dataframe['deltat_degC'] = deltaT
dataframe['tau'] = tau
return dataframe
if __name__ == '__main__':
climate_data = load_demo_climatedata()
solar_rad = calc_solar_rad(
lon_dd=-76.4687209,
lat_dd=56.5213541,
alt=100,
climate_data=climate_data)
print(solar_rad)