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

# Dependencies | Version
# --- | ---
# satpy | 0.10.0
# pyresample | 1.10.1
# metpy | 0.9.1
# siphon | 0.8.0
# xarray | 0.10
# fiona | 1.7.13
# 
# # Hurricane Florence Plot with SatPy, MetPy, and Siphon with THREDDS Server
# 
# This notebook provides a complex example of how SatPy, MetPy, and Siphon can be used to load data from different sources and combine them in to one image using cartopy and matplotlib.
# 
# This example is based on an [example notebook](https://gist.github.com/dopplershift/48001c3102b1583b78c2c6542618beac) originally created by Ryan May (@dopplershift).

# In[1]:


import metpy.calc as mpcalc
from metpy.plots import StationPlot, add_metpy_logo, add_unidata_logo, add_timestamp
from metpy.units import units
from siphon.catalog import TDSCatalog
from siphon.simplewebservice.ndbc import NDBC

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.patheffects as mpatheffects
import matplotlib.pyplot as plt
import numpy as np

from satpy import Scene
from satpy.writers import get_enhanced_image

from urllib.request import urlopen
import fiona

def get_zip(url):
    data = urlopen(url).read()
    return fiona.BytesCollection(data)


# ## Step 1: Load GOES-16 ABI from a remote THREDDS server using SatPy and Siphon

# In[2]:


base_cat_url = 'http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/{platform}/{sector}/{channel}/current/catalog.xml'
urls = []
# Access ABI channels 1-3 to use in making a true color image (ABI has 16 total channels)
for channel in ['Channel{:02d}'.format(x) for x in range(1, 4)]:
    cat_url = base_cat_url.format(platform='GOES16', sector='CONUS', channel=channel)
    cat = TDSCatalog(cat_url)
    # get the latest dataset from each
    ds = cat.datasets[-1]
    # get the opendap url for this dataset
    urls.append(ds.access_urls['OPENDAP'])    

# Access the files to figure out what is available and load a true_color RGB
scn = Scene(reader='abi_scmi', filenames=urls)
scn.load(['true_color']) 
new_scn = scn.resample(resampler='native')
# Scale the reflectance data to look better as an RGB on a 0 to 1 scale
var = get_enhanced_image(new_scn['true_color']).data
# Get true color data to use later and reorder the dimensions so matplotlib can use the image
# Sadly, this operation is not lazy (bad performance) in xarray at the time of writing
var = var.transpose('y', 'x', 'bands')


# ## Step 2: Load other data using MetPy and Siphon

# In[3]:


best_track = get_zip('https://www.nhc.noaa.gov/gis/best_track/al062018_best_track.zip')
forecast = get_zip('https://www.nhc.noaa.gov/gis/forecast/archive/al062018_5day_latest.zip')
track_x, track_y = np.array(list(zip(*([item['geometry']['coordinates'] for _, item in best_track.items()]
                                       + list(forecast.items())[0][1]['geometry']['coordinates']))))


# In[4]:


buoy = NDBC.latest_observations()
buoy.dropna(subset=['pressure', 'wind_speed', 'wind_direction'], inplace=True)
buoy_u, buoy_v = mpcalc.wind_components(buoy['wind_speed'].values, buoy['wind_direction'].values * units.degree)


# In[5]:


sst_cat = TDSCatalog('https://www.ncei.noaa.gov/thredds/catalog/OisstBase/NetCDF/AVHRR/201809/catalog.xml')
sst = sst_cat.datasets[-1].remote_access(use_xarray=True)
sst_data = sst.metpy.parse_cf('sst')


# ## Step 3: Plot all of the data using Cartopy

# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')
fig = plt.figure(figsize=(20, 10), dpi=200)
abi_crs = var.attrs['area'].to_cartopy_crs()
ax = fig.add_subplot(1, 1, 1, projection=abi_crs)

xy = abi_crs.transform_points(ccrs.PlateCarree(), buoy['longitude'].values, buoy['latitude'].values)
mask = mpcalc.reduce_point_density(xy[..., 0:2], 30000, priority=buoy['pressure'].values)

kwargs = dict(path_effects=[mpatheffects.withStroke(foreground='white', linewidth=2)],
              clip_on=True)
sp = StationPlot(ax, buoy['longitude'].values[mask], buoy['latitude'].values[mask], transform=ccrs.PlateCarree())
sp.plot_parameter('NW', buoy['air_temperature'].values[mask], color='red', **kwargs)
sp.plot_parameter('SW', buoy['dewpoint'].values[mask], color='darkgreen', **kwargs)
sp.plot_parameter('NE', buoy['pressure'].values[mask], color='black', **kwargs)
sp.plot_barb(buoy_u[mask], buoy_v[mask], edgecolor='white')

ax.imshow(var.data, extent=(var.x[0], var.x[-1], var.y[-1], var.y[0]), origin='upper')
ax.contour(sst_data.lon, sst_data.lat, sst_data.squeeze(),
           [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], transform=sst_data.metpy.cartopy_crs)
ax.plot(track_x, track_y, marker='o', color='tab:blue', transform=ccrs.Geodetic())
ax.add_feature(cfeature.COASTLINE.with_scale('10m'), edgecolor='orange')
ax.add_feature(cfeature.STATES.with_scale('10m'), edgecolor='orange')
ax.set_extent([-85, -69, 27, 38], crs=ccrs.PlateCarree())

ax.set_title('GOES-16 Visible, SST (contours), NDBC Buoy Observations, NHC Best and Forecast Track')
add_unidata_logo(fig, y=1375, x=2350, size='large')
add_metpy_logo(fig, y=1375, x=25, size='large')
add_timestamp(ax, high_contrast=True, y=0.01)
plt.savefig('florence.png', dpi=200, bbox_inches='tight')


# <img src="https://media.githubusercontent.com/media/pytroll/pytroll-examples/main/satpy/florence.png" width="1000">

# In[ ]: