#!/usr/bin/env python
# coding: utf-8

# <a name="top"></a>
# <div style="width:1000 px">
# 
# <div style="float:right; width:98 px; height:98px;">
# <img src="https://raw.githubusercontent.com/Unidata/MetPy/master/src/metpy/plots/_static/unidata_150x150.png" alt="Unidata Logo" style="height: 98px;">
# </div>
# 
# <h1>python-awips: Working with Satellite Data</h1>
# <h3>Unidata AMS 2021 Student Conference</h3>
# 
# <div style="clear:both"></div>
# </div>
# 
# ---
# 
# <div style="float:right; width:250 px"><img src="../../instructors/images/python-awips-satellite-preview.png" alt="Colorized GOES East Mesoscale image" style="height: 300px;"></div>
# 
# ### Focuses
# 
# * Investigate what satellite data is available from an EDEX server.
# * Define map properties in a function that can be used to plot multiple images.
# * Retreive Mesoscale GOES satellite grid data from an EDEX server.
# * Use matplotlib pcolormesh to plot the colorized images with a colorbar.
# 
# 
# ### Objectives
# 
# 1. [Define Data Request](#1.-Define-Data-Request)
# 1. [View Optional Identifiers](#2.-View-Optional-Identifiers)
# 1. [View Sources](#3.-View-Sources)
# 1. [View Creating Entities](#4.-View-Creating-Entities)
# 1. [View Sector IDs](#5.-View-Sector-IDs)
# 1. [Create a Satellite Product Tree](#6.-Create-a-Satellite-Product-Tree)
# 1. [Define Map Properties](#7.-Define-Map-Properties)
# 1. [Plot Image Data!](#8.-Plot-Image-Data!)
# 
# ---

# ### Imports

# In[ ]:


from awips.dataaccess import DataAccessLayer
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.pyplot as plt
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import numpy as np
import datetime


# ---
# 
# ## 1. Define Data Request
# 
# If you read through the [python-awips: How to Access Data](https://nbviewer.jupyter.org/github/Unidata/pyaos-ams-2021/blob/master/notebooks/dataAccess/python-awips-HowToAccessData.ipynb) training, you will know that we need to set an EDEX url to access our server, and then we create a data request.  In this example we use *satellite* as the data type to define our request.

# In[ ]:


# Create an EDEX data request
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("satellite")

# Take a look at our request
print(request)


# <a href="#top">Top</a>
# 
# ---
# 
# ## 2. View Optional Identifiers
# 
# There can be many sources of satellite data stored in an EDEX server.  First, let's take a look at the available "identifiers".  We will then look a little more in depth at a few of the identifiers.

# In[ ]:


# get optional identifiers for satellite datatype
identifiers = set(DataAccessLayer.getOptionalIdentifiers(request))

print("Available Identifiers:")
for id in identifiers:
    if id.lower() == 'datauri':
        continue
    print(" - " + id)


# <a href="#top">Top</a>
# 
# ---
# ## 3. View Sources
# 
# Here, use *source* as the identifier, and take a look at the available sources are.  The *source* is a parameter that is set in the netcdf file, and gives some insight to where the product is coming from.

# In[ ]:


# Show available sources
identifier = "source"
sources = DataAccessLayer.getIdentifierValues(request, identifier)
print(identifier + ":")
print(list(sources))


# <a href="#top">Top</a>
# 
# ---
# ## 4. View Creating Entities
# 
# Next, use *creatingEntity* as the identifier, and take a look at what creating entities are.  The *creatingEntity* is related to what instrument or institute created the product.

# In[ ]:


# Show available creatingEntities
identifier = "creatingEntity"
creatingEntities = DataAccessLayer.getIdentifierValues(request, identifier)
print(identifier + ":")
print(list(creatingEntities))


# <a href="#top">Top</a>
# 
# ---
# ## 5. View Sector IDs
# 
# Next, use *sectorID* as the identifier and take a look at what sectors are available.  The *sectorID* is what AWIPS uses to know what geographic section the data is related to.

# In[ ]:


# Show available sectorIDs
identifier = "sectorID"
sectorIDs = DataAccessLayer.getIdentifierValues(request, identifier)
print(identifier + ":")
print(list(sectorIDs))


# <a href="#top">Top</a>
# 
# ---
# ## 6. Create a Satellite Product Tree
# 
# By cycling through all the identifiers, a detailed overview of all available products can be created.  In this example, first *creatingEntity* is used, and then *availableLocationNames*, and *availableParameters* are used to build the product list further. 
# 
# > Note: The identifieres *source* and *sectorID* are not used in this tree, but this is only one way to construct such a product overview.

