#!/usr/bin/env python3 ''' Name: global_det_atmos_plots_nohrsc_spatial_map.py Contact(s): Mallory Row (mallory.row@noaa.gov) Abstract: This script generates a spatial map for 24 hour NOHRSC snowfall. (lat-lon plots; contours: snowfall) (EVS Graphics Naming Convention: nohrsc.vYYYYmmdd12.024h.conus.[gif][png]) ''' import netCDF4 as netcdf import pyproj import numpy as np import matplotlib matplotlib.use('agg') import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import os import logging import sys import datetime import subprocess import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter from cartopy import config import global_det_atmos_util as gda_util from global_det_atmos_plots_specs import PlotSpecs class NOHRSCSpatialMap: """ Create a NOHRSC spatial map graphic """ def __init__(self, logger, input_dir, DATA_output_dir, COMOUT_output_dir, date_info_dict, plot_info_dict, logo_dir): """! Initalize NOHRSCSpatialMap class Args: logger - logger object input_dir - path to input directory (string) DATA_output_dir - path to DATA output directory (string) COMOUT_output_dir - path to COMOUT output directory (string) plot_info_dict - plot information dictionary (strings) date_info_dict - date information dictionary (strings) logo_dir - directory with logo images (string) Returns: """ self.logger = logger self.input_dir = input_dir self.DATA_output_dir = DATA_output_dir self.COMOUT_output_dir = COMOUT_output_dir self.date_info_dict = date_info_dict self.plot_info_dict = plot_info_dict self.logo_dir = logo_dir def make_nohrsc_spatial_map(self): """! Create the NOHRSCn spatial map graphic Args: Returns: """ self.logger.info(f"Creating NOHRSC spatial map...") self.logger.debug(f"Input directory: {self.input_dir}") self.logger.debug(f"Output directory: {self.DATA_output_dir}") self.logger.debug(f"Date information dictionary: " +f"{self.date_info_dict}") self.logger.debug(f"Plot information dictionary: " +f"{self.plot_info_dict}") # Set valid date valid_date_dt = datetime.datetime.strptime( self.date_info_dict['end_date'] +self.date_info_dict['valid_hr_end'], '%Y%m%d%H' ) # Set contour levels and color map clevs_in = [0, 0.1, 1, 2, 3, 4, 6, 8, 12, 18, 24, 30, 36, 48, 60, 72, 96, 120] colorlist_in = ['#e4eef5', '#bdd7e7', '#6baed6', '#3282bd', '#09519c', '#092694', '#fffc96', '#ffc400', '#ff8700', '#db1f00', '#9d1300', '#690900', '#360200', '#ccccff', '#9f8cd8', '#7c52a5', '#561c72'] cmap_under_color_in = '#ffffff' cmap_over_color_in = '#2e0533' clevs_m = [0.1, 0.5, 1, 2, 3, 4, 5, 10, 15, 20, 25, 35, 50] colorlist_m = ['#a8e7e9', '#70e4ea', '#69b5d5', '#669ce5', '#5065d0', '#4739ce', '#561ec1', '#7300be', '#b705c6', '#970983', '#a01c60', '#c03b60', '#924b4a'] cmap_under_color_m = '#ffffff' cmap_over_color_m = '#5d2c2e' # Set Cartopy shapefile location config['data_dir'] = config['repo_data_dir'] # Convert NOHRSC grib2 to netCDF and read in data self.logger.info(f"Reading in NOHRSC file from {self.input_dir}") nohrsc_grib2_file = os.path.join( self.input_dir, 'nohrsc_snow_24hrAccum_valid'+valid_date_dt.strftime('%Y%m%d%H') +'.grb2' ) nohrsc_netcdf_file = os.path.join( self.input_dir, 'nohrsc_snow_24hrAccum_valid'+valid_date_dt.strftime('%Y%m%d%H') +'.nc' ) check_wgrib2 = subprocess.run( ['which', 'wgrib2'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL ) if check_wgrib2.returncode == 0: self.logger.info(f"Converting {nohrsc_grib2_file} GRIB2 file to " +f"{nohrsc_netcdf_file} netCDF file") run_convert = subprocess.run( ['wgrib2', nohrsc_grib2_file, '-netcdf', nohrsc_netcdf_file] ) else: self.logger.error("wgrib2 executable not in PATH, " "cannot converting file to netCDF") sys.exit(1) make_png = False make_gif = False DATA_png_name = os.path.join( self.DATA_output_dir, 'nohrsc.v'+valid_date_dt.strftime('%Y%m%d%H')+'.024h.' +self.plot_info_dict['vx_mask']+'.png' ) COMOUT_png_name = DATA_png_name.replace( self.DATA_output_dir, self.COMOUT_output_dir ) if not os.path.exists(DATA_png_name): if os.path.exists(nohrsc_netcdf_file): make_png = True else: make_png = False self.logger.info(f"{nohrsc_netcdf_file} does not exist, " +"not making plot") else: make_png = False if make_png and os.path.exists(COMOUT_png_name): gda_util.copy_file(COMOUT_png_name, DATA_png_name) make_png = False if make_png: self.logger.debug("Plotting data from "+nohrsc_netcdf_file) nohrsc_data = netcdf.Dataset(nohrsc_netcdf_file) nohrsc_lat = nohrsc_data.variables['latitude'][:] nohrsc_lon = nohrsc_data.variables['longitude'][:] nohrsc_ASNOW_surface = nohrsc_data.variables['ASNOW_surface'][0,:,:] var_name = (nohrsc_data.variables['ASNOW_surface']\ .getncattr('long_name')) var_level = (nohrsc_data.variables['ASNOW_surface'] .getncattr('level')) var_units = (nohrsc_data.variables['ASNOW_surface']\ .getncattr('units')) if nohrsc_lat.ndim == 1 and nohrsc_lon.ndim == 1: x, y = np.meshgrid(nohrsc_lon, nohrsc_lat) elif nohrsc_lat.ndim == 2 and nohrsc_lon.ndim == 2: x = nohrsc_lon y = nohrsc_lat if var_units == 'm': self.logger.info(f"Converting from {var_units} to " +"inches") nohrsc_ASNOW_surface = nohrsc_ASNOW_surface * 39.3701 var_units = 'inches' nohrsc_ASNOW_surface = np.ma.masked_equal(nohrsc_ASNOW_surface, 0) self.logger.info(f"Doing plot set up") plot_specs_nsm = PlotSpecs(self.logger, 'nohrsc_spatial_map') plot_specs_nsm.set_up_plot() plot_title = ( f'NOHRSC 24 hour {var_name} ({var_units})\n' +'valid ' +(valid_date_dt-datetime.timedelta(hours=24))\ .strftime('%d%b%Y %H')+'Z to ' +valid_date_dt.strftime('%d%b%Y %H')+'Z' ) plot_left_logo = False plot_left_logo_path = os.path.join(self.logo_dir, 'noaa.png') if os.path.exists(plot_left_logo_path): plot_left_logo = True left_logo_img_array = matplotlib.image.imread( plot_left_logo_path ) left_logo_xpixel_loc,left_logo_ypixel_loc,left_logo_alpha = ( plot_specs_nsm.get_logo_location( 'left', plot_specs_nsm.fig_size[0], plot_specs_nsm.fig_size[1], plt.rcParams['figure.dpi'] ) ) plot_right_logo = False plot_right_logo_path = os.path.join(self.logo_dir, 'nws.png') if os.path.exists(plot_right_logo_path): plot_right_logo = True right_logo_img_array = matplotlib.image.imread( plot_right_logo_path ) right_logo_xpixel_loc,right_logo_ypixel_loc,right_logo_alpha = ( plot_specs_nsm.get_logo_location( 'right', plot_specs_nsm.fig_size[0], plot_specs_nsm.fig_size[1], plt.rcParams['figure.dpi'] ) ) if var_units == 'inches': clevs = clevs_in cmap = matplotlib.colors.ListedColormap(colorlist_in) cmap_under_color = cmap_under_color_in cmap_over_color = cmap_over_color_in elif var_units == 'm': clevs = clevs_m cmap = matplotlib.colors.ListedColormap(colorlist_m) cmap_under_color = cmap_under_color_m cmap_over_color = cmap_over_color_m norm = matplotlib.colors.BoundaryNorm(clevs, cmap.N) # Create plot self.logger.info(f"Creating plot for NOHRS") fig = plt.figure(figsize=(plot_specs_nsm.fig_size[0], plot_specs_nsm.fig_size[1])) gs_hspace, gs_wspace = 0, 0 gs_bottom, gs_top = 0.125, 0.85 gs = gridspec.GridSpec(1,1, bottom=gs_bottom, top=gs_top, hspace=gs_hspace, wspace=gs_wspace) fig.suptitle(plot_title) if plot_left_logo: left_logo_img = fig.figimage( left_logo_img_array, left_logo_xpixel_loc, left_logo_ypixel_loc, zorder=1, alpha=right_logo_alpha ) left_logo_img.set_visible(True) if plot_right_logo: right_logo_img = fig.figimage( right_logo_img_array, right_logo_xpixel_loc, right_logo_ypixel_loc, zorder=1, alpha=right_logo_alpha ) if self.plot_info_dict['vx_mask'] == 'conus': extent = [-124,-70,18.0,50.0] central_lon = -97.6 central_lat = 35.4 elif self.plot_info_dict['vx_mask'] == 'alaska': extent = [-180,-110,45.0,75.0] central_lon = -145 central_lat = 60 myproj=ccrs.LambertConformal(central_longitude=central_lon, central_latitude=central_lat, false_easting=0.0, false_northing=0.0, globe=None) ax1 = fig.add_subplot(gs[0], projection=myproj) ax1.set_extent(extent) ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), zorder=2, linewidth=1) ax1.add_feature(cfeature.BORDERS.with_scale('50m'), zorder=2, linewidth=1) ax1.add_feature(cfeature.STATES.with_scale('50m'), zorder=2, linewidth=1) CF1 = ax1.contourf(x, y, nohrsc_ASNOW_surface, transform=ccrs.PlateCarree(), levels=clevs, norm=norm, cmap=cmap, extend='max') CF1.cmap.set_over(cmap_under_color) CF1.cmap.set_over(cmap_over_color) cbar_left = gs.get_grid_positions(fig)[2][0] cbar_width = (gs.get_grid_positions(fig)[3][-1] - gs.get_grid_positions(fig)[2][0]) cbar_bottom = 0.075 cbar_height = 0.03 cbar_ax = fig.add_axes( [cbar_left, cbar_bottom, cbar_width, cbar_height] ) cbar = fig.colorbar(CF1, cax=cbar_ax, orientation='horizontal', ticks=CF1.levels) cbar.ax.set_xlabel(f'Total Snowfall ({var_units})', labelpad = 0) cbar.ax.xaxis.set_tick_params(pad=0) cbar_tick_labels_list = [] for tick in cbar.get_ticks(): if str(tick).split('.')[1] == '0': cbar_tick_labels_list.append( str(int(tick)) ) else: cbar_tick_labels_list.append( str(round(tick,3)).rstrip('0') ) cbar.ax.set_xticklabels(cbar_tick_labels_list) self.logger.info("Saving image as "+DATA_png_name) plt.savefig(DATA_png_name) plt.clf() plt.close('all') gda_util.copy_file(DATA_png_name, COMOUT_png_name) DATA_gif_name = DATA_png_name.replace('.png', '.gif') COMOUT_gif_name = COMOUT_png_name.replace('.png', '.gif') if os.path.exists(COMOUT_gif_name): gda_util.copy_file(COMOUT_gif_name, DATA_gif_name) make_gif = False elif os.path.exists(DATA_png_name) \ and not os.path.exists(DATA_gif_name): make_gif = True if make_gif: # Convert png to gif, if possible check_convert = subprocess.run( ['which', 'convert'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL ) if check_convert.returncode == 0: self.logger.info(f"Converting {DATA_png_name} to " +f"{DATA_gif_name}") run_convert = subprocess.run( ['convert', DATA_png_name, DATA_gif_name] ) gda_util.copy_file(DATA_gif_name, COMOUT_gif_name) else: self.logger.warning("convert executable not in PATH, " "not creating gif of image") def main(): # Need settings INPUT_DIR = os.environ['HOME'] DATA_OUTPUT_DIR = os.environ['HOME'] COMOUT_OUTPUT_DIR = os.environ['HOME'] LOGO_DIR = os.environ['HOME'] DATE_INFO_DICT = { 'date_type': 'DATE_TYPE', 'start_date': 'START_DATE', 'end_date': 'END_DATE', 'valid_hr_start': 'VALID_HR_START', 'valid_hr_end': 'VALID_HR_END', 'valid_hr_inc': 'VALID_HR_INC', } PLOT_INFO_DICT = { 'line_type': 'LINE_TYPE', 'stat': 'STAT', 'vx_mask': 'VX_MASK', 'obs_var_name': 'OBS_VAR_NAME', 'obs_var_level': 'OBS_VAR_LEVEL', } # Create OUTPUT_DIR gda_util.make_dir(OUTPUT_DIR) # Set up logging logging_dir = os.path.join(OUTPUT_DIR, 'logs') gda_util.make_dir(logging_dir) job_logging_file = os.path.join(logging_dir, os.path.basename(__file__)+'_runon' +datetime.datetime.now()\ .strftime('%Y%m%d%H%M%S')+'.log') logger = logging.getLogger(job_logging_file) logger.setLevel('DEBUG') formatter = logging.Formatter( '%(asctime)s.%(msecs)03d (%(filename)s:%(lineno)d) %(levelname)s: ' + '%(message)s', '%m/%d %H:%M:%S' ) file_handler = logging.FileHandler(job_logging_file, mode='a') file_handler.setFormatter(formatter) logger.addHandler(file_handler) logger_info = f"Log file: {job_logging_file}" print(logger_info) logger.info(logger_info) p = NOHRSCSpatialMap(logger, INPUT_DIR, DATA_OUTPUT_DIR, COMOUT_OUTPUT_DIR, DATE_INFO_DICT, PLOT_INFO_DICT, LOGO_DIR) p.make_nohrsc_spatial_map() if __name__ == "__main__": main()