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

# # A Demonstration of Python's Power
# 
# Here's just a quick demonstration of how to accomplish a pretty interesting task, plotting a satellite image from a remote server, in Python. We won't explain much of what's going on here, but just want to show how much you can accomplish in Python.

# First we just bring in some tools from a variety of Python libraries, using the `import` command. Python is really powerful at assembling various tools together like lego bricks.

# In[1]:


# A whole bunch of imports
get_ipython().run_line_magic('matplotlib', 'inline')
from urllib.request import urlopen

import matplotlib.pyplot as plt
from matplotlib import patheffects
import cartopy.crs as ccrs
import cartopy.feature as cfeat

from metpy.io import GiniFile
from metpy.plots import colortables

from siphon.catalog import TDSCatalog

from xarray import open_dataset


# Now we download the latest water vapor satellite imagery from a remote server and pull out a variety of needed information. This makes use Unidata's MetPy and Siphon libraries.

# In[2]:


# Scan the catalog and download the data

cat = TDSCatalog(
    'http://thredds.ucar.edu/thredds/catalog/satellite/WV/WEST-CONUS_4km/current/catalog.xml'
)
dataset = cat.datasets[list(cat.datasets)[0]]
gini_file = GiniFile(urlopen(dataset.access_urls['HTTPServer']))
gini_ds = open_dataset(gini_file)
# Pull parts out of the data file
data_var = gini_ds.variables['WV']
x = gini_ds.variables['x'][:]
y = gini_ds.variables['y'][:]
time_var = gini_ds.variables['time']
proj_var = gini_ds.variables['projection']

# Set up the projection
globe = ccrs.Globe(
    ellipse='sphere',
    semimajor_axis=proj_var.attrs['earth_radius'],
    semiminor_axis=proj_var.attrs['earth_radius'])
proj = ccrs.LambertConformal(
    central_longitude=proj_var.attrs['longitude_of_central_meridian'],
    central_latitude=proj_var.attrs['latitude_of_projection_origin'],
    standard_parallels=[proj_var.attrs['standard_parallel']],
    globe=globe)


# Now we plot the satellite data along with some mapping information; this plotting is provided by the matplotlib library, while the mapping capabilities come from a library called cartopy.

# In[3]:


# Create the figure
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)

wv_norm, wv_cmap = colortables.get_with_range('WVCIMSS', 100, 260)
wv_cmap.set_under('k')
im = ax.imshow(
    data_var[:], cmap=wv_cmap, norm=wv_norm,
    extent=(x.min(), x.max(), y.min(), y.max()),
    origin='upper')

# Add mapping information
ax.coastlines(resolution='50m', color='black')
state_boundaries = cfeat.NaturalEarthFeature(
    category='cultural',
    name='admin_1_states_provinces_lakes',
    scale='50m',
    facecolor='none')

ax.add_feature(state_boundaries, edgecolor='black', linestyle=':')
ax.add_feature(cfeat.BORDERS, linewidth=2, edgecolor='black')

text = ax.text(
    0.99, 0.01, gini_file.prod_desc.datetime, horizontalalignment='right',
    transform=ax.transAxes, color='white', fontsize='x-large', weight='bold')

text.set_path_effects([
    patheffects.Stroke(linewidth=2, foreground='black'),
    patheffects.Normal()
])