stoporad to Python

[scotthavens]

SMRF uses stoporad from IPW to estimate the clear sky radiation over snow. This PR implements a full Python version of 3 main IPW utilities: elevrad, toporad and stoporad.

TL;DR

• elevrad is within the bit resolution noise of the IPW elevrad
• toporad is within the bit resolution noise of the IPW toporad
• stoporad has differences of +/- a few W/m2 due to slightly different calculations for horizon and albedo
• New file structure in envphys to split up the monster radiation.py into manageable pieces that are also more intuiative (see below)
• Fixed an issue with the loadTopo, the new basin_setup has the y direction decreasing which creates a negative dy. This negative dy is passed to the gradient function which flips everything up/down.
• The threaded aspect of distribute.solar was simplified as seperate threads for the IPW calls are no longer needed.
• The Python version is 20x faster than the IPW version

New file structure in envphys

Takes the old radiation.py and seperates into a solar and thermal package with organized files.

smrf/envphys/
├── albedo.py
├── constants.py
├── core
│   ├── dewpt.c
│   ├── envphys_c.c
│   ├── envphys_c.h
│   ├── envphys_c.pyx
│   ├── envphys.h
│   ├── iwbt.c
│   └── topotherm.c
├── precip.py
├── snow.py
├── solar
│   ├── cloud.py
│   ├── ipw.py
│   ├── irradiance.py
│   ├── model.py
│   ├── toporad.py
│   ├── twostream.py
│   └── vegetation.py
├── storms.py
├── sunang.py
├── thermal
│   ├── clear_sky.py
│   ├── cloud.py
│   ├── topotherm.py
│   └── vegetation.py
└── vapor_pressure.py

stoporad validation

Validation of the new functions were performed against the IPW versions. Everything was held constant between the two as much as possible. A notebook containing the comparisons are in notebooks/toporad_comparison.ipynb. All comparisons were done in the Lakes.

elevrad

elevrad is the spatial version of twostream which was implemented in PR #123. Comparison with the IPW version show bit resolution noise from the 8-bit IPW output image. All values are in W/m2.

toporad

toporad uses the output of elevrad and adjusts for topography affects. Comparison with the IPW version show bit resolution noise from the 8-bit IPW output image. All values are in W/m2.

stoporad

stoporad in IPW is a shell script that runs multiple functions, including elevrad and toporad. The stoporad function calls are specific for either the visible or infrared wavelengths. There are differences in both the visible and infrared images but are minimal. These differences are caused by using a different horizon and albedo function before calling elevrad and toporad. For the visible, the largest differences are single pixels that are different from the two horizon outputs. All values are in W/m2.

Gold changes

RME stations

All variables are 0 expect net solar with minimal changes.

RME HRRR

All variables are 0 expect net solar with minimal changes.

Lakes HRRR

All variables are 0 except net solar and thermal. The cause for the large changes were an issue with loadTopo as the topo.y variable is decreasing, hence dy is negative which is passed to the gradient calculation. The RME y increases so this wasn't noticed.

[jomey]

Wondering if there is an opportunity to use much of the constants or standard solar calculations from the astropy package?
https://www.astropy.org

[jomey]

Really nice improvement to the readability of the radiation component!

One thought I want to throw out there is to add a ticket to improve performance by vectorizing a lot of those functions? Without looking much into it, most calculations look like they can be easily converted with numba for instance.

[scotthavens]

Wondering if there is an opportunity to use much of the constants or standard solar calculations from the astropy package?
https://www.astropy.org

Perhaps. At one time I did look around at other packages but they seem to either be tailored to astronomy or solar panel applications. I currently haven't found a package that is specific the the problem that the the toporad functions solve. However, the constants should almost all be the same but should we have another dependency just for constants?

[jomey]

However, the constants should almost all be the same but should we have another dependency just for constants?

I agree that a dependency for constants is not worthwhile. Mostly wondering more about the standard solar angle calculations , since those seem standard.

[scotthavens]

I think it would be something to look at in the future and see how they compare. Since you'll be working with albedo which is a calculation dependent on the solar angles, might be worth looking into it.

[jomey]

I think it would be something to look at in the future and see how they compare. Since you'll be working with albedo which is a calculation dependent on the solar angles, might be worth looking into it.

#166

[scotthavens]

Really nice improvement to the readability of the radiation component!

One thought I want to throw out there is to add a ticket to improve performance by vectorizing a lot of those functions? Without looking much into it, most calculations look like they can be easily converted with numba for instance.

I've tried to keep them as vectorized as possible to take advantage of numpy. Numba could definitely help on the larger basins or high resolution DEM's perhaps. As it stands, the straight python versions of these functions are 20x faster than the IPW versions. A 4 month simulations I did for a 500x1000 domain completed in like an hour instead of 3-4 hours (and not even threaded). So we're already seeing huge speed ups from this.

The tests freaked me out a little bit because the tests went from ~90 seconds to ~13 seconds. I thought that some of the tests were being skipped but they weren't.

[jomey]

Really nice improvement to the readability of the radiation component!
One thought I want to throw out there is to add a ticket to improve performance by vectorizing a lot of those functions? Without looking much into it, most calculations look like they can be easily converted with numba for instance.

I've tried to keep them as vectorized as possible to take advantage of numpy. Numba could definitely help on the larger basins or high resolution DEM's perhaps. As it stands, the straight python versions of these functions are 20x faster than the IPW versions. A 4 month simulations I did for a 500x1000 domain completed in like an hour instead of 3-4 hours (and not even threaded). So we're already seeing huge speed ups from this.

The tests freaked me out a little bit because the tests went from ~90 seconds to ~13 seconds. I thought that some of the tests were being skipped but they weren't.

That is a really good boost. Looking forward to apply that to ERW.
Part of the improvement probably is from removing all the calls to a subprocess, which can be expensive with many at once.

[jomey]

General question: Is 50 degrees consider sunset?

[scotthavens]

This is 50 W/m2. There can be some unintended artifacts at low solar values for calculating cloud factor (measured / model) especially if the measured doesn't quite go to 0 at the same time as the modeled.

[scotthavens]

And many of the IPW functions have a lot of file IO. For stoporad there are probably 6 or so temporary files that are generated and removing the IO can be fairly significant.