<jupyter_start><jupyter_text>Visualization_of_station_location_and_BC_map<jupyter_code>import sys
import glob
import h5py

import numpy as np
import pandas as pd
import netCDF4 as nc
sys.path.insert(0, '/glade/u/home/ksha/WORKSPACE/QC_OBS/')
sys.path.insert(0, '/glade/u/home/ksha/WORKSPACE/utils/')

import data_utils as du
import graph_utils as gu
from namelist import *
# graph tools
import cmaps
import cartopy.crs as ccrs
import cartopy.mpl.geoaxes
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import matplotlib.colors as colors
import matplotlib.patches as patches
from matplotlib.collections import PatchCollection

from matplotlib import ticker
import matplotlib.ticker as mticker
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

%matplotlib inline
RB = [0.25, 0.41, 0.88]  # Royal blue
AM = [0.90, 0.17, 0.31]  # Amathran

JET = []
JET.append(gu.xcolor('indian red'))
JET.append(gu.xcolor('gold'))
JET.append(gu.xcolor('dark sea green'))
JET.append(gu.xcolor('deep sky blue'))
JET.append(gu.xcolor('royal blue'))
need_publish = False

# True: publication quality figures
# False: low resolution figures in the notebook

if need_publish:
    dpi_ = fig_keys['dpi']
else:
    dpi_ = 75<jupyter_output><empty_output><jupyter_text>Elevation Data<jupyter_code>with h5py.File(BACKUP_dir+'CDEM_original.hdf', 'r') as h5io:
    CDEM = h5io['CDEM'][...]
    nav_lon = h5io['CDEM_lon'][...]
    nav_lat = h5io['CDEM_lat'][...]

CDEM[CDEM<=0]=np.nan<jupyter_output><empty_output><jupyter_text>Number of samples to be QC'd per station<jupyter_code>with pd.HDFStore(BACKUP_dir+'BCH_metadata.hdf', 'r') as h5io:
    metadata = h5io['metadata']

stn_code = metadata['code'].values.tolist()
stn_lat = metadata['lat'].values
stn_lon = metadata['lon'].values
stn_elev = metadata['elev [m]'].values
N_sample = {
    'ILL': 5893, 'GOL': 5180, 'HFG': 3114, 'KEY': 2643, 'MSK': 1064, 'PKA': 691, 'PMD': 2508, 'SKI': 176, 'TNS': 1118,
    'WCC': 3339, 'WRU': 63, 'AKI': 2925, 'AKN': 2006, 'ALU': 13756, 'ASH': 6816, 'BAR': 2045, 'BIR': 3808, 'BLN': 6719,
    'BMN': 2062, 'BRI': 4228, 'BUL': 2794, 'BVR': 13368, 'CHK': 4662, 'CHW': 3705, 'CLO': 5359, 'CLW': 7921, 'CMX': 4476,
    'COQ': 14311, 'COX': 4595, 'CQM': 8588, 'DAI': 3870, 'DBC': 5183, 'DCN': 3507, 'DLU': 2724, 'DON': 3288, 'DOW': 10725,
    'EAC': 3457, 'ELK': 4811, 'ELR': 14599, 'ERC': 4426, 'FDL': 4846, 'FER': 4307, 'FIN': 2809, 'FLK': 2954, 'FST': 10357,
    'GLD': 6974, 'GOC': 6265, 'GOL': 5180, 'GRN': 4831, 'GRP': 6854, 'GRT': 13424, 'HEB': 5086, 'HFF': 2844, 'HRN': 3694,
    'HUR': 6596, 'ILL': 5893, 'JHT': 6239, 'KEY': 2643, 'MCQ': 4870, 'MIS': 2200, 'MOB': 4247, 'MOL': 3269, 'MOR': 2811,
    'MTR': 3935, 'NTY': 3461, 'PAK': 5974, 'PAR': 3253, 'PNK': 2849, 'PRS': 3813, 'PUL': 3480, 'PYN': 3522, 'QBY': 1926,
    'SGL': 14298, 'SHH': 3670, 'SLK': 4370, 'STA': 9187, 'STV': 5777, 'UCE': 2865, 'WAH': 5390, 'WOL': 5350, 'WON': 2070,
    'YGS': 1041,'KWA': 4124}
with pd.HDFStore(BACKUP_dir+'BCH_clean/BCH_combine.hdf', 'r') as hdf_temp:
    keys = hdf_temp.keys()
    
keys = du.del_slash(keys)
inds = du.str_search(stn_code, keys, verbose=False)
L = len(stn_code); flag_pick = du.ind_to_flag(inds, L)<jupyter_output><empty_output><jupyter_text>Separate stations by watershed groups<jupyter_code>#rocky = [-119, -114, 48, 54]
south = [-130, -121, 48.75, 50.25]
north = [-127.5, -110, 53, 60]

loc_id = [] # 0 van isl, 1 south, 2 rocky, 3 north

for i in range(len(stn_code)):
    temp_lat = stn_lat[i]
    temp_lon = stn_lon[i]
    if du.check_bounds(temp_lon, temp_lat, south):
        loc_id.append(1)
    elif du.check_bounds(temp_lon, temp_lat, north):
        loc_id.append(3)
    else:
        loc_id.append(2)
        
loc_id = np.array(loc_id)
loc_pick = loc_id[flag_pick]

# ===== Save identified zones
# metadata['zone'] = np.array(loc_id)
# with pd.HDFStore(base_dir+'BCH_metadata_zoned.hdf', 'w') as hdf_temp:
#     hdf_temp['metadata'] = metadata

stn_pick = np.array(stn_code)[flag_pick]
lon_pick = np.array(stn_lon)[flag_pick]
lat_pick = np.array(stn_lat)[flag_pick]

loc_pick = loc_id[flag_pick]

zone1 = {}
zone2 = {}
zone3 = {}
for i, key in enumerate(stn_pick):
    if loc_pick[i] == 1:
        zone1[key] = N_sample[key]
    elif loc_pick[i] == 2:
        zone2[key] = N_sample[key]        
    elif loc_pick[i] == 3:
        zone3[key] = N_sample[key]<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code># Cartopy map settings
scale_param = '10m' # 10m for publication quality
# BC and Ca-US boundary
PROVINCE = cfeature.NaturalEarthFeature(
    category='cultural',
    name='admin_1_states_provinces_lines',
    scale=scale_param,
    facecolor='none')
shape_dir = fig_dir+'wshed_hires/MajorHydroWatershedsProject.shp'
name_list = ['United States of America']
geom_US = gu.get_country_geom(name_list)
topo_level=np.arange(0, 3200, 200)
edge = [-127.5, -114, 48, 59]
s = 1
fig = plt.figure(figsize=(13, 7))
ax1 = plt.subplot2grid((5, 3), (0, 0), rowspan=5, colspan=2, projection=ccrs.PlateCarree())
ax2 = plt.subplot2grid((5, 3), (0, 2))
ax3 = plt.subplot2grid((5, 3), (1, 2))
ax4 = plt.subplot2grid((5, 3), (2, 2))
ax5 = plt.subplot2grid((5, 3), (3, 2))

plt.subplots_adjust(0, 0, 1, 1, hspace=2, wspace=0.05)
# ===== #
ax = ax1
ax.set_extent(edge, ccrs.Geodetic())
ax.add_feature(cfeature.LAND.with_scale(scale_param), facecolor='0.75')
ax.add_feature(cfeature.COASTLINE.with_scale(scale_param), zorder=3)
ax.add_feature(cfeature.BORDERS.with_scale(scale_param), linestyle='--', linewidth=3.5, zorder=3)
ax.add_feature(PROVINCE, edgecolor='k', linestyle='-', linewidth=3.5, zorder=3)
ax.add_feature(cfeature.LAKES.with_scale(scale_param), linewidth=0, facecolor='k', zorder=5)
ax.add_feature(cfeature.RIVERS.with_scale(scale_param), linewidth=0.5, edgecolor='k', zorder=5)
ax.spines['geo'].set_linewidth(2.5)
ax.add_geometries(Reader(shape_dir).geometries(), ccrs.PlateCarree(),
                  facecolor='none', edgecolor='0.25', linewidth=2.5, hatch='///')

ax.add_geometries(geom_US, ccrs.PlateCarree(),
                  facecolor='0.75', edgecolor='none', linewidth=0, zorder=6)
#
GL = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0)
GL.top_labels = False
GL.right_labels = False
GL.xformatter = LONGITUDE_FORMATTER
GL.yformatter = LATITUDE_FORMATTER
GL.xlabel_style = {'size': 14}
GL.ylabel_style = {'size': 14}
GL.ylocator = mticker.FixedLocator([50, 52, 54, 56, 58, 60])
GL.xlocator = mticker.FixedLocator([-126, -124, -122, -120, -118, -116])
CS = ax.contourf(nav_lon[::s, ::s], nav_lat[::s, ::s], CDEM[::s, ::s], topo_level, cmap=plt.cm.terrain, extend='max')

handle1 = ax.plot(lon_pick[loc_pick==1], lat_pick[loc_pick==1], 'ro', mfc=JET[3], ms=12, mec='k', 
                  mew=2.5, zorder=4, label='South Coast stations')
handle2 = ax.plot(lon_pick[loc_pick==2], lat_pick[loc_pick==2], 'ro', mfc=JET[0], ms=12, mec='k', 
                  mew=2.5, zorder=4, label='Southern Interior stations')
handle3 = ax.plot(lon_pick[loc_pick==3], lat_pick[loc_pick==3], 'ro', mfc=JET[1], ms=12, mec='k', 
                  mew=2.5, zorder=4, label='Northeast stations')

