Viewtime and conversion to UTC#
This notebook demonstrates the process of:
Accessing the viewtime asset from the coarse resolution fusion input
Converting the viewtime in local solar time to UTC
We’ll start with the necessary inputs.
import json
import requests
import pystac
from pystac_client import Client
from shapely.geometry import Point
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import base64
from FH_Hydrosat import FH_Hydrosat
from datetime import datetime as DT
import datetime as DT
The next few cells will execute the STAC search as usual. In this case, we’ll get the MOD21A1D coarse resolution fusion input as well as the HDST products.
with open('creds.json') as f:
creds = json.loads(f.read())
userpass = f"{creds['username']}:{creds['password']}"
b64 = base64.b64encode(userpass.encode()).decode()
headers = {'Authorization':'Basic ' + b64}
cat_url = 'https://fusion-stac.hydrosat.com'
catalog = Client.open(cat_url, headers)
geom = {'type': 'Point', 'coordinates': [-120.21102905273436,36.535019209789]} # Point for Hanford, CA, USA
start_date = "2021-08-17T00:00:00Z"
end_date = "2021-08-30T00:00:00Z"
search = catalog.search(
collections = ["prepped_inputs_mod21a1d"], #The coarse resolution MOD21A1D LST input to fusion
intersects = geom,
datetime = [start_date, end_date],
max_items = 500
)
# mod_items = list(search.items()) # for pystac-client >= 0.4.0
mod_items = list(search.get_all_items()) # for pystac-client < 0.4.0
mod_items.reverse() # make the results ascending in time
search = catalog.search(
collections = ["starfm_predictions_modis_landsat", "pydms_sharpened_landsat"], # The HDST products
intersects = geom,
datetime = [start_date, end_date],
max_items = 500
)
# hdst_items = list(search.items()) # for pystac-client >= 0.4.0
hdst_items = list(search.get_all_items()) # for pystac-client < 0.4.0
hdst_items.reverse() # make the results ascending in time
print(len(mod_items))
print(len(hdst_items))
14
13
mgrs_tiles = []
for i in hdst_items:
for l in i.to_dict()['links']:
if 'element84' in l['href']:
mgrs_tiles.append(l['href'].split(r'/')[-1].split('_')[1])
hdst_items = [i for i in hdst_items if mgrs_tiles[0] in i.id]
Next we’ll use the FH_Hydrosat
helper functions to create a list of the “viewtimes” in the timeseries for a given point geometry. These times come from the lst_view_time
asset for MOD21A1D items.
pt_wgs84 = Point(geom['coordinates'][0], geom['coordinates'][1])
viewtime_res = FH_Hydrosat(mod_items, asset='lst_view_time')
viewtime_ts = viewtime_res.point_time_series_from_items(pt_wgs84, tol=2000, nproc=6)
viewtime_ts = [v*0.1 for v in viewtime_ts]
using 6 processes to sample 14 assets
viewtime_ts
[10.4,
11.100000000000001,
10.200000000000001,
10.9,
10.0,
10.700000000000001,
11.4,
10.5,
11.200000000000001,
10.3,
11.0,
10.100000000000001,
10.8,
11.5]
Now we’ll plot the viewtimes. These are in local solar time. You’ll notice that they vary from day to day, but are generally within 1.5 hours of eachother. The viewtime can also vary across a single MOD21A1D product. For more information, see the MOD21A1D ATBD: https://lpdaac.usgs.gov/documents/1399/MOD21_ATBD.pdf
plt.plot(viewtime_ts)
plt.show()
![../_images/8ca6ac052ea06955701afae80187c30a1b122c5fb817236a8aafd2fbc8804e96.png](../_images/8ca6ac052ea06955701afae80187c30a1b122c5fb817236a8aafd2fbc8804e96.png)
Next we’ll create a timeseries of LST from the HDST product for the same point.
lst_res = FH_Hydrosat(hdst_items, asset='lst')
lst_ts = lst_res.point_time_series_from_items(pt_wgs84, tol=2000, nproc=6)
using 6 processes to sample 13 assets
You’ll notice that some of the HDST products only have a date and do not show a nominal observation time.
hdst_dt = lst_res.datetime
hdst_df = pd.DataFrame({'fused_lst': lst_ts,
'datetime': pd.to_datetime(hdst_dt)}).sort_values(by='datetime')
hdst_df['fused_lst'] = hdst_df['fused_lst'].astype(float)
hdst_df.dropna(inplace=True)
hdst_df
fused_lst | datetime | |
---|---|---|
0 | 319.647186 | 2021-08-17 00:00:00+00:00 |
1 | 311.785156 | 2021-08-18 00:00:00+00:00 |
2 | 301.110474 | 2021-08-19 00:00:00+00:00 |
3 | 309.528900 | 2021-08-20 00:00:00+00:00 |
4 | 309.813385 | 2021-08-21 18:40:04.472248+00:00 |
5 | 317.071136 | 2021-08-22 00:00:00+00:00 |
6 | 318.495239 | 2021-08-23 00:00:00+00:00 |
7 | 317.482910 | 2021-08-24 00:00:00+00:00 |
8 | 316.709290 | 2021-08-25 00:00:00+00:00 |
9 | 317.916840 | 2021-08-26 00:00:00+00:00 |
10 | 321.382233 | 2021-08-27 00:00:00+00:00 |
11 | 319.674347 | 2021-08-28 00:00:00+00:00 |
12 | 322.872528 | 2021-08-29 00:00:00+00:00 |
We can derive the nominal observation time by using the viewtime (local solar time) from the coarse resolution assets above and converting it to UTC.
# adjust the datetime
viewtime_ts = [x for x in viewtime_ts if str(x) != 'nan']
local_solar_time = []
data_type = []
for x,dt in zip(viewtime_ts, hdst_df['datetime']):
if '00:00:00' in str(dt):
local_solar_time.append(dt + DT.timedelta(seconds=x*60*60.0))
data_type.append('fused')
else:
local_solar_time.append(dt)
data_type.append('landsat')
#local_solar_time = [dt + DT.timedelta(seconds=x*60*60.0) if '00:00:00' in str(dt) else dt for x,dt in zip(viewtime_ts, hdst_df['datetime'])]
# adjust for longitude and UTC if not associated with landsat
pt_lon = geom['coordinates'][0]
pt_utc_offset = -pt_lon/15
pt_utc = []
for dt, _type in zip(local_solar_time, data_type):
if _type == 'landsat':
pt_utc.append(dt)
else:
pt_utc.append(dt + DT.timedelta(seconds=pt_utc_offset*60*60.0))
#pt_utc = [dt + DT.timedelta(seconds=pt_utc_offset*60*60.0) for dt in local_solar_time]
# pv equations for site
gmt_offset = 8
d = np.array([d.timetuple().tm_yday for d in hdst_df['datetime']])
B = (360/365)*(d-81)
EoT = 9.87*np.sin(np.deg2rad(2*B)) - 7.53*np.cos(np.deg2rad(B)) - 1.5*np.sin(np.deg2rad(B))
TC = 4*(pt_lon + gmt_offset*15) + EoT
LT = []
for _dtz, _tc, _type in zip(local_solar_time, TC, data_type):
if _type == 'landsat':
LT.append(_dtz)
else:
LT.append(_dtz - DT.timedelta(hours=_tc/60) + DT.timedelta(hours=gmt_offset))
# LT = [_dtz - DT.timedelta(hours=_tc/60) + DT.timedelta(hours=gmt_offset) for _dtz, _tc in zip(hdst_df['datetime'], TC)]
hdst_df['UTC_time'] = LT
hdst_df
fused_lst | datetime | UTC_time | |
---|---|---|---|
0 | 319.647186 | 2021-08-17 00:00:00+00:00 | 2021-08-17 18:28:35.861866+00:00 |
1 | 311.785156 | 2021-08-18 00:00:00+00:00 | 2021-08-18 19:10:22.330436+00:00 |
2 | 301.110474 | 2021-08-19 00:00:00+00:00 | 2021-08-19 18:16:08.255089+00:00 |
3 | 309.528900 | 2021-08-20 00:00:00+00:00 | 2021-08-20 18:57:53.647627+00:00 |
4 | 309.813385 | 2021-08-21 18:40:04.472248+00:00 | 2021-08-21 18:40:04.472248+00:00 |
5 | 317.071136 | 2021-08-22 00:00:00+00:00 | 2021-08-22 18:45:22.887148+00:00 |
6 | 318.495239 | 2021-08-23 00:00:00+00:00 | 2021-08-23 19:27:06.761276+00:00 |
7 | 317.482910 | 2021-08-24 00:00:00+00:00 | 2021-08-24 18:32:50.157547+00:00 |
8 | 316.709290 | 2021-08-25 00:00:00+00:00 | 2021-08-25 19:14:33.091216+00:00 |
9 | 317.916840 | 2021-08-26 00:00:00+00:00 | 2021-08-26 18:20:15.578181+00:00 |
10 | 321.382233 | 2021-08-27 00:00:00+00:00 | 2021-08-27 19:01:57.634964+00:00 |
11 | 319.674347 | 2021-08-28 00:00:00+00:00 | 2021-08-28 18:07:39.278696+00:00 |
12 | 322.872528 | 2021-08-29 00:00:00+00:00 | 2021-08-29 18:49:20.527096+00:00 |