# In[ ]:


# Construct a full satellite product tree
for entity in creatingEntities:
    print(entity)
    # Create a new request each time through so only one Identifer is set per request
    request = DataAccessLayer.newDataRequest("satellite")
    request.addIdentifier("creatingEntity", entity)
    # Group by available locations
    availableSectors = DataAccessLayer.getAvailableLocationNames(request)
    availableSectors.sort()
    for sector in availableSectors:
        print(" - " + sector)
        request.setLocationNames(sector)
        # Get all available products
        availableProducts = DataAccessLayer.getAvailableParameters(request)
        availableProducts.sort()
        for product in availableProducts:
            print("    - " + product)


# <a href="#top">Top</a>
# 
# ---
# ## 7. Define Map Properties
# 
# In order to plot more than one image, it's easiest to define common logic in a function.  Here, a new function called **make_map** is defined.  This function uses the [matplotlib.pyplot package (plt)](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html) to create a figure and axis.  The coastlines (continental boundaries) are added, along with lat/lon grids.  

# In[ ]:


def make_map(bbox, projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(10,12),
            subplot_kw=dict(projection=projection))
    if bbox[0] is not np.nan:
        ax.set_extent(bbox)
    ax.coastlines(resolution='50m')
    gl = ax.gridlines(draw_labels=True)
    gl.top_labels = gl.right_labels = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax


# <a href="#top">Top</a>
# 
# ---
# ## 8. Plot Image Data!
# 
# For this example, use Channel 13 on the two mesoscale sectors from GOES-East satellite. Create a figure to contain the plot.  Create a new *Data Request* for each sector and set the location and parameters on it.  Limit the data to the most recently acquired image using the *getAvailableTimes* function.  Then use pcolormesh to create a plot from the data in the response in the *GridData* object.
# 
# > Note: You may see a warning appear with a red background, this is okay, and will go away with subsequent runs of the cell.

# In[ ]:


# Specify the sectors - GOES East Mesocales
sectors = ["EMESO-1","EMESO-2"]
# Create a figure to contain all subplots
fig = plt.figure(figsize=(16,7*len(sectors)))

# Cycle through the sectors to create and plot recent data from each one
for i, sector in enumerate(sectors):

    # Create the Ch 13 data request
    request = DataAccessLayer.newDataRequest()
    request.setDatatype("satellite")
    request.setLocationNames(sector)
    request.setParameters("CH-13-10.35um")

    # Get the available times
    utc = datetime.datetime.utcnow()
    times = DataAccessLayer.getAvailableTimes(request)

    # Get the grid data using the latest time
    response = DataAccessLayer.getGridData(request, [times[-1]])
    grid = response[0]
    data = grid.getRawData()
    lons,lats = grid.getLatLonCoords()
    
    # Create the bounding box from the grid data
    bbox = [lons.min(), lons.max(), lats.min(), lats.max()]

    # Draw a new subplot based on the bounding box
    fig, ax = make_map(bbox=bbox)
    
    # Add state boundaries where available
    states = cfeat.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none')
    ax.add_feature(states, linestyle=':')
    
    # Create the color scale
    cs = ax.pcolormesh(lons, lats, data, cmap='coolwarm')
    
    #Create the colorbar and add a label
    cbar = fig.colorbar(cs, shrink=0.6, orientation='horizontal')
    cbar.set_label(sector + " " + grid.getParameter() + " " \
    + str(grid.getDataTime().getRefTime()))


# <a href="#top">Top</a>
# 
# ---
# ## See also
# 
# Documentation for:
# 
# * [awips.DataAccessLayer](http://unidata.github.io/python-awips/api/DataAccessLayer.html#)
# * [awips.PyGridData](http://unidata.github.io/python-awips/api/PyGridData.html)
# * [matplotlib.pyplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html)
# * [matplotlib.pyplot.pcolormesh](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.pcolormesh.html)
# * [matplotlib.pyplot.subplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.subplot.html)
# 
# 
# ### Related Notebooks
# 
# * [python-awips: How to Access Data](https://nbviewer.jupyter.org/github/Unidata/pyaos-ams-2021/blob/master/notebooks/dataAccess/python-awips-HowToAccessData.ipynb)
# 
# <a href="#top">Top</a>