handle_patch = patches.Patch(facecolor='none', edgecolor='0.25', linewidth=2.5, hatch='///', label='Watersheds')
#
ax_lg = fig.add_axes([0.663, 1.00325, 0.01, 0.01])
ax_lg.set_axis_off()
LG = ax_lg.legend(handles=[handle1[0], handle2[0], handle3[0], handle_patch], 
               bbox_to_anchor=(0.5, 0.975), ncol=1, prop={'size':14}, fancybox=False)
LG.get_frame().set_facecolor('white')
LG.get_frame().set_edgecolor('k')
LG.get_frame().set_linewidth(2.5)
LG.get_frame().set_alpha(1.0)

ax.set_title('(a) Location of BC Hydro stations', fontsize=14)

ax_base = fig.add_axes([0.6775, 0.035, 0.325, 0.07]) # 0.52, 0.2535, 0.4615, 0.07
[j.set_linewidth(0) for j in ax_base.spines.values()]
ax_base.tick_params(axis='both', left=False, top=False, right=False, bottom=False, \
                    labelleft=False, labeltop=False, labelright=False, labelbottom=False)
cax = inset_axes(ax_base, width='100%', height='40%', borderpad=0, loc=9)
CBar = plt.colorbar(CS, orientation='horizontal', cax=cax) 
CBar.ax.tick_params(axis='x', labelsize=14, direction='in', length=0)
CBar.set_label('Elevation [m]', fontsize=14)
tick_locator = ticker.MaxNLocator(nbins=7)
CBar.locator = tick_locator
CBar.update_ticks()
CBar.outline.set_linewidth(2.5)

# ===== #
AX = [ax2, ax3, ax4, ax5]
N = [len(list(zone1.keys())), len(list(zone2.keys())), len(list(zone3.keys())), 3]
for i, ax in enumerate(AX):
    ax = gu.ax_decorate(ax, True, True)
    ax.grid(False)
    ax.spines["bottom"].set_visible(True)
    ax.spines["left"].set_visible(False)
    ax.spines["right"].set_visible(True)

    ax.xaxis.set_tick_params(labelsize=10)
    ax.ticklabel_format(style='sci', scilimits=(3, 4), axis='y')
    ax.tick_params(axis="both", which="both", bottom=False, top=False, labelbottom=True, 
                   left=False, labelleft=False, right=False, labelright=True)
    ax.set_xlim([-0.5, N[i]-0.25])

ax2.bar(range(N[0]), list(zone1.values()), color=JET[3], linewidth=1.5, edgecolor='k', align='center', width=0.65)
ax3.bar(range(N[1]), list(zone2.values()), color=JET[0], linewidth=1.5, edgecolor='k', align='center', width=0.65)
ax4.bar(range(N[2]), list(zone3.values()), color=JET[1], linewidth=1.5, edgecolor='k', align='center', width=0.65)
ax5.bar(range(3), [sum(list(zone1.values())), sum(list(zone2.values())), sum(list(zone3.values()))], 
        color=[JET[3], JET[0], JET[1]], linewidth=1.5, edgecolor='k', align='center')

ax2.set_xticks(range(N[0]));
ax3.set_xticks(range(N[1]));
ax4.set_xticks(range(N[2]));
ax5.set_xticks(range(3));

ax2.set_ylim([0, 2e4])
ax3.set_ylim([0, 2e4])
ax4.set_ylim([0, 1e4])
ax5.set_ylim([0, 2e5])

ax2.set_xticklabels(list(zone1.keys()), rotation=90);
ax3.set_xticklabels(list(zone2.keys()), rotation=90);
ax4.set_xticklabels(list(zone3.keys()), rotation=90);
ax5.set_xticklabels(['South Coast', 'Southern Interior', 'Northeast']);

ax2.text(0.065, 1.075, ', Sample Number', fontsize=8, transform=ax2.transAxes)

ax2.set_title('(b) South Coast stations', ha='right', x=1, fontsize=14)
ax3.set_title('(c) Southern Interior stations', ha='right', x=1, fontsize=14)
ax4.set_title('(d) Northeast stations', ha='right', x=1, fontsize=14)
ax5.set_title('(e) Sample number by regions', ha='right', x=1, fontsize=14);

if need_publish:
    # Save figure
    fig.savefig(fig_dir+'QC_stn_locs.png', format='png', **fig_keys)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18786127167630057, 'fraction_numerical': 0.07327988066380757, 'mean_word_length': 5.3435836782968655, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '238589', 'license': 'MIT', 'n_tokens_mistral': 4982, 'n_tokens_neox': 4471, 'n_words': 857, 'path': 'publish/PLOT_QC_map_stations.ipynb', 'repo_name': 'yingkaisha/rainbow', 'size': '10726'}