Simplify creating this combo map

[dopplershift]

This combined map with satellite, radar, and surface data takes entirely too much code to create. Let's simplify this:

# coding: utf-8

# In[1]:

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import num2date


# In[2]:

lat = 39
lon = -89
north = lat + 5
south = lat - 5
west = lon - 5
east = lon + 5


# Get the latest visible satellite image and pull it apart using MetPy's GINI file support.

# In[3]:

from siphon.catalog import TDSCatalog
from metpy.io import GiniFile
from urllib.request import urlopen

sat_cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/satellite/VIS/EAST-CONUS_1km/current/catalog.xml')
ds = list(sat_cat.datasets.values())[0]
gini_file = urlopen(ds.access_urls['HTTPServer'])
gini_ds = GiniFile(gini_file).to_dataset()

data_var = gini_ds.variables['Visible']
x = gini_ds.variables['x'][:]
y = gini_ds.variables['y'][:]
proj_var = gini_ds.variables[data_var.grid_mapping]
time_var = gini_ds.variables['time']
timestamp = num2date(time_var[:].squeeze(), time_var.units)


# In[4]:

import cartopy.crs as ccrs

# Create a Globe specifying a spherical earth with the correct radius
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.earth_radius,
                   semiminor_axis=proj_var.earth_radius)
proj = ccrs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,
                             central_latitude=proj_var.latitude_of_projection_origin,
                             standard_parallels=[proj_var.standard_parallel],
                             globe=globe)


# Using the time from the satellite image, get the NIDS reflectivity product for the radar closest to the center of our box; we query the data from the THREDDS server.

# In[5]:

from siphon.radarserver import RadarServer
from metpy.io import Level3File

rs_cat = TDSCatalog('http://thredds.ucar.edu/thredds/radarServer/catalog.xml')
nexrad3 = RadarServer(rs_cat.catalog_refs['NEXRAD Level III Radar from IDD'].href)

query = nexrad3.query()
query.variables('N0Q')
#query.lonlat_box(east=east, west=west, north=north, south=south)
query.lonlat_point(lon=lon, lat=lat)
query.time(timestamp)

radar_catalog = nexrad3.get_catalog(query)
radar_files = [Level3File(urlopen(ds.access_urls['HTTPServer']))
               for ds in radar_catalog.datasets.values()]


# In[6]:

from siphon.ncss import NCSS

sfc_cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/nws/metar/ncdecoded/catalog.xml')
sfc_ds = sfc_cat.datasets['Feature Collection']
sfc_ncss = NCSS(sfc_ds.access_urls['NetcdfSubset'])

query = sfc_ncss.query()
query.variables('air_temperature', 'dew_point_temperature',
                'cloud_area_fraction', 'inches_ALTIM', 'weather',
                'wind_from_direction', 'wind_speed')
query.lonlat_box(east=east, west=west, north=north, south=south)
query.time(timestamp)
query.accept('netcdf4')

sfc_data = sfc_ncss.get_data(query)


# In[7]:

from metpy.units import units
def with_units(ds, mask, varname):
    var = ds[varname]
    unit_str = var.units
    if unit_str == 'Celsius':
        unit_str = 'degC'
    elif unit_str == 'inches':
        unit_str = 'inHg'
    return units.Quantity(var[mask], unit_str)


# In[8]:

from scipy.spatial import cKDTree

def thin_points(xy, radius, sort_key=None):
    # All points masked initially
    mask = np.ones(station_lons.shape, dtype=np.bool)

    if sort_key:
        # Need in decreasing priority
        sorted_indices = np.argsort(sort_key)[::-1]
    else:
        sorted_indices = np.arange(len(xy))

    # Make our tree
    tree = cKDTree(xy)

    # Loop over all the potential points
    for sort_ind in sorted_indices:
        val = mask[sort_ind]
        # Only proceed if we haven't already excluded this point
        if val:
            # Loop over all the neighbors within the radius
            for neighbor in tree.query_ball_point(xy[sort_ind], radius):
                # Mask them out, but don't mask ourselves
                if neighbor != sort_ind:
                    mask[neighbor] = False

    return mask


# In[9]:

code_map = {'':0, 'HZ':4, 'BR':10,
            '-DZ':51, 'DZ':52, '+DZ':53,
            '-RA':61, 'RA':62, '+RA':63,
            '-TSRA':92, 'TSRA': 92, '+TSRA': 95}
def to_codes(strings):
    for s in strings:
        s = s.decode('ascii').strip()
        if ' ' in s:
            s = s.split()[0]
        if s.startswith('VC'):
            s = s[2:]
        yield code_map[s]


# In[10]:

from metpy.calc import get_wind_components

