import requests

request_url = (
    "https://modis.ornl.gov/rst/api/v1/MOD13Q1/subsetOrder?"
    "latitude={lat}&"
    "longitude={lon}&"
    "email=jjmcnelis%40outlook.com&"
    "uid={uid}&"
    "startDate={stdt}&"
    "endDate={enddt}&"
    "kmAboveBelow={ab}&"
    "kmLeftRight={lr}"
).format(
    lat   =  28.27149,
    lon   = -81.9668,
    uid   = "beta",
    stdt  = "A2012001",
    enddt = "A2012111",
    ab    =  0,
    lr    =  0
)

print("Request URL:\n"+request_url)

response = requests.get(request_url)

if response.status_code==200:
    orderid = response.text
    print("Order ID:"); print(orderid)
else:
    print(response.text)

import numpy as np
import pandas as pd
import xarray as xr
from io import StringIO
from datetime import datetime
from datetime import timedelta
from pyproj import Proj, transform

order_nir = {
    "nir_black": None, 
    "nir_white": None, 
    "nir_actual": None}


def mdate2dt(modis_date):
    """AYYYYDOY::A2019001 >> yyyy-mm-dd"""
    idate = datetime(int(modis_date[1:5]),1,1)
    idelt = timedelta(int(modis_date[5:])-1)    
    return((idate+idelt).strftime('%Y-%m-%d'))

colnames = ["id", "prod", "mdate", "loc", "pdate", "var", "data"]
for csv in list(order_nir.keys()):               # iterate over list of nir
    
    csv_url = (
        "https://modis.ornl.gov/subsetdata/"     # global tool orders dir
        "03Jul2019_13:43:49_346167250"           # order id
        "L28.27149"                              # latitude of center
        "L-81.9668"                              # longitude of center
        "S1L1_MCD43A/"+csv+".csv")               
    r = requests.get(csv_url)                    # download with requests.get
    
    order_nir[csv] = pd.read_csv(                # parse response text to pd 
        StringIO(r.text),                        #  StringIO pseudofile
        names=colnames,
        header=None, 
        sep=",") 
    
    # set the data column to mask the fill values (32767.0)
    order_nir[csv].loc[order_nir[csv].data==32767.0, "data"] = np.nan
    
    # get dates from the third column (base zero); set to index
    order_nir[csv].index = [mdate2dt(md) for md in order_nir[csv]["mdate"]] 

    
order_nir["nir_black"].head(5)               # print five rows of nir black

import glob

bsa_result = glob.glob("result/black_albedo*.nc")
wsa_result = glob.glob("result/white_albedo*.nc")
alb_result = glob.glob("result/blue_albedo*.nc")

print("Black sky albedo outputs for 2018:"); bsa_result

bsa_result = sorted(bsa_result, key=lambda x: int(x.split("_")[-1][:-3]))
wsa_result = sorted(wsa_result, key=lambda x: int(x.split("_")[-1][:-3]))
alb_result = sorted(alb_result, key=lambda x: int(x.split("_")[-1][:-3]))

bsa_result

ll = order_nir["nir_black"].iloc[0]["loc"].split("Samp")[0][3:].split("Lon")
ll

# EPSG code for geographic WGS84 is our input projection this time
inproj = Proj(init="epsg:4326")

# build the MODIS sinusoidal proj4 string again
getpar = lambda a: str(bsa.crs.attrs[a])
outproj = Proj(" ".join([
    "+proj=sinu",
    "+lon_0="+getpar("longitude_of_central_meridian"),
    "+x_0="+getpar("false_easting"),
    "+y_0="+getpar("false_northing"),
    "+a="+getpar("semi_major_axis"),
    "+b="+getpar("semi_minor_axis"),
    "+units="+"meter +no_defs"]))

# pass the two projection definitions to transform along with the lat,lon
x,y = transform(inproj, outproj, float(ll[1]), float(ll[0]))

print("Sinusoidal equivalent for the lat, lon from the DAAC subset:"); x,y

import numpy as np

find_nearest = lambda array, value: array[(np.abs(array - value)).argmin()]
x_near = find_nearest(bsa.x, x).data.item()
y_near = find_nearest(bsa.y, y).data.item()

print("The nearest x,y coordinate in our dataset:"); x_near, y_near

our_nir_bsa  = bsa["BRDF_Albedo_Parameters_nir"]  # get nir BSA for our ds
ornl_nir_bsa = order_nir["nir_black"]["data"]     # get nir BSA for subset

