search

USDA-ARS-NWRC / smrf

SMRF was designed to increase the flexibility of taking measured weather data, or atmospheric models, and distributing the data across a watershed.
Other
12 stars 4 forks source link

Solar radiation calculations to python: sun angle, solar and twostream #123

Closed scotthavens closed 4 years ago

scotthavens commented 5 years ago

There were three major functions in the radiation module that called IPW functions for sun angle, solar and twostream. This PR brings them all into the Python world with some extra support.

This text (csv) file show the difference for every hour in this test for water year 2016 in Reynolds Mountain East. The code used to generate the figures and text file are in tests/test_radiation.py. radiation_comparison.txt

TL;DR

Sun angle

The sun angle calculation was done with a subprocess call to the IPW sunang function. The calculations were extracted from IPW and converted to Python. Tests were added to ensure that the sun angle functions were returning correct values and improved on the functionality. If latitude and longitude are provided as a numpy array, sunang will calculate the sun angle for every point.

At a point, single date

The first test replicates the example in the IPW man page for sunang with no difference between the example and the Python version. The values must be truncated in Python to match the print to stdout that SMRF reads from. The truncation can be turned off if needed or if a numpy array is passed.

At a point, whole water year

Comparing over a full water year for RME shows very little differences, most likely in the precision of the print statements. The azimuth below shows a difference in azimuth (degrees) of less that 0.01 degrees. smrf_sunang_azimuth

The RME zenith angle in degrees is similar, with absolute values less than 0.0002 degrees. smrf_sunang_zenith

Affects to net solar radiation

With the small differences, the net solar radiation is not the same as the gold standard for RME causing the tests to fail. However, looking at the differences most values are about 0 with a maximum of 0.175 W/m2 which is insignificant for the magnitude of solar radiation. Accepting this PR will require an update in the gold standard. smrf_net_solar_diff

Sun angle for a whole basin

With the sun angle in Python, we can now pass latitude and longitude as an array. This will return an array of cosz and azimuth that can be provided to functions when they are ready. This will require rework of illumination angle and solar modules. Below is an example from the Tuolumne River Basin. cosz does not vary widely across the basin but the azimuth does have some differences of around 1 degree from one side to the other. smrf_sunang_cosz_tuolumne smrf_sunang_azimuth_tuolumne

Solar

The exoatmospheric direct solar irradiance at the top of the atmosphere for a wavelength range. The radiation.solar is the Python version and radiation.solar_ipw is the IPW version, both still in SMRF. The solar calculation requires fitting an Akima spline to NASA solar data and integrating between two wavelengths. While the IPW versions of solar performed these spline fittings and integrals, the python solar using scipy.interpolate.Akima1DInterpolator to fit the spline and scipy.integrate.quad to do the integral. Because these are different function, the resulting values will be different.

Below is the difference between IPW solar and python solar, which is less than +/- 0.2 W/m2 over the whole 2016 water year. smrf_solar

Twostream

twostream provides a solution for a single-layer attest_twostreammosphere over horizontal surface. There is a specific test in tests.test_radiation.test_twostream that compares the outputs from the IPW command twostream to the python twostream with identical results. This means that given the same inputs to twostream it replicates the original IPW function. However, because the inputs will be coming from the python version of sunang and solar there will be slight differences. Again, over water year 2016 in Reynolds Mountain East has differences less than 0.2 W/m2.

smrf_twostream

HRRR cloud factor

radiation.model_solar is only used at the moment for HRRR when converting the incoming solar radiation to a cloud factor, the station runs do not use twostream although twostream may be used outside of SMRF to calculate the cloud factor. The net radiation difference in the HRRR test is minimal with most values less that 0.05 W/m2. smrf_hrrr_net_solar_diff

The story is similar in thermal for HRRR, a very small difference of less than 0.016 W/m2. smrf_hrrr_thermal_diff

micahjohnson150 commented 4 years ago

@etrujil and @Hedrick-ARS this PR is on deck for being merged in. If you could take a look at this from a scientific stand point and provide any feedback, I will be going over code in the couple days.

Hedrick-ARS commented 4 years ago

@micahjohnson150 - I'm out of the office today after the science meeting, then at RME Wednesday and Thursday so won't be able to get to this until Friday.

micahjohnson150 commented 4 years ago

@USDA-ARS-NWRC/snow-modelers please take a look when you get a chance. Looking for scientific assessment on this PR.

scotthavens commented 4 years ago

Over the whole upper Tuolumne, the azimuth to the sun can vary up to 1 degree from the west side to the east. Given that the wide of the Tuolumne domain is ~70km, it makes sense that the solar azimuth is slightly different.

etrujil commented 4 years ago

Over the whole upper Tuolumne, the azimuth to the sun can vary up to 1 degree from the west side to the east. Given that the wide of the Tuolumne domain is ~70km, it makes sense that the solar azimuth is slightly different.

@scotthavens I see, I was interpreting as a map of differences between the ipw and python functions. Yes, it all makes sense then. Thanks.

dannymarks commented 4 years ago

Hey Guys, We knew about this way back when the ipw topo and rad programs where being developed; As most of the work done before I met you guys was over relative small areas, I did not worry about it - if the area is small enough - say less than 0.5 degree by 0.5 degree, you can get away with assuming that sun angles are constant over the domain grid; however if the domain is larger (like the whole Sierra) then you have to calculate sun angles, and corresponding azimuths and horizons more carefully. I wrote a series of programs called "global" that uses Lat/Long rather than a conformal projection (Mercator, UTM, etc), and used that for analysis of larger areas (in particular the Columbia River Basin in the Northwest). The IPW utility names all start with a "g":

Topographic Calculations

ggradient - slope and aspect from elevation (global-scale)

ghorizon - horizon for one direction (global-scale)

gshade - cosine of illumination angle (global-scale)

gviewf - sky view and terrain configuration factors (global-scale)

gviewf-mp - multi-processor version of 'gviewf' command

Radiation Modeling: Global-Scale

gelevrad - beam and diffuse irradiance from elevation

gsunlight - sun angles: azimuth and zenith angles of sun's position

gsunweights - Kronrod quadrature weights and sun angles

gtoporad - topographically-corrected solar radiation

gtoporad.24 - topographically-corrected solar radiation for 24-hour period

gtopquad - daily integrated radiation over a topographic grid

If you look at these you should be able to substitute them in as needed...

However, It may not be necessary; for the BRB work (as we had pretty slow computers) we simply adjusted the UTM grid so that the area of interest was in the center of the "zone", so that the sun angle differences from east to west and north to south were not that far off, and the sun angle and azimuth differences were not that off either.

The problem is the Tuo is long and skinny, so the problem occurs even though the actual area is not that large...

Ultimately you will have to do this, but it is much much slower - however, it only has to be done once if you save it. I would do it for the whole domain, or all of California, or the entire western US and be done with it. Or maybe your new code is fast enough that you can deal with it on the fly?

Call me if you want to discuss this...

Danny

On Mon, Mar 16, 2020 at 2:52 PM Scott Havens notifications@github.com wrote:

Over the whole upper Tuolumne, the azimuth to the sun can vary up to 1 degree from the west side to the east. Given that the wide of the Tuolumne domain is ~70km, it makes sense that the solar azimuth is slightly different.

— You are receiving this because you are on a team that was mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#issuecomment-599774339, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4C5MG3H2VLOLUUVIKU3RH2NRPANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

dannymarks commented 4 years ago

Note: when I say much much slower, that is because in UTM you calculate sun angles and azimuths once, and then move on, assuming those are uniform at all points in the domain. However in Lat/Long, you have to do that calculation for every grid point, so for the TUO the calculations take about 1 million times longer...

On Mon, Mar 16, 2020 at 3:46 PM Danny Marks ars.danny@gmail.com wrote:

Hey Guys, We knew about this way back when the ipw topo and rad programs where being developed; As most of the work done before I met you guys was over relative small areas, I did not worry about it - if the area is small enough - say less than 0.5 degree by 0.5 degree, you can get away with assuming that sun angles are constant over the domain grid; however if the domain is larger (like the whole Sierra) then you have to calculate sun angles, and corresponding azimuths and horizons more carefully. I wrote a series of programs called "global" that uses Lat/Long rather than a conformal projection (Mercator, UTM, etc), and used that for analysis of larger areas (in particular the Columbia River Basin in the Northwest). The IPW utility names all start with a "g":

Topographic Calculations

ggradient - slope and aspect from elevation (global-scale)

ghorizon - horizon for one direction (global-scale)

gshade - cosine of illumination angle (global-scale)

gviewf - sky view and terrain configuration factors (global-scale)

gviewf-mp - multi-processor version of 'gviewf' command

Radiation Modeling: Global-Scale

gelevrad - beam and diffuse irradiance from elevation

gsunlight - sun angles: azimuth and zenith angles of sun's position

gsunweights - Kronrod quadrature weights and sun angles

gtoporad - topographically-corrected solar radiation

gtoporad.24 - topographically-corrected solar radiation for 24-hour period

gtopquad - daily integrated radiation over a topographic grid

If you look at these you should be able to substitute them in as needed...

However, It may not be necessary; for the BRB work (as we had pretty slow computers) we simply adjusted the UTM grid so that the area of interest was in the center of the "zone", so that the sun angle differences from east to west and north to south were not that far off, and the sun angle and azimuth differences were not that off either.

The problem is the Tuo is long and skinny, so the problem occurs even though the actual area is not that large...

Ultimately you will have to do this, but it is much much slower - however, it only has to be done once if you save it. I would do it for the whole domain, or all of California, or the entire western US and be done with it. Or maybe your new code is fast enough that you can deal with it on the fly?

Call me if you want to discuss this...

Danny

On Mon, Mar 16, 2020 at 2:52 PM Scott Havens notifications@github.com wrote:

Over the whole upper Tuolumne, the azimuth to the sun can vary up to 1 degree from the west side to the east. Given that the wide of the Tuolumne domain is ~70km, it makes sense that the solar azimuth is slightly different.

— You are receiving this because you are on a team that was mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#issuecomment-599774339, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4C5MG3H2VLOLUUVIKU3RH2NRPANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

scotthavens commented 4 years ago

Thanks for the input @dannymarks moving to lat/long may be a good way to go in the future or have another option. A great example is the Tuolumne where it overlaps UTM zones so having that as lat/long would be much more accurate. I've opened up an issue #141 so we can bring in the lat/long as an option for users in larger basins.

dannymarks commented 4 years ago

Good job guys! Let’s hope 1) it works, & 2) it isn’t incredibly slow... Danny

On Wed, Mar 18, 2020 at 10:36 AM Scott Havens notifications@github.com wrote:

@scotthavens commented on this pull request.

In smrf/envphys/sunang.py https://github.com/USDA-ARS-NWRC/smrf/pull/123#discussion_r394524975:

+def sunang(date_time, latitude, longitude, truncate=True):

  • """
  • Calculate the sun angle (the azimuth and zenith angles of
  • the sun's position) for a given geodetic location for a single
  • date time and coordinates. The function can take either latitude
  • longitude position as a single point or numpy array.
  • Args:
  • date_time: python datetime object
  • latitude: value or np.ndarray (in degrees)
  • longitude: value or np.ndarray (in degrees)
  • truncate: True will replicate the IPW output precision, not applied
  • if position is an array
  • Returns
  • cosz - cosine of the zeinith angle, same shape as input position

Fixed

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#discussion_r394524975, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4CYLPXVG5V4JNSISJT3RIEBC3ANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

micahjohnson150 commented 4 years ago

@robertson-mark we're still looking for your review on this.

dannymarks commented 4 years ago

If its my code, the the correspondence to python is coincidence. One thing, I believe the M_PI is a C code defined value... Not sure how its defined, but I will try to track it down...

Danny

On Wed, Mar 18, 2020 at 2:02 PM Micah Johnson notifications@github.com wrote:

@robertson-mark https://github.com/robertson-mark we're still looking for your review on this.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#issuecomment-600856631, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4C3J6WZ32A4GFAPJFU3RIEZG5ANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

scotthavens commented 4 years ago

@micahjohnson150 Does Python 3 still have that float times integer problem? I thought they did away with that, but changing it doesn't hurt.

@dannymarks The M_PI is straight from the C code and is defined somewhere in SMRF. It's just a convenience thing that could be changed in favor of Python built-ins but I don't see it as a big hindrance to the code.

micahjohnson150 commented 4 years ago

Nope you're correct. It was a clarity thing thats all.

dannymarks commented 4 years ago

Well I’m impressed, Micah Let know how it compares Danny

On Wed, Mar 18, 2020 at 2:21 PM Micah Johnson notifications@github.com wrote:

Nope you're correct. It was a clarity thing thats all.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#issuecomment-600865125, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4C24T7WLQVVGVVRUMRTRIE3PDANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon

micahjohnson150 commented 4 years ago

@dannymarks I would love to take credit for all this work here but this completely Scott's work and I am just providing my review which is a part of our new code review protocol.

dannymarks commented 4 years ago

I know, but it is a team effort - success will require everyone.

On Wed, Mar 18, 2020 at 3:35 PM Mark Robertson notifications@github.com wrote:

@robertson-mark approved this pull request.

This is the first time I have looked at much of this code, but from what I understand it looks good to me.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/USDA-ARS-NWRC/smrf/pull/123#pullrequestreview-377284526, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHUC4C42YOEUGI5OMGSL67LRIFEC5ANCNFSM4I7AZPLA .

-- Dr. Danny Marks Retired Snow Scientist Bend, Oregon