indexes = sfc_data['stationIndex'][:]
station_lons = sfc_data['longitude'][indexes]
station_lats = sfc_data['latitude'][indexes]

wx = np.ascontiguousarray(sfc_data['weather'][:]).view('S16').squeeze()
wx_codes = list(to_codes(wx))

proj_pts = proj.transform_points(ccrs.PlateCarree(), station_lons,
                                 station_lats)[..., :-1]
mask = thin_points(proj_pts, 100000, sort_key=wx_codes)
curr_wx = np.array(wx_codes)[mask]

station_lons = station_lons[mask]
station_lats = station_lats[mask]

tair = with_units(sfc_data, mask, 'air_temperature')
dewp = with_units(sfc_data, mask, 'dew_point_temperature')
slp = with_units(sfc_data, mask, 'inches_ALTIM')

cloud_cover = 8 * sfc_data['cloud_area_fraction'][mask]
cloud_cover[np.isnan(cloud_cover)] = 9
cloud_cover = cloud_cover.astype(np.int)

speed = with_units(sfc_data, mask, 'wind_speed')
direc = with_units(sfc_data, mask, 'wind_from_direction')
u, v = get_wind_components(speed, direc)

sfc_datad = {'air_temperature': tair, 'dew_point_temperature': dewp,
             'eastward_wind': u, 'northward_wind': v,
             'cloud_coverage': cloud_cover, 'present_weather': curr_wx,
             'air_pressure_at_sea_level': slp}


# In[11]:

import cartopy.feature as cfeat

# Set up a feature for the state/province lines. Tell cartopy not to fill in the polygons
state_boundaries = cfeat.NaturalEarthFeature(category='cultural',
                                             name='admin_1_states_provinces_lakes',
                                             scale='50m', facecolor='none')


# In[12]:

from metpy.plots import ctables
ref_norm, ref_cmap = ctables.registry.get_with_steps('NWSReflectivity', 5, 5)

def plot_nids_file(ax, proj, file):
    datadict = file.sym_block[0][0]
    data = np.ma.array(file.map_data(datadict['data']))
    data[np.isnan(data)] = np.ma.masked

    az = np.array(datadict['start_az'] + [datadict['end_az'][-1]])
    rng = np.linspace(0, file.max_range, data.shape[-1] + 1) * 1000.
    x_proj, y_proj = proj.transform_point(file.lon, file.lat, ccrs.PlateCarree())
    xlocs = rng * np.sin(np.deg2rad(az[:, np.newaxis])) + x_proj
    ylocs = rng * np.cos(np.deg2rad(az[:, np.newaxis])) + y_proj

    return ax.pcolormesh(xlocs, ylocs, data, norm=ref_norm, cmap=ref_cmap,
                         alpha=0.3)


# In[13]:

from metpy.plots import StationPlotLayout
from metpy.plots.wx_symbols import sky_cover, current_weather
import matplotlib.patheffects as mpatheffects
outline = [mpatheffects.Stroke(linewidth=1, foreground='black'),
           mpatheffects.Normal()]

text_style = dict(path_effects=outline, clip_on=True)
#text_style = dict()

simple_layout = StationPlotLayout()
simple_layout.add_barb('eastward_wind', 'northward_wind', 'knots')
simple_layout.add_value('NW', 'air_temperature', units='degC',
                        color='red', **text_style)
simple_layout.add_value('SW', 'dew_point_temperature', units='degC',
                        color='lightgreen', **text_style)
simple_layout.add_value('NE', 'air_pressure_at_sea_level', units='mbar',
                        fmt=lambda v: format(10 * v, '03.0f')[-3:],
                        color='lightblue', **text_style)
simple_layout.add_symbol('C', 'cloud_coverage', sky_cover,
                         color='black')
simple_layout.add_symbol('W', 'present_weather', current_weather,
                         color='black')


# In[14]:

from metpy.plots import StationPlot
from metpy.io.nexrad import is_precip_mode

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

im = ax.imshow(data_var[:], extent=(x[0], x[-1], y[0], y[-1]), origin='upper',
               cmap='Greys_r')
ax.coastlines(resolution='50m', color='black')

# Add country borders with a thick line.
ax.add_feature(cfeat.BORDERS, linewidth='2', edgecolor='black')
ax.add_feature(state_boundaries, linestyle=':')

for nids in radar_files:
    if is_precip_mode(nids.prod_desc.vcp):
        mesh = plot_nids_file(ax, proj, nids)

station_plot = StationPlot(ax, station_lons, station_lats,
                           transform=ccrs.PlateCarree())
simple_layout.plot(station_plot, sfc_datad)

ax.set_extent([west, east, south, north])