Merge branch 'main' of github-work:RENCI/adda_for_floodwater

adda/adda.py

@@ -350,10 +350,10 @@ def main(args):
     utilities.log.debug('Wrote ADCIRC dwlc field PKL {}'.format(adcirc_pkl))
 
     # Test using the generic grid and plot to see the generated offset surface. Use the same model as previously generated
-    adc_plot_grid = interpolate_scaled_offset_field.generic_grid()
-    df_plot_transformed = interpolate_scaled_offset_field.interpolation_model_transform(adc_plot_grid, model=model, input_grid_type='grid',pathpoly=pathpoly) 
-    tempdf=pd.DataFrame.from_dict(df_plot_transformed)
-    print(tempdf)
+    plot_grid = interpolate_scaled_offset_field.generic_grid()
+    df_plot_transformed = interpolate_scaled_offset_field.interpolation_model_transform(plot_grid, model=model, input_grid_type='grid',pathpoly=pathpoly) 
+    #tempdf=pd.DataFrame.from_dict(df_plot_transformed)
+    #print(tempdf)
 
     # Write out the model 
     newfilename = io_utilities.get_full_filename_with_subdirectory_prepended(rootdir, iosubdir, 'interpolate_linear_model.h5')
@@ -365,11 +365,16 @@ def main(args):
         utilities.log.error('Could not dump model file to disk '+ newfilename)
 
 ##
-## Apply the model to a 400x500 grid and plot the surface, stations, clamps
+## Apply the model to the regular grid and plot the surface, stations, clamps
 ##
-    newfilename = io_utilities.get_full_filename_with_subdirectory_prepended(rootdir, iosubdir, 'surface.png')
-    adda_visualization_plots.save_plot_model( adc_plot_grid=adc_plot_grid, df_surface=df_plot_transformed, df_stations=df_stations, df_land_control=df_land_controls, df_water_control=df_water_controls, filename=newfilename, plot_now=False)
-    utilities.log.info('Saved debug image to {}'.format(newfilename))
+    #newfilename = io_utilities.get_full_filename_with_subdirectory_prepended(rootdir, iosubdir, 'surface.png')
+    plotfilename = 'surface.png'
+    plottitle=args.current_time # f'{obs_starttime} : {obs_endtime}'
+    adda_visualization_plots.save_plot_model(plot_grid=plot_grid, df_surface=df_plot_transformed, 
+                                             df_stations=df_stations, df_land_control=df_land_controls,
+                                             df_water_control=df_water_controls, filename=plotfilename, 
+                                             plot_now=False, title=plottitle)
+    utilities.log.info('Saved image to {}'.format(plotfilename))
 
     utilities.log.info('Finished')
 

adda/adda_visualization_plots.py

@@ -10,92 +10,114 @@
 
 vmax=0.5
 vmin=-vmax
-xlim=[-99, -58]
-ylim=[7, 50]
+xlim=[-100, -50]
+ylim=[5, 55]
+N=20
+#base_cmap='tab20c' # Or simply use None tab20c is also okay
+#cmap= plt.cm.get_cmap('jet', N)
 
-def save_plot_model(adc_plot_grid=None, df_surface=None, df_stations=None, df_land_control=None, df_water_control=None, 
+def save_plot_model(plot_grid=None, df_surface=None, df_stations=None, df_land_control=None, df_water_control=None, 
                     filename=None, plot_now=False, title=None):
     """
     Wrapper to take generated plt file and save to disk
     """
     if filename is None:
         utilities.log.error('save_plot_model: No filename provided to save the plot')
         return
-    plt_data = plot_model(adc_plot_grid=adc_plot_grid, 
+
+    plt_data = plot_model(plot_grid=plot_grid, 
                           df_surface=df_surface, df_stations=df_stations, 
                           df_land_control=df_land_control, df_water_control=df_water_control, 
                           plot_now=plot_now, title=title)
     plt_data.savefig(filename, bbox_inches='tight')
 
-def plot_model(adc_plot_grid=None, df_surface=None, df_stations=None, df_land_control=None, df_water_control=None, plot_now=True, title=None):
+def plot_model(plot_grid=None, df_surface=None, df_stations=None, df_land_control=None, df_water_control=None, plot_now=True, title=None):
     """
-    Basic plotter to display a 2D extrapolation field. The best images will include
-    not only the surface, but also the stations, and control points as a reference 
+    Basic plotter to display the error field. 
 
-    adc_plot_grid carries the (x,y) AXIS values spanning the data: df_surface
+    adc_plot_grid carries the (x,y) values spanning the data: df_surface
     So if, len(LAT) = 300, len(LON)=400, len(dim(df_surface)) is (400*300)
 
     Parameters:
-        adc_plot_grid: dict with Keys of 'LON','LAT'
+        plot_grid: dict with Keys of 'LON','LAT'
         df_surface: (DataFrame) (optional) with headers 'LON','LAT','VAL'
         df_stations: (DataFrame) (optional) with headers 'LON','LAT','VAL'
         df_land_control: (DataFrame) (optional) with headers 'LON','LAT','VAL'
         df_water_control: (DataFrame) (optional) with headers 'LON','LAT','VAL'
         plot_now: bool. True display the plot. Else not
     Results:
-        A plot in the USA East Coast region
+        A plot 
     """
     coastline=np.loadtxt(os.path.join(os.path.dirname(__file__), "misc", "coarse_us_coast.dat"))
-    #N=16
-    #base_cmap='tab20c' # Or simply use None tab20c is also okay
-    #cmap= plt.cm.get_cmap(base_cmap, N)
-    cmap=plt.cm.jet
+
     #
     fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,10), dpi=144) #, sharex=True)
-    # Set up surface
+    cmap=discrete_cmap(N,'jet')
+
+    # draw surface
     if df_surface is not None:
-        x = adc_plot_grid['LON']
-        y = adc_plot_grid['LAT']
+        x = plot_grid['LON']
+        y = plot_grid['LAT']
         v = df_surface['VAL'].values
         v = v.reshape(-1,len(x)) # Reshapes to LAT column-major format
         mesh = ax.pcolormesh(x, y, v, shading='nearest', cmap=cmap, vmin=vmin, vmax=vmax)
-    # Merge control points
+
+    # add in points
     if df_stations is not None:
         stations_X=df_stations['LON'].values
         stations_Y=df_stations['LAT'].values
         stations_V=df_stations['VAL'].values
-        ax.scatter(stations_X, stations_Y, s=60, marker='o',
+        ax.scatter(stations_X, stations_Y, s=50, marker='o',
                    c=stations_V, cmap=cmap, edgecolor='black',
                    vmin=vmin, vmax=vmax)
+
     if df_land_control is not None:
         land_X=df_land_control['LON'].values
         land_Y=df_land_control['LAT'].values
         land_V=df_land_control['VAL'].values
-        ax.scatter(land_X, land_Y, s=30, marker='o',
-                   c=land_V, cmap=cmap, edgecolor='black',
-                   vmin=vmin, vmax=vmax)
+        ax.scatter(land_X, land_Y, s=30, marker='o', edgecolor='black')
+                   #c=land_V, edgecolor='black',
+                   #vmin=vmin, vmax=vmax)
+
     if df_water_control is not None:
         water_X=df_water_control['LON'].values
         water_Y=df_water_control['LAT'].values
         water_V=df_water_control['VAL'].values
         ax.scatter(water_X, water_Y, s=30, marker='x',c='black')
-                   #,edgecolor='black',
-                   #vmin=vmin, vmax=vmax)
 
-    ax.plot(coastline[:,0],coastline[:,1],color='black',linewidth=.25)
+    ax.plot(coastline[:,0],coastline[:,1],color='black',linewidth=1.5)
+
     ax.axis('equal')
     ax.set_xlim(xlim)
     ax.set_ylim(ylim)
-    ax.set_ylabel('Latitude')
-    ax.set_xlabel('Longitude')
+    ax.set_ylabel('Latitude', fontsize=16)
+    ax.set_xlabel('Longitude', fontsize=16)
+    plt.grid(True)
+
     if title is not None:
-        ax.set_title(title)
+        ax.set_title(title, fontsize=20)
+
+    v1 = np.linspace(vmin, vmax, 11, endpoint=True)
+    cbar=plt.colorbar(mesh,ticks=v1,orientation='vertical')
+    cbar.ax.set_yticklabels(["{:4.2f}".format(i) for i in v1]) # add the labels
 
-    plt.colorbar(mesh,orientation='vertical' )
     if (plot_now):
         plt.show()
     return plt
 
+def discrete_cmap(N, base_cmap='jet_r'):
+    """
+    Create an N-bin discrete colormap from the specified input map
+    """
+
+    # Note that if base_cmap is a string or None, you can simply do
+    #    return plt.cm.get_cmap(base_cmap, N)
+    # The following works for string, None, or a colormap instance:
+
+    base = plt.cm.get_cmap(base_cmap)
+    color_list = base(np.linspace(0, 1, N))
+    cmap_name = base.name + str(N)
+    return base.from_list(cmap_name, color_list, N)
 
 def main(args):
     # Grab the test data

processing/interpolate_scaled_offset_field.py

@@ -323,14 +323,14 @@ def generic_grid():
     """
     utilities.log.debug('Generating points for surface figure.')
     lowerleft_x = -100
-    lowerleft_y =7 
+    lowerleft_y =5 
     res = .1  # resolution in deg
-    nx = 420
-    ny = 450
+    nx = 500
+    ny = 500
     x = np.arange(lowerleft_x, lowerleft_x + nx * res, res)
     y = np.arange(lowerleft_y, lowerleft_y + ny * res, res)
-    adc_coords = {'LON':x.tolist(), 'LAT':y.tolist()}
-    return adc_coords 
+    coords = {'LON':x.tolist(), 'LAT':y.tolist()}
+    return coords 
 
 def test_interpolation_fit(df_source, df_land_controls=None, df_water_controls=None, cv_splits=5, nearest_neighbors=3, stratified_header_name=None):
     """