# The problem

Rough estimates of global solar radiation are useful in many cases, like for example checking a radiometer’s response, or predicting evaporations-transpiration.

Computing these estimates is possible, using one of the many existing models. The pbl_met library, for example, allows doing this in two different ways (three, if we consider the old MPDA method).

Existing estimation software is most often written in Fortran (the pbl_met!), in C/C++ or other compiled language, which is fine if performance is at stake. There exist use cases, however, where performance is not that big issue, and maybe functionality is of greater interest; this is the area of use of Python and other interpreted languages. In them, support for numerical computations may be very strong, but specific physical modelling may be less available.

Library radest attempt filling this specific gap.

Using radest is a simple task: you just have to import it. Wel’ll see soon how to do, using a tiny example.

For the library (and the example code) to work, your Python installation should support version 3.6 or more, and have a working installation of NumPy. For using the library on your own, an elementary knowledge of NumPy, in addition to plain Python, is recommended. In particular, an understanding of ‘datetime64’ and ‘timedelta64’, the NumPy equivalents of ‘datetime’ standard Python types, is highly helpful.

Nothing really difficult, by the way.

And here’s the example.

#!/usr/bin/env python3

import numpy as np
import os
import sys

if __name__ == "__main__":

# Get parameters
if len(sys.argv) != 9:
print()
print("Usage:")
print()
print("  ./radlist <Lat> <Lon> <Fuse> <Height> <Initial_ISO_DateTime> <Delta_Time> <Num_Steps> <Out_File>")
print()
print("where")
print()
print("  <Lat>, the latitude, is positive and increasing northwards (decimal degrees)")
print("  <Lon>, longitude, is positive and increasing eastwards (decimal degrees)")
print("  <Fuse> is an integer indicating the hours displacement from GMT")
print("  <Height> is the height above ground (m)")
print("  <Initial_ISO_DateTime> is a string like 2019-03-08T12:00:00, indicating a date and time")
print("  <Delta_Time> is the (intgeer!) time step (s)")
print("  <Num_Steps> is the number of time steps (greater than 0)")
print("  <Out_File> is the output file name")
print()
print("This is open-source code, covered by the MIT license")
print()
sys.exit(1)
lat   = float(sys.argv[1])
lon   = float(sys.argv[2])
fuse  = int(sys.argv[3])
z     = float(sys.argv[4])
first = np.datetime64(sys.argv[5], 's')
delta = int(sys.argv[6])
n     = int(sys.argv[7])
out   = sys.argv[8]

# Generate the time sequence corresponding to the desired time range
tm = np.array([first + np.timedelta64(i*delta, 's') for i in range(n)])

# Print results
f = open(out, "w")
f.write("date, Rg\n")
for i in range(len(tm)):
f.write("%s, %f\n" % (str(tm[i]), rg[i]))
f.close()

# Some examples from Milan

Using radest allows to estimate the clear-sky global radiation. What would we then expect on a year, 2019 for example?

The preceding plot shows quite well the envelope of maximal radiation one could expect on a day. But, what happens on a specific day? We can see in the next plot, where data from approximate season beginnings are shown.

# Is the time equation visible in radest?

The preceding plot suggests something is in fact discernible in radest data. But, what its magnitude is, and what’s its time evolution? In other term: if we measure the shift of the maximum radiation occurrence from the local “noon”, do we see something really similar to the Time Equation?

