A Python-only library for estimating global radiation

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 the radest library

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 radest
import os
import sys

if __name__ == "__main__":

    # Get parameters
    if len(sys.argv) != 9:
        print("radlist - Program generating a time series of global radiation estimates")
        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("Copyright 2019 by Mauri Favaron")
        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)])

    # Estimate global radiation
    ra = radest.ExtraterrestrialRadiation(tm,delta, lat, lon, fuse)
    rg = radest.GlobalRadiation(ra, z)

    # 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?

Here’s an answer:

Global radiation, estimated using radest at Milan position (Latitude 45.5°, Longitude 9.16°, Altitude 123m) on year 2019

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.

Global radiation estimated by radest in Milan at four representative days: 21 June (red), 21 Mar and 21 Sep (pink) and 21 Dec (blue) 2019; notice the maxima do not occur exactly at 12:00:00, because of the combined effects of Milan longitude differing by 5.5° by the time zone central longitude, and the Time Equation

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 radest
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("Copyright 2019 by Mauri Favaron")
        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)])
        ra         = radest.ExtraterrestrialRadiation(time_stamp, 1, lat, lon, fuse)

        # 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:

Time equation, as visible in Milan; the mean value, -1402 s, accounts for the difference between the time zone central meridian and site actual meridian. Changes over year take into account the Earth orbit eccentricity and the acceleration around Winter solstice.

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 radest
import os
import sys

if __name__ == "__main__":

    # Get parameters
    if len(sys.argv) != 1:
        print("radtst - Program generating various interesting tables for global radiation estimates")
        print()
        print("Usage:")
        print()
        print("  ./radtst")
        print()
        print("Copyright 2019 by Mauri Favaron")
        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')
        ra  = radest.ExtraterrestrialRadiation(local03, 3600, lat, lon, fuse)
        rg03[ang] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(local06, 3600, lat, lon, fuse)
        rg06[ang] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(local09, 3600, lat, lon, fuse)
        rg09[ang] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(local12, 3600, lat, lon, fuse)
        rg12[ang] = radest.GlobalRadiation(ra, z)

    # 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:

Global radiation, estimated on June 21 (red), December 21 (blue), March 21 and September 21 (pink) at Equator, mean sea level, 00:00:00 UTC, and different longitudes. Note the maxima do not occur all at 180°: the shift, of seasonal nature, is due to the Equation of Time.

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 radest
import os
import sys

if __name__ == "__main__":

    # Get parameters
    if len(sys.argv) != 1:
        print("radtst2 - Program generating various interesting tables for global radiation estimates")
        print()
        print("Usage:")
        print()
        print("  ./radtst2.py")
        print()
        print("Copyright 2019 by Mauri Favaron")
        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)
        ra  = radest.ExtraterrestrialRadiation(time03, 3600, flat, lon, fuse)
        rg03[lat+89] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(time06, 3600, flat, lon, fuse)
        rg06[lat+89] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(time09, 3600, flat, lon, fuse)
        rg09[lat+89] = radest.GlobalRadiation(ra, z)
        ra  = radest.ExtraterrestrialRadiation(time12, 3600, flat, lon, fuse)
        rg12[lat+89] = radest.GlobalRadiation(ra, z)

    # 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:

Radiation change with latitude (longitude = 0°, height = 0m; red line for 21. 06. 2019, blue for 21. 12. 2019, pink lines for 21. 03. 2019 and 21. 09. 2019)

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”

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *