Daylight Savings Time and UTC Offset for s.geometry.from_time_and_location

[Luke Brown]

Dear Robin,

attached is my code. I call s.geometry.from_time_and_location in the method called "setup".

My issue is that I have imagery flown during UTC minus 4 due to DST, but the most accurate match between measured radiance and 6S output is for UTC minus 5. May we discuss my code's use of the method, and also clarify my understanding of the timestamp in Py6S? 

All my best,

Luke Brown

[Robin Wilson]

Hi Luke,

Hmm, that's a good question - I haven't been very clear about how the from_time_and_location method deals with timezones.

What times are you passing the method at the moment? When you say you're passing times with UTC-4 and UTC-5, are you passing the times in UTC or in the local timezone?

Internally, Py6S uses the Pysolar library (https://pysolar.readthedocs.io/en/latest/) to calculate solar angles for a given time and location. For various reasons we're using an older version of Pysolar (v0.6), which requires that the time is given in UTC. It doesn't matter whether the datetime object that you pass to from_time_and_location (which is then passed to Pysolar.GetAltitude) has timezone information attached to it. If it doesn't, it will be assumed to be UTC. However, if it does have timezone information attached then this will be ignored and it will still be assumed to be UTC. Basically, it's treated as UTC no matter what.

So, hopefully your use case should work if you convert all times to UTC before passing to from_time_and_location.

I'm aware that the documentation for this needs improving, and Py6S should also give a warning if a timestamp is passed with a non-UTC timezone - I will add this to my todo list. In fact, Py6S could be updated to use a later version of Pysolar which has some of this checking built-in - the reason for not doing this earlier was that later versions of Pysolar don't support Python 2, and we still had some users on Python 2. It is probably reasonable now to only support Python 3, so I will add this to my todo list too.

Hope that helps,

Robin

--
You received this message because you are subscribed to the Google Groups "Py6S" group.
To unsubscribe from this group and stop receiving emails from it, send an email to py6s+uns...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/ab922487-0cc1-413e-8de5-57da61c63950n%40googlegroups.com.

[Luke Brown]

Robin, 

First, thank you very much for clarifying the internals of how the program is handling timestamps; it is good information to have. 

However, my issue persists and hopefully this concrete example below can help.

The 6S values are much larger than what were measured by our CASI 1500H (VNIR, 48 bands):

the code to produce this:

______________________________________________________________ 

 

from Py6S import *
import numpy as np
import rasterio as rio
# import pandas as pd
# from datetime import datetime, timedelta
import matplotlib.pyplot as plt


# Setup 6S model
img6s = SixS()

img6s.altitudes = Altitudes()
img6s.altitudes.set_target_custom_altitude(0.309)
img6s.altitudes.set_sensor_custom_altitude(1.289, aot=0.1, ozone=0.9, water=2.0)

img6s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.45)

img6s.atmos_profile = AtmosProfile.FromLatitudeAndDate(40.665442, '2020-07-18')

img6s.aero_profile = AeroProfile.Continental

acq_datetime_string = "07/18/2020 15:48:30"

img6s.geometry = Geometry.User()
img6s.geometry.from_time_and_location(40.665442, -77.946006, '2020-07-18 15:48:30+00:00', 357.040699, 5.474)  # here I am using the UTC timestamp, which because of DST is ~11:46 + 4 hours

wv, results = SixSHelpers.Wavelengths.run_wavelengths(img6s, np.linspace(0.38238, 1.053674, 48), output_name='apparent_radiance')


# Getting our measured tarp spectrum:
image = r'D:\CASI_2020_07_18_114614_rad.bsq'  # Filename contains local time of flight (11:46 am)
x, y = 626, 8086 # image geometry pixel coordinates of white calibration tarp with 45% reflectance 

 
with rio.open(image, 'r') as src:
    for val in src.sample([(x, y)]):
        tarp_spectra = val / 1000 # extracting tarp spectrum from image and converting mW to W

# Plotting estimated and measured radiance together
plt.plot(wv, results, label='6S ' + acq_datetime_string)
plt.plot(wv, tarp_spectra, label='measured tarp radiance')
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Radiance (W/m^2/sr/nm)")
plt.legend()
plt.show()

______________________________________________________________

Please let me know if any of this is unclear. I have not been able to figure out the mismatch between measured and estimated radiance. I hope it's something I'm doing incorrectly!

Again thank you,

Luke Brown

Technical Specialist - Spectral analytics

Quantum Spatial, Inc.

[Luke Brown]

[Robin Wilson]

Hi Luke,

I haven't been able to spend long on this tonight - but I have an initial thought which you could investigate further. The units you're using on your y axis are W/m^2/sr/nm, but 6S provides output as W/m^2/sr/um (ie. per micrometre rather than per nanometre). To convert, you'll need to divide each data point by the bandwidth in nanometres.

From your code, np.diff(np.linspace(0.38238, 1.053674, 48)) gives an array where all values are 0.01428285 (ie. 14.2nm). Dividing all radiance values by 14.2 brings them far closer to the measured tarp radiance (though I can't work out exactly how close, as I don't have the tarp data).

Hope that helps,

Robin

To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/b635183d-2c4e-4546-9ebf-4e9b1706ea87n%40googlegroups.com.

[Luke Brown]

Robin, 

Thank you for your suggestion. I wish I had recognized this when looking over the other responses in the group as it appears to be a common oversight. The magnitudes are now in the same neighborhood, but my RMSE is still ~ 3 units. This is workable for my purposes, but I am interested to know if you have other suggestions for improvement. Would it be helpful to have any other information?

Luke