new_contents: from pyiem.plot import MapPlot, nwsprecip
  from shapely.wkb import loads
  import psycopg2
  import sys
  import numpy as np
  from matplotlib.patches import Polygon
  from matplotlib.collections import PatchCollection
  import matplotlib.colors as mpcolors
  import matplotlib.pyplot as plt
  
  DBCONN = psycopg2.connect(database='idep', host='iemdb', user='nobody')
  cursor = DBCONN.cursor()
  
  scenario = int(sys.argv[1])
  titles2 = {10: '4',
             11: '3',
             7: '4',
             9: '3'}
  nt = {10: '(no-till)', 11: '(no-till)'}
  threshold = 2
  m = MapPlot(sector='iowa', axisbg='white', nologo=True,
              subtitle='1 Jan 2008 thru 31 Dec 2015',
              title=('UCS %s Year %s Scenario Percent Erosion Loss Change from Baseline'
                     ) % (titles2[scenario], nt.get(scenario, '')))
  
  cursor.execute("""
  with baseline as (
      SELECT huc_12, sum(avg_loss) * 4.463 as loss from results_by_huc12
      where valid between '2008-01-01' and '2016-01-01' and
      scenario = 0 GROUP by huc_12),
  scenario as (
      SELECT huc_12, sum(avg_loss) * 4.463 as loss from results_by_huc12
      where valid between '2008-01-01' and '2016-01-01' and
      scenario = %s GROUP by huc_12),
  agg as (
      SELECT b.huc_12, b.loss as baseline_loss, s.loss as scenario_loss from
      baseline b LEFT JOIN scenario s on (b.huc_12 = s.huc_12))
  
   SELECT ST_Transform(simple_geom, 4326),
   (scenario_loss  - baseline_loss) / baseline_loss * 100. as val
   from huc12 i JOIN agg d on (d.huc_12 = i.huc_12) ORDER by val DESC
  
  """, (scenario, ))
  
  # bins = np.arange(0, 101, 10)
  bins = [-100, -75, -50, -25, 0, 25, 50, 75, 100]
  cmap = plt.get_cmap("BrBG_r")  # nwsprecip()
  cmap.set_under('white')
  cmap.set_over('black')
  norm = mpcolors.BoundaryNorm(bins, cmap.N)
  patches = []
  
  for row in cursor:
      # print "%s,%s" % (row[2], row[1])
      polygon = loads(row[0].decode('hex'))
      a = np.asarray(polygon.exterior)
      x, y = m.map(a[:, 0], a[:, 1])
      a = zip(x, y)
      # c = cmap( norm([float(row[1]),]) )[0]
      c = cmap(norm([float(row[1]), ]))[0]
      p = Polygon(a, fc=c, ec='None', zorder=2, lw=.1)
      patches.append(p)
  
  m.ax.add_collection(PatchCollection(patches, match_original=True))
  m.draw_colorbar(bins, cmap, norm, units='%')
  
  m.drawcounties()
  m.postprocess(filename='test.png')
  
subject: add huc12 mapping app