add map spatial analogs

eggshell/plot/plt_ncdata.py

@@ -229,6 +229,72 @@ def uncertainty(resouces, variable=None, ylim=None, title=None, file_extension='
     return output_png
 
 
+def map_spatial_analog(ncfile, variable='dissimilarity', cmap='viridis', title='Spatial analog'):
+    """Return a matplotlib Figure instance showing a map of the dissimilarity measure.
+    """
+    import netCDF4 as nc
+    from eggshell.nc import nc_utils
+    from mpl_toolkits.axes_grid import make_axes_locatable
+    import matplotlib.axes as maxes
+    from cartopy.util import add_cyclic_point
+
+
+    try:
+        var = nc_utils.get_values(ncfile, variable)
+        LOGGER.info('Data loaded')
+
+        lats, lons = nc_utils.get_coordinates(ncfile, variable=variable, unrotate=False)
+
+        if len(lats.shape) == 1:
+            cyclic_var, cyclic_lons = add_cyclic_point(var, coord=lons)
+
+            lons = cyclic_lons.data
+            var = cyclic_var
+
+        with nc.Dataset(ncfile) as D:
+            V = D.variables[variable]
+            lon, lat = map(float, V.target_location.split(','))
+
+        LOGGER.info('Lat and lon loaded')
+
+    except Exception as e:
+        msg = 'Failed to get data for plotting: {0}\n{1}'.format(ncfile, e)
+        LOGGER.exception(msg)
+        raise Exception(msg)
+
+    try:
+        fig = plt.figure(facecolor='w', edgecolor='k')
+        fig.subplots_adjust(top=.95, bottom=.05, left=.03, right=.95)
+
+        ax = plt.axes(
+            projection=ccrs.Robinson(central_longitude=int(np.mean(lons))))
+
+        divider = make_axes_locatable(ax)
+        cax = divider.new_horizontal("4%", pad=0.15, axes_class=maxes.Axes)
+        fig.add_axes(cax)
+
+        ax.plot(lon, lat, marker='o', mfc='#292421', ms=13, transform=ccrs.PlateCarree())
+        ax.plot(lon, lat, marker='o', mfc='#ffffff', ms=7, transform=ccrs.PlateCarree())
+
+        cs = ax.contourf(lons, lats, var, 60,
+                         transform=ccrs.PlateCarree(),
+                         cmap=cmap, interpolation='nearest')
+
+        ax.coastlines(color='k', linewidth=.8)
+        ax.set_title(title)
+
+        cb = plt.colorbar(cs, cax=cax, orientation='vertical')
+        cb.set_label(u"–            Dissimilarity             +")  # ha='left', va='center')
+        cb.set_ticks([])
+
+    except:
+        msg = 'failed to plot graphic'
+        LOGGER.exception(msg)
+
+    LOGGER.info('Plot created and figure saved')
+    return fig
+
+
 # def map_robustness(signal, high_agreement_mask, low_agreement_mask,
 #                    variable=None, cmap='seismic', title=None,
 #                    file_extension='png'):