# select only the numpy arrays
sza = bsa["solar_zenith_angle"].sel(x=x_near, y=y_near).data.compute()
ours  = our_nir_bsa.sel(x=x_near, y=y_near).data.compute()
ornls = ornl_nir_bsa.values
diff  = ours-ornls

# make a df
df = pd.DataFrame({
    "ORNLDAAC": ornls, 
    "OURWORKFLOW": ours, 
    "DIFFERENCE": diff,
    "SZA": sza
}, index=bsa.time.dt.dayofyear)

df.head(5)

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14 

fig, axs = plt.subplots(3, 1, sharex=True, figsize=(14,8))
fig.suptitle("Computed black sky albedo (MCD43A1): comparison "
             "ORNL DAAC versus Workflow")

df["ORNLDAAC"].plot(ax=axs[0], color="black")
axs[0].set_ylabel("ORNL DAAC")
axs[0].grid(True)

df["OURWORKFLOW"].plot(ax=axs[1], color="black")
axs[1].set_ylabel("OUR WORKFLOW")
axs[1].grid(True)

df["DIFFERENCE"].plot(ax=axs[2])
axs[2].set_ylabel("DIFFERENCE")
axs[2].grid(True)

axs3 = axs[2].twinx()
df["SZA"].plot(ax=axs3, color="orange")
axs3.set_ylabel("ZENITH ANGLE")

plt.show()

# we will plot the NIR band
band_name = "BRDF_Albedo_Parameters_nir"
xyd = dict(x=x_near, y=y_near)

# get pixel from black sky albedo results
with xr.open_mfdataset(bsa_result) as bsa:
    sza = bsa["solar_zenith_angle"].sel(**xyd).data.compute()
    ours_nir_bsa  = bsa[band_name].sel(**xyd).data.compute()
    ornls_nir_bsa  = order_nir["nir_black"]["data"].values
    
# get pixel from white sky albedo results
with xr.open_mfdataset(wsa_result) as wsa:
    ours_nir_wsa  = wsa[band_name].sel(**xyd).data.compute()

# get pixel from blue sky albedo results
with xr.open_mfdataset(alb_result) as alb: 
    ours_nir_alb  = alb[band_name].sel(**xyd).data.compute()

# select only the numpy arrays
ornls = ornl_nir_bsa.values
diff  = ours-ornls

# make a df
df = pd.DataFrame({
    "SZA": sza,
    "ORNLDAAC_BLACK": ornls_nir_bsa, 
    "WORKFLOW_BLACK": ours_nir_bsa, 
    "WORKFLOW_WHITE": ours_nir_wsa, 
    "WORKFLOW_BLUE": ours_nir_alb, 
    "DIFFERENCE": ours_nir_bsa-ornls_nir_bsa
}, index=bsa.time.dt.dayofyear)

# and display the first five rows
df.head(5)

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14  
plt.style.use('seaborn-muted')

fig = plt.figure(num=1, figsize=(14, 4))

lnstyle = dict(color="black", alpha=0.75,)
df["WORKFLOW_BLUE"].plot(label="blue", linewidth=2.5)
df["WORKFLOW_BLACK"].plot(linestyle='dashed', label="black", **lnstyle)
df["WORKFLOW_WHITE"].plot(linestyle='dotted', label="white", **lnstyle)
location = "[ %s, %s ]" % (ll[0], ll[1])

plt.grid(alpha=0.5)
plt.legend(loc="upper right")
plt.xlabel("DAY OF THE YEAR")
plt.ylabel("ALBEDO [ 0.0 - 1.0 ]")
plt.title("Blue, black, white sky albedo for 2018 at "+location)

plt.show()

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14  
plt.style.use('seaborn-muted')

fig = plt.figure(num=2, figsize=(14, 4))

df["DIFFERENCE"].plot(label="difference")
plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
plt.xlabel("DAY OF THE YEAR", labelpad=10)
plt.ylim(-0.000001, 0.000005)
plt.ylabel("DIFFERENCE")
ax = plt.gca()

ax2 = ax.twinx()
plt.ylim(0, 55)
ax2.fill_between(df.index, df["SZA"], color="gray", alpha=0.1)
ax2.set_ylabel("ZENITH ANGLE [ DEG ]", labelpad=10)
location = "[ %s, %s ]" % (ll[0], ll[1])

plt.grid(alpha=0.5)
plt.title("Albedo difference [ WORKFLOW - ORNLDAAC ] at "+location)

plt.show()