%matplotlib notebook
import numpy as np
from dh_py_access import package_api
import dh_py_access.lib.datahub as datahub
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.basemap import Basemap
from matplotlib import colors
import warnings
warnings.filterwarnings("ignore")

"""
Explanation: What We Saw from Satellite-Retrieved Snow Cover Data
In this demo we are going to show snow cover data for the Christmas Eve using MODIS/Terra Snow Cover Daily L3 Global, Version 6 dataset. 
This dataset reports the percentage of snow-covered land in 0.05° (approx. 5 km) resolution Climate Modeling Grid (CMG) cells. Percentages are computed from snow cover observations in the MODIS/Terra Snow Cover Daily L3 Global 500m Grid (MOD10A1) dataset. Cloud cover percentages are also provided. 
We are going to use daily snow extent variable and we will be looking at the Northern part of the globe, from latitude 10N to 71N. 
MODIS  (Moderate Resolution Imaging Spectroradiometer) is a key instrument aboard the Terra and Aqua satellites. Terra's orbit around the Earth is timed so that it passes from north to south across the equator in the morning, while Aqua passes south to north over the equator in the afternoon. Terra MODIS and Aqua MODIS are viewing the entire Earth's surface every 1 to 2 days, acquiring data in 36 spectral bands. 
As MODIS uses visible light from Earth, data is not collected during the night. Also, clouds are influencing the results. In this dataset there are some values that are referring to different kind of conditions. We have created a map to illustrate why there is no data collected from some of these places.
In this demo we will:
1) use the Planet OS package API to fetch data;
2) show snow cover map from Christmas 2017;
3) show map with masked value meanings.
End of explanation
"""


server = 'api.planetos.com'
API_key = open('APIKEY').readlines()[0].strip() #'<YOUR API KEY HERE>'
version = 'v1'

"""
Explanation: <font color='red'>Please put your datahub API key into a file called APIKEY and place it to the notebook folder or assign your API key directly to the variable API_key!</font>
End of explanation
"""


dh=datahub.datahub(server,version,API_key)
dataset='nsidc_modis_snow_cover_global'
variable_name1 = 'Day_CMG_Snow_Cover'
time_start = '2017-12-24T00:00:00'
time_end = '2017-12-24T23:00:00'

"""
Explanation: At first, we need to define the dataset name and a variable we want to use.
End of explanation
"""


area_name = 'north'
latitude_north = 71; longitude_west = -160
latitude_south = 10; longitude_east = 180

"""
Explanation: Then we define spatial range. At this time, we will look into Northern Hemisphere from 10N to 71N.
End of explanation
"""


package = package_api.package_api(dh,dataset,variable_name1,longitude_west,longitude_east,latitude_south,latitude_north,time_start,time_end,area_name=area_name)

package.make_package()

package.download_package()

"""
Explanation: Download the data with package API

Create package objects
Send commands for the package creation
Download the package files
End of explanation
"""


dd1 = xr.open_dataset(package.local_file_name)
dd1['masked_snow_cover'] = dd1.Day_CMG_Snow_Cover

"""
Explanation: Work with downloaded files
We start by opening the files with xarray and by making the data column names masked_snow_cover. At first, we will have the same data as Day_CMG_Snow_Cover, but later, we will mask out values that are explained under the Day_CMG_Snow_Cover Key attribute.
End of explanation
"""


masks = dd1.Day_CMG_Snow_Cover.Key.split(',')
print (masks)
masks = masks[1:]

for mask in masks:
    mask_value = float(mask[:4])
    dd1['masked_snow_cover'].data = np.ma.masked_where(dd1['masked_snow_cover'].data == mask_value,dd1['masked_snow_cover'].data)

"""
Explanation: Here we are finding the mask values and we will print them out. However, we will ingnore the data values.
End of explanation
"""


m_demo=Basemap(projection='merc',llcrnrlat=latitude_south,urcrnrlat=latitude_north,\
            llcrnrlon=-160,urcrnrlon=longitude_east,lat_ts=(latitude_north+latitude_south)/2,resolution='l')
Dx,Dy = np.meshgrid(dd1.lon.data,dd1.lat.data)
dx,dy = m_demo(Dx,Dy)

"""
Explanation: Now we are making Basemap of the area on what we will show the data.
End of explanation
"""


fig=plt.figure(figsize=(10,5))
mp = m_demo.pcolormesh(dx,dy,dd1['masked_snow_cover'].data[0,:,:],vmin=0,vmax=100,cmap = 'Blues')
m_demo.colorbar(mp,label = dd1.Day_CMG_Snow_Cover.long_name + '[' + dd1.Day_CMG_Snow_Cover.units + ']')
plt.title('Snow Cover')
m_demo.drawcoastlines()
m_demo.drawcountries()

"""
Explanation: Then it is time to plot all the data,so that we can see where the snow occured during December 24th. However, we will  not able to see the data from some of the locations in the North. Reasons to why are explained below.
On this map, we can see that Canada had plenty of snow during Christmas, also some places in the North of US, mostly in the mountains. While in Europe we didn't see much snow. Although, people in the Alps had a white Christmas, as did people in some parts of Norway, Sweden, Finland and Estonia. The eastern part of the world, especially Russia also had snow.
If you want to have closer look, you can use tools below the map to zoom in. We definitely encourage you to try that.
End of explanation
"""


masked_values = dd1.Day_CMG_Snow_Cover.data[0,:,:]
masked_values = np.ma.masked_where(masked_values < 101,masked_values)
masked_values[np.isnan(masked_values)] = 255

"""
Explanation: As we briefly explained above we couldn't show data for some of the places. We will now create a map that highlights these places.For doing that we need to get rid off the real values and also make fill values back to 255 as they occur as the NaN at the moment.
End of explanation
"""


#107=lake ice', ' 111=night', ' 237=inland water', ' 239=ocean', ' 250=cloud obscured water', ' 253=data not mapped', ' 255=fill'
fig=plt.figure(figsize=(10,2))
cmap = colors.ListedColormap(['grey','black','pink','cyan','olive','red','blue'])
bounds = [106,108,112,238,240,251,254,256]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
d = m_demo.pcolormesh(dx,dy,masked_values,cmap=cmap, norm=norm)
cbar = m_demo.colorbar(d,ticks=[107,110,175,239,245.5,252.5,255])
cbar.ax.tick_params(labelsize=7) 
cbar.ax.set_yticklabels(masks) 
m_demo.drawcoastlines()
m_demo.drawcountries()
plt.show()

"""
Explanation: On the new map below we can see that in some of the northern parts there is no data because MODIS instrument requires sunlight to operate (it's like a digital camera in many respects). During Christmas most of that region has a polar night or very little daylight. Therefore, the data about snow can't be shown. In some places between 10N to 29N latitude, the data wasn't mapped, because those areas were out of the sensor range .
End of explanation
"""