In order to check this, I’ve built a tiny application (which, incidentally, shows some data processing in NumPy/Python.

This is the code:

#!/usr/bin/env python3

import numpy as np
import os
import sys

if __name__ == "__main__":

# Get parameters
if len(sys.argv) != 6:
print("timeq - Program generating an estimate of the Time Equation")
print()
print("Usage:")
print()
print("  ./timeq <Lat> <Lon> <Fuse> <Year> <Out_File>")
print()
print("where")
print()
print("  <Lat>, the latitude, is positive and increasing northwards (decimal degrees)")
print("  <Lon>, longitude, is positive and increasing eastwards (decimal degrees)")
print("  <Fuse> is an integer indicating the hours displacement from GMT")
print("  <Year> is an integer designating the year for which processing should occur")
print("  <Out_File> is the output file name")
print()
print("This is open-source code, covered by the MIT license")
print()
sys.exit(1)
lat   = float(sys.argv[1])
lon   = float(sys.argv[2])
fuse  = int(sys.argv[3])
year  = int(sys.argv[4])
out   = sys.argv[5]

# Main loop: Iterate over all days in dsired year, and process them
initial_date = np.datetime64("%4.4d-01-01T00:00:00" % year, 's')
final_date   = np.datetime64("%4.4d-01-01T00:00:00" % (year+1), 's')
num_days     = ((final_date - initial_date) / np.timedelta64(1, 'D')).astype(int)
days         = np.array([np.datetime64('2000-01-01', 'D') for i in range(num_days)])
time_eq      = np.array([0.0 for i in range(num_days)])
for day in range(num_days):

start_of_day = initial_date + np.timedelta64(day, 'D')
print(start_of_day)

# Compute extraterrestrial radiation for all seconds in this day
time_stamp = np.array([start_of_day + np.timedelta64(sec, 's') for sec in range(86400)])

# Find the index corresponding to the maximum extraterrestrial radiation
# and use it to calculate the time stamp at which it occurs
idx_max_ra  = np.argmax(ra)
time_max_ra = time_stamp[idx_max_ra]

# Compute difference from noon, and save it
noon = start_of_day + np.timedelta64(12, 'h')
days[day]    = start_of_day
time_eq[day] = (time_max_ra - noon) / np.timedelta64(1, 's')

# Print results
f = open(out, "w")
f.write("date, tm.eqn\n")
for i in range(len(days)):
f.write("%s, %f\n" % (str(days[i]), time_eq[i]))
f.close()

And this is the plot:

As we can see, there is something similar to the Time Equation.

# Other interesting effects

One of the thing the radiation estimate is able to do is, as with any mathematical model, eliciting questions, and providing an instrument for answering them – to the degree of accuracy permitted by the model.

And these two are questions, maybe useful in a high school sciences course:

1. How does radiation change with latitude, as we moe instantaneously over a given meridian?
2. Let’s go back to Milan: how does radiation change with altitude on a given day?

The first question can be solved using code which scans all longitudes on a given hour. The code is a bit trickier than usual, to cope with change of local hour on different fuses:

#!/usr/bin/env python3

import numpy as np
import os
import sys

if __name__ == "__main__":

# Get parameters
if len(sys.argv) != 1:
print()
print("Usage:")
print()
print()
print("This is open-source code, covered by the MIT license")
print()
sys.exit(1)

# Build data sets and print them, in sequence
lat    = 0.
lon    = 0.
fuse   = 0
z      = 0.
time03 = np.array([np.datetime64('2019-03-21T00:00:00', 's')])
time06 = np.array([np.datetime64('2019-06-21T00:00:00', 's')])
time09 = np.array([np.datetime64('2019-09-21T00:00:00', 's')])
time12 = np.array([np.datetime64('2019-12-21T00:00:00', 's')])
rg03   = np.array([0. for ang in range(360)])
rg06   = np.array([0. for ang in range(360)])
rg09   = np.array([0. for ang in range(360)])
rg12   = np.array([0. for ang in range(360)])
out    = 'Equator_2019-03-21.csv'
for ang in range(360):
lon  = ang
fuse = int(ang/24)
if fuse == 360:
fuse = 0
local03 = time03 + np.timedelta64(fuse, 'h')
local06 = time06 + np.timedelta64(fuse, 'h')
local09 = time09 + np.timedelta64(fuse, 'h')
local12 = time12 + np.timedelta64(fuse, 'h')

# Print results
f = open(out, "w")
f.write("date, rg.03, rg.06, rg.09, rg.12\n")
for i in range(len(rg03)):
f.write("%d, %f, %f, %f, %f\n" % (i, rg03[i], rg06[i], rg09[i], rg12[i]))
f.close()

Here is the plotted result:

The second question is somewhat simpler to answer, as the fuse is now just the same. Here is the code:

#!/usr/bin/env python3

import numpy as np
import os
import sys

if __name__ == "__main__":

# Get parameters
if len(sys.argv) != 1:
print()
print("Usage:")
print()
print()
print("This is open-source code, covered by the MIT license")
print()
sys.exit(1)

# Build data sets and print them, in sequence
lat    = 0.
lon    = 0.
fuse   = 0
z      = 0.
time03 = np.array([np.datetime64('2019-03-21T12:00:00', 's')])
time06 = np.array([np.datetime64('2019-06-21T12:00:00', 's')])
time09 = np.array([np.datetime64('2019-09-21T12:00:00', 's')])
time12 = np.array([np.datetime64('2019-12-21T12:00:00', 's')])
lon = 0
fuse = 0
lats   = [i-89 for i in range(179)]
rg03   = np.array([0. for ang in range(179)])
rg06   = np.array([0. for ang in range(179)])
rg09   = np.array([0. for ang in range(179)])
rg12   = np.array([0. for ang in range(179)])
out    = 'Equator_Latitudes.csv'
for lat in lats:
flat = float(lat)

# Print results
f = open(out, "w")
f.write("lat, rg.03, rg.06, rg.09, rg.12\n")
for i in range(len(rg03)):
f.write("%d, %f, %f, %f, %f\n" % (i-89, rg03[i], rg06[i], rg09[i], rg12[i]))
f.close()

And this is the corresponding plot:

# Conclusions

So, the radest library and its routines are here!

For you to use, if you desire.

The MIT license I’ve adopted is really permissive: you are free to use my software in both open-source and commercial applications, basically with no obligation to thank me.

But of course, a citation, or a public mention, of me and the library, is welcome and, from your side, a sign of kindness.

If you have questions, both related to the physics and use, you may ask to me, maybe even by posting a comment to this article.

Have a nice (and, of course, radiation-filled) day.

## 2 risposte a “A Python-only library for estimating global radiation”

1. Marta Salvati ha detto:

Thanks (as always) for sharing your work!