def main():

    grid_config_to_dir_file = OrderedDict([
        #(default_domains.bc_mh_044, Path("/RESCUE/skynet3_rech1/huziy/directions_for_ManitobaHydro/directions_mh_0.44deg.nc")),
        (default_domains.bc_mh_022, Path("/RESCUE/skynet3_rech1/huziy/directions_for_ManitobaHydro/directions_mh_0.22deg.nc")),
        # (default_domains.bc_mh_011, Path("/RESCUE/skynet3_rech1/huziy/directions_for_ManitobaHydro/directions_bc-mh_0.11deg_new.nc")),
    ])


    stations = stfl_stations.load_stations_from_csv(selected_ids=None)

    print(stations)

    gs = GridSpec(1, len(grid_config_to_dir_file))

    plot_utils.apply_plot_params(width_cm=25, height_cm=20, font_size=8)
    fig = plt.figure()

    for col, (grid_config, dir_path) in enumerate(grid_config_to_dir_file.items()):
        ax = fig.add_subplot(gs[0, col])
        plot_station_positions(directions_file_path=dir_path, stations=stations, ax=ax, grid_config=grid_config)



    img_file = img_folder / "{}_{}.png".format("mh", "_".join([str(gc.dx) for gc in grid_config_to_dir_file]))
    fig.savefig(str(img_file), bbox_inches="tight", dpi=300)

def plot_point_positions_with_upstream_areas(processed_stations, processed_model_points,
                                             basemap, cell_manager, lake_fraction_field=None):
    # plot point positions with upstream areas


    # TODO: plot zonal average lake fractions and plot it at the side of the map
    rc_params_backup = plt.rcParams.copy()
    plot_utils.apply_plot_params(font_size=10, width_pt=None, width_cm=18, height_cm=8)
    fig = plt.figure()

    gs = gridspec.GridSpec(1, 1, width_ratios=[1, 1], wspace=0.01)

    ax = fig.add_subplot(gs[0, 0])
    plot_positions_of_station_list(ax, processed_stations, processed_model_points, basemap, cell_manager)
    # ax.set_aspect("auto")

    # ax = fig.add_subplot(gs[0, 1])
    # ydata = range(lake_fraction_field.shape[1])
    # ax.plot(lake_fraction_field.mean(axis=0) * 100, ydata, lw=2)
    # ax.set_xlabel("%")
    # ax.set_ylim(min(ydata), max(ydata))

    for tl in ax.yaxis.get_ticklabels():
        tl.set_visible(False)

    impath = os.path.join(images_folder, "station_positions.png")
    fig.savefig(impath, bbox_inches="tight", dpi=cpp.FIG_SAVE_DPI, transparent=True)
    plt.close(fig)
    plt.rcParams.update(rc_params_backup)

def main():

    plot_utils.apply_plot_params(width_cm=20, height_cm=20, font_size=10)

    high_hles_years = [1993, 1995, 1998]
    low_hles_years = [1997, 2001, 2006]
    data_path = "/BIG1/skynet1_rech1/diro/sample_obsdata/eraint/eraint_uvslp_years_198111_201102_NDJmean_ts.nc"





    with xr.open_dataset(data_path) as ds:
        print(ds)


        u = get_composit_for_name(ds, "u10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        v = get_composit_for_name(ds, "v10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        msl = get_composit_for_name(ds, "msl", high_years_list=high_hles_years, low_years_list=low_hles_years)

        lons = ds["longitude"].values
        lats = ds["latitude"].values

        print(lats)
        print(msl.shape)
        print(lons.shape, lats.shape)


        lons2d, lats2d = np.meshgrid(lons, lats)


    fig = plt.figure()

    map = Basemap(llcrnrlon=-130, llcrnrlat=22, urcrnrlon=-28,
                  urcrnrlat=65, projection='lcc', lat_1=33, lat_2=45,
                  lon_0=-95, resolution='i', area_thresh=10000)


    clevs = np.arange(-11.5, 12, 1)
    cmap = cm.get_cmap("bwr", len(clevs) - 1)
    bn = BoundaryNorm(clevs, len(clevs) - 1)

    x, y = map(lons2d, lats2d)
    im = map.contourf(x, y, msl / 100, levels=clevs, norm=bn, cmap=cmap) # convert to mb (i.e hpa)
    map.colorbar(im)


    stride = 2
    ux, vy = map.rotate_vector(u, v, lons2d, lats2d)
    qk = map.quiver(x[::stride, ::stride], y[::stride, ::stride], ux[::stride, ::stride], vy[::stride, ::stride],
               scale=10, width=0.01, units="inches")
    plt.quiverkey(qk, 0.5, -0.1, 2, "2 m/s", coordinates="axes")


    map.drawcoastlines(linewidth=0.5)
    map.drawcountries()
    map.drawstates()
    #plt.show()

    fig.savefig("hles_wind_compoosits.png", bbox_inches="tight", dpi=300)

def main():
    # Show selected domains, basins, and/or flow directions or flow accumulations


    mh_gc044 = default_domains.gc_cordex_na_044.subgrid(20, 60, di=130, dj=110)
    mh_gc022 = mh_gc044.double_resolution_keep_free_domain_same()

    test_bc_011 = default_domains.gc_cordex_na_011.subgrid(12, 244, di=404, dj=380)
    test_bc_044 = test_bc_011.decrease_resolution_keep_free_domain_same(4)



    print(test_bc_044)

    plot_utils.apply_plot_params()
    # show_domain(mh_gc044)
    # show_domain(mh_gc022)
    # show_domain(mh_gc011)


    print(test_bc_011)

    # fig, ax, bmp = show_domain(default_domains.gc_cordex_011, draw_rivers=False)
    # show_domain(test_bc_011, draw_rivers=False, show_GRDC_basins=True)
    # show_domain(test_bc_044, draw_rivers=False, show_GRDC_basins=True)


    show_domain(default_domains.bc_mh_044, show_GRDC_basins=True, show_Churchil_Nelson_basins=False)
    # show_domain(default_domains.bc_mh_011)

    plt.show()

    def plot_comparisons_of_seasonal_sst_with_homa_obs(self, start_year=None, end_year=None, season_to_months=None,
                                                       exp_label=""):

        model_data = self.get_seasonal_mean_lst(season_to_months=season_to_months,
                                                start_year=start_year, end_year=end_year)

        obs_sst_path = os.path.expanduser("~/skynet3_rech1/nemo_obs_for_validation/GreatLakes_2003_5km-2/sst-glk.nc")

        obs_data = self.read_and_interpolate_homa_data(path=obs_sst_path, start_year=start_year, end_year=end_year,
                                                       season_to_months=season_to_months)

        plot_utils.apply_plot_params(font_size=10, width_pt=None, width_cm=20, height_cm=10)
        # calculate climatologic differences
        diff = {}
        for season in list(season_to_months.keys()):
            diff[season] = np.mean(
                [model_data[y][season] - obs_data[y][season] for y in range(start_year, end_year + 1)], axis=0)
            diff[season] = np.ma.masked_where(~self.lake_mask, diff[season])
            the_field = diff[season]
            print("diff stats({}): min={}; max={}; avg={}".format(
                season, the_field.min(), the_field.max(), the_field.mean()))


        # plot seasonal biases
        xx, yy = self.basemap(self.lons.copy(), self.lats.copy())


        # calculate difference ranges
        diff_max = 0
        for season, the_diff in diff.items():
            diff_max = max(np.percentile(np.abs(the_diff[~the_diff.mask]), 90), diff_max)
        diff_max = 5

        locator = MaxNLocator(nbins=12, symmetric=True)
        bounds = locator.tick_values(-diff_max, diff_max)
        bn = BoundaryNorm(bounds, len(bounds) - 1)
        cmap = cm.get_cmap("RdBu_r", len(bounds) - 1)

        im = None
        fig = plt.figure()
        ncols = 2
        # fig.suptitle(r"LST $\left({\rm ^\circ C}\right)$", font_properties=FontProperties(weight="bold"))
        gs = GridSpec(len(season_to_months) // ncols, ncols + 1, width_ratios=[1.0, ] * ncols + [0.05, ])
        for i, season in enumerate(season_to_months.keys()):
            ax = fig.add_subplot(gs[i // ncols, i % ncols])
            im = self.basemap.pcolormesh(xx, yy, diff[season][:], ax=ax, cmap=cmap, norm=bn)
            ax.set_title(season)
            self.basemap.drawcoastlines(ax=ax, linewidth=0.5)
            if not i:
                ax.set_ylabel("NEMO - Obs")

        cb = plt.colorbar(im, ticks=locator, cax=fig.add_subplot(gs[:, -1]), extend="both")

        nemo_img_dir = "nemo"
        if not os.path.isdir(nemo_img_dir):
            os.mkdir(nemo_img_dir)

        plt.tight_layout()
        fig.savefig(os.path.join(nemo_img_dir, "sst_homa_validation_{}.pdf".format(exp_label)))
        plt.show()

def main():

    plot_utils.apply_plot_params()
    diag_folder = "/RECH2/huziy/BC-MH/bc_mh_044deg/Diagnostics"

    gc = default_domains.bc_mh_044
    plot_monthly_clim_in_a_panel(diag_folder=diag_folder, grid_config=gc, basins_of_interest_shp=default_domains.MH_BASINS_PATH)

def main():
    plot_utils.apply_plot_params(font_size=20, width_pt=None, width_cm=20, height_cm=20)
    i_start, j_start = 0, 0
    i_end, j_end = -1, -1

    print(get_section_hor_indices(i_start=5, j_start = 0, i_end = 10, j_end = 20))


    # #Superior
    # params = dict(
    #     i_start = 25, j_start = 65,
    #     i_end = 70, j_end = 65,
    #     var_name = "votemper"
    #
    # )
    # plot_cross_section_for_seasons(data_path=T_FILE_PATH, **params)
    #
    # #Michigan
    # params = dict(
    #     i_start = 55, j_start = 55,
    #     i_end = 55, j_end = 5,
    #     var_name = "votemper"
    #
    # )
    # plot_cross_section_for_seasons(data_path=T_FILE_PATH, **params)

    #Huron
    params = dict(
        i_start = 10, j_start = 30,
        i_end = 30, j_end = 10,
        var_name = "votemper"

    )
    plot_cross_section_for_seasons(data_path=T_FILE_PATH, **params)
    plt.show()

def main():
    dh = DataHolder(var_name = "tos")
    plot_utils.apply_plot_params(width_pt=None, font_size=9)

    current_start_date = datetime(1970,1,1)
    current_end_date = datetime(1999,12,31)

    future_start_date = datetime( 2041, 1, 1 )
    future_end_date = datetime( 2070, 12, 31)

    plt.subplot(2,1,1)
    plt.title("SST: current CV")

    cv = dh.get_cv_for_seasonal_mean(start_date = current_start_date,
                                     end_date = current_end_date, months=range(3,8))
    dh.plot_2d_field(field_2d= cv,
                        color_levels=np.arange(0, 0.001, 0.0001)
                    )

    plt.subplot(2,1,2)
    plt.title("SST: future CV")
    cv = dh.get_cv_for_seasonal_mean(current = False,
                                        start_date=future_start_date,
                                        end_date = future_end_date, months=range(3,8))
    dh.plot_2d_field(field_2d = cv,
                        color_levels=np.arange(0, 0.001, 0.0001)
                    )




    plt.savefig("sst_cv.png")
    pass

def plot_alt_from_monthly_climatologies():
    plot_utils.apply_plot_params(width_pt=None, height_cm=20, width_cm=16, font_size=12)
    figure = plt.figure()
    b, lons2d, lats2d = draw_regions.get_basemap_and_coords(llcrnrlat=40.0, llcrnrlon=-145, urcrnrlon=-10)
    x, y = b(lons2d, lats2d)
    permafrost_mask = draw_regions.get_permafrost_mask(lons2d, lats2d)
    dm = CRCMDataManager(data_folder="data/CORDEX")

    year_ranges = [list(range(1981, 1985)), list(range(2041, 2045)), list(range(2071, 2075))]

    gs = gridspec.GridSpec(len(year_ranges), 1)

    pf_mask = (permafrost_mask == 1) | (permafrost_mask == 2)
    pf_mask = ~pf_mask

    permafrost_mask = np.ma.masked_where(permafrost_mask <= 0, permafrost_mask)
    for i, year_range in enumerate(year_ranges):
        ax = figure.add_subplot(gs[i, 0])
        alt = dm.get_alt_using_monthly_mean_climatology(year_range)
        alt = np.ma.masked_where(pf_mask, alt)

        img = b.contourf(x, y, alt, levels=range(11), cmap=cm.get_cmap("jet", ), ax=ax)
        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", "5%", pad="3%")
        cb = plt.colorbar(img, cax=cax)

        b.contour(x, y, permafrost_mask, levels=range(5), linewidth=0.1, colors="k", ax=ax)
        b.drawcoastlines(ax=ax, linewidth=0.5)
        ax.set_title("period: {0} - {1}".format(year_range[0], year_range[-1]))
    plt.savefig("alt_from_clim.png")

def plot_alt_for_different_e_scenarios():
    labels = ("E1", "E2", "E3", "E4")
    p_format = "pmNorthAmerica_0.44deg_CanHistoE{0}"
    prefixes = [p_format.format(x) for x in range(1, 5)]
    b, lons2d, lats2d = draw_regions.get_basemap_and_coords(llcrnrlat=40.0, llcrnrlon=-145, urcrnrlon=-10)
    x, y = b(lons2d, lats2d)
    permafrost_mask = draw_regions.get_permafrost_mask(lons2d, lats2d)

    dm = CRCMDataManager(data_folder="data/CORDEX")  # needed here only for verticla levels

    plot_utils.apply_plot_params(width_pt=None, height_cm=12, width_cm=16, font_size=12)
    fig = plt.figure()

    gs = gridspec.GridSpec(2, 2)

    scen_index = 0
    for row in range(2):
        for col in range(2):
            sc = labels[scen_index]
            ax = fig.add_subplot(gs[row, col])
            h = dm.get_alt_using_files_in(folder="data/CORDEX/na/means_month", file_name_prefix=prefixes[scen_index])
            h = np.ma.masked_where((permafrost_mask == 0) |
                                   (permafrost_mask >= 3) | (h < 0), h)

            plot_for_scenario(sc, ax, basemap=b, x=x, y=y, alt=h, permafrost_mask=permafrost_mask,
                              start_year=1950, end_year=1954
            )

            scen_index += 1
    gs.tight_layout(fig, h_pad=0.9, w_pad=16)
    fig.savefig("alt_diff_scenarios.png")
    pass

def validate_seasonal_mean_atm_fields():
    plot_utils.apply_plot_params(font_size=10, width_pt=None, width_cm=20, height_cm=20)
    from crcm5.analyse_hdf import validate_model_fields

    p = Process(target=validate_model_fields.do_4_seasons, kwargs=dict(
        start_year=1980, end_year=2010))
    p.start()

def plot_temperature_biases():

    seasons = ["(a) Annual", " (b) Winter (DJF)", "(c) Spring (MAM)", "(d) Summer (JJA)", "(e) Fall (SON)"]
    months =  [range(1, 13), [12, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]  ]
    season_to_months = dict(zip(seasons, months))

    cru_var_name = "tmp"
    cru_data_store, crcm4_data_store = _get_comparison_data(cru_var_name= cru_var_name,
        season_to_months=season_to_months)

    x, y = polar_stereographic.xs, polar_stereographic.ys
    i_array, j_array = _get_routing_indices()
    x_min, x_max, y_min, y_max = plot_utils.get_ranges(x[i_array, j_array], y[i_array, j_array])



    plot_utils.apply_plot_params(width_pt= None, font_size=9, aspect_ratio=2.5)
    fig = plt.figure()
    assert isinstance(fig, Figure)

    basemap = polar_stereographic.basemap
    assert isinstance(basemap, Basemap)
    gs = gridspec.GridSpec(3,2)

    color_map = my_cm.get_red_blue_colormap(ncolors = 14, reversed=True)
    clevels = xrange(-8, 9, 2)
    all_plot_axes = []
    for i, season in enumerate(seasons):
        if not i:
            ax = fig.add_subplot(gs[0,:])
        else:
            row, col = (i - 1)  // 2 + 1, (i - 1) % 2
            ax = fig.add_subplot(gs[row, col])
        all_plot_axes.append(ax)
        assert isinstance(ax, Axes)
        delta = crcm4_data_store[season] - cru_data_store[season]
        if cru_var_name == "tmp": delta -= 273.15
        #delta = maskoceans(polar_stereographic.lons, polar_stereographic.lats, delta)
        save = delta[i_array, j_array]
        delta[:, :] = np.ma.masked
        delta[i_array, j_array] = save
        img = basemap.pcolormesh(x, y, delta, cmap = color_map, vmin = -7, vmax = 7)
        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", "8%", pad="3%")
        int_ticker = LinearLocator(numticks = color_map.N + 1)
        fig.colorbar(img, cax = cax, ticks = MultipleLocator(base = 2))
        ax.set_title(season)

    for the_ax in all_plot_axes:
        the_ax.set_xlim(x_min, x_max)
        the_ax.set_ylim(y_min, y_max)
        basemap.drawcoastlines(ax = the_ax, linewidth = 0.1)
        plot_utils.draw_meridians_and_parallels(basemap, step_degrees=30.0, ax = the_ax)
        plot_basin_boundaries_from_shape(basemap, axes = the_ax, linewidth=0.4)
        put_selected_stations(the_ax, basemap, i_array, j_array)

    #gs.tight_layout(fig)
    fig.suptitle("T(2m), degrees, CRCM4 - CRU")
    fig.savefig("seasonal_{0}_ccc.png".format(cru_var_name))

def compare_averages_over_basins():
    print getBasinNames()
    plot_utils.apply_plot_params(width_pt = 400, font_size = 15)

    basinNames = getSelectedBasinNames()
    for basinName in basinNames:
        compare_means(basinName = basinName)
        print basinName

def plot_seasonal_circulation_as_subplots(start_year=1995, end_year=2010,
                                          data_dir="/HOME/huziy/skynet3_rech1/NEMO_OFFICIAL/dev_v3_4_STABLE_2012/NEMOGCM/CONFIG/GLK_LIM3/EXP_GLK_LIM3_1980/zdf_gls_dt_and_sbc_30min",
                                          season_to_months=None, level_index=0):




    # get the data
    lons, lats, season_to_fields = get_seasonal_flows(data_dir=data_dir, start_year=start_year, end_year=end_year, season_to_months=season_to_months, level=level_index)



    img_folder = Path("nemo/circulation_plots/" + Path(data_dir).name)

    if not img_folder.exists():
        img_folder.mkdir(parents=True)

    img_name = "circ_{}-{}_lev{}.png".format(start_year, end_year, level_index)


    nrows = 3
    nsubplots = len(season_to_months)
    ncols = ceil(nsubplots / (nrows - 1))
    gs = GridSpec(nrows, ncols, wspace=0.01, hspace=0, height_ratios=[1,] * (nrows - 1) + [0.05, ])

    plot_utils.apply_plot_params(font_size=8, width_cm=8 * ncols, height_cm=min(4.5 * (nrows - 1), 25))

    fig = plt.figure()

    plot_ind = 0
    for season in season_to_months:

        row = plot_ind // ncols
        col = plot_ind % ncols

        uu, vv = season_to_fields[season]
        flow_speed = (uu ** 2 + vv ** 2) ** 0.5

        uu1, vv1 = uu / flow_speed, vv / flow_speed

        ax = fig.add_subplot(gs[row, col])

        assert isinstance(ax, Axes)
        ax.set_frame_on(False)

        ax.text(0.01, 0.1, season, va="bottom", ha="left", fontsize=12, transform=ax.transAxes)
        im = plot_flow_vectors_basemap(lons=lons, lats=lats, uu=uu1, vv=vv1, flow_speed=flow_speed, ax=ax,
                                       draw_colorbar=(col == 0) and (row == nrows - 2),
                                       streamplot=False)

        plot_ind += 1


    # plt.colorbar(im, cax=fig.add_subplot(gs[nrows - 1, 0]), orientation="horizontal")

    img_path = img_folder / img_name
    fig.savefig(str(img_path), bbox_inches="tight", dpi=300)

def prepare():
    import application_properties

    application_properties.set_current_directory()

    if not img_folder.is_dir():
        img_folder.mkdir(parents=True)

    plot_utils.apply_plot_params(font_size=10, width_cm=20, height_cm=18)

def plot_comparison_hydrographs(basin_name_to_out_indices_map, rea_config=None, gcm_config=None):
    """

    :type basin_name_to_out_indices_map: dict
    """
    assert isinstance(rea_config, RunConfig)
    assert isinstance(gcm_config, RunConfig)

    assert hasattr(rea_config, "data_daily")
    assert hasattr(gcm_config, "data_daily")

    bname_to_indices = OrderedDict([item for item in sorted(basin_name_to_out_indices_map.items(),
                                                            key=lambda item: item[1][1], reverse=True)])

    print(bname_to_indices)

    plot_utils.apply_plot_params(font_size=12, width_pt=None, width_cm=25, height_cm=12)
    fig = plt.figure()
    ncols = 3
    nrows = len(bname_to_indices) // ncols + int(len(bname_to_indices) % ncols != 0)
    gs = GridSpec(nrows, ncols)

    ax_last = None
    for pl_index, (bname, (i, j)) in enumerate(bname_to_indices.items()):
        row = pl_index // ncols
        col = pl_index % ncols
        ax = fig.add_subplot(gs[row, col])

        ax.plot(rea_config.data_daily[0], rea_config.data_daily[1][:, i, j], color="b", lw=2,
                label=rea_config.label)
        ax.plot(gcm_config.data_daily[0], gcm_config.data_daily[1][:, i, j], color="r", lw=2,
                label=gcm_config.label)


        ax.xaxis.set_major_locator(MonthLocator())
        ax.xaxis.set_minor_locator(MonthLocator(bymonthday=15))
        ax.xaxis.set_minor_formatter(FuncFormatter(lambda d, pos: num2date(d).strftime("%b")[0]))
        plt.setp(ax.xaxis.get_majorticklabels(), visible=False)
        ax.grid()

        sfmt = ScalarFormatter(useMathText=True)
        sfmt.set_powerlimits([-2, 2])
        ax.yaxis.set_major_formatter(sfmt)

        bbox_props = dict(boxstyle="round,pad=0.3", fc="cyan", ec="b", lw=1, alpha=0.5)
        ax.annotate(bname, (0.9, 0.1), xycoords="axes fraction", bbox=bbox_props, zorder=10,
                    alpha=0.5, horizontalalignment="right", verticalalignment="bottom")

        ax_last = ax

    ax_last.legend(loc="upper right", bbox_to_anchor=(1, -0.2), borderaxespad=0, ncol=2)

    img_file = img_folder.joinpath("bfe_hydrographs.eps")
    with img_file.open("wb") as f:
        fig.tight_layout()
        fig.savefig(f, bbox_inches="tight", format="eps")

def plot_area_avg(data_dict: dict, bias_dict:dict, img_dir: Path, obs_label_hint="DAYMET",
                  panel_titles=()):

    img_dir.mkdir(parents=True, exist_ok=True)

    # calculate are avg
    ar_avgs = {}
    ar_avg_bias = {}

    times = None
    for data_key, data in data_dict.items():
        if times is None:
            times = data.coords["t"].values
            times = dates.date2num(pd.to_datetime(times).to_pydatetime())

        good_i, good_j = np.where(~data[0].to_masked_array().mask)
        ar_avgs[data_key] = data.values[:, good_i, good_j].mean(axis=1)
        if obs_label_hint not in data_key:
            ar_avg_bias[data_key] = bias_dict[data_key].values[:, good_i, good_j].mean(axis=1)

    # plotting
    plot_utils.apply_plot_params(font_size=14)

    fig = plt.figure()

    ax_list = []
    gs = GridSpec(2, 1)
    ax = fig.add_subplot(gs[0, 0])
    for data_key, data in ar_avgs.items():
        ax.plot(times, data, label=data_key.split("_ndrw")[0])


    ax_list.append(ax)

    if len(panel_titles) > 0:
        ax.set_title(panel_titles[0])

    ax = fig.add_subplot(gs[1, 0])
    for data_key, data in ar_avgs.items():
        if obs_label_hint in data_key:
            continue
        ax.plot(times, ar_avg_bias[data_key], label="$\Delta$" + data_key.split("_ndrw")[0])

    ax_list.append(ax)

    for pt, ax in zip(panel_titles, ax_list):
        ax.set_title(pt)

    for ax in ax_list:
        ax.xaxis.set_major_formatter(DateFormatter("%b"))
        ax.xaxis.set_major_locator(MonthLocator(bymonthday=1))
        ax.legend()

    imfile = img_dir / ("_".join([dl for dl in data_dict if obs_label_hint not in dl]) + ".png")
    fig.savefig(str(imfile), dpi=400, bbox_inches="tight")
    plt.close(fig)

def main_for_cc_paper(start_date=None, end_date=None):

    # Station ids to get from the CEHQ database
    ids_with_lakes_upstream = [
        "104001", "093806", "093801", "081002", "081007", "080718"
    ]

    # for a presentation
    # ids_with_lakes_upstream = [
    #     "104001", "093806", "081002", "081007",
    # ]



    selected_ids = ids_with_lakes_upstream

    sim_labels = [
        # "CRCM5-R",
        # "CRCM5-L2",
        "ERAI-CRCM5-L",
        "CanESM2-CRCM5-L"
        # "CRCM5-HCD-RL-INTF-a"
    ]

    sim_file_names = [
        "/RESCUE/skynet3_rech1/huziy/hdf_store/quebec_0.1_crcm5-hcd-rl.hdf5",
        "/RESCUE/skynet3_rech1/huziy/hdf_store/cc-canesm2-driven/quebec_0.1_crcm5-hcd-rl-cc-canesm2-1980-2010.hdf5"
    ]

    color_list = ["dodgerblue", "r", "g"]
    sim_name_to_color = OrderedDict([
        ("Obs.", "k")
    ])

    for sim_name, the_color in zip(sim_labels, color_list):
        sim_name_to_color[sim_name] = the_color

    sim_name_to_file_name = OrderedDict(zip(sim_labels, sim_file_names))

    # Get the list of stations to do the comparison with
    stations = cehq_station.read_station_data(
        start_date=start_date, end_date=end_date, selected_ids=selected_ids
    )

    print("Initial list of stations:")
    for s in stations:
        print("{0}".format(s))

    plot_utils.apply_plot_params(font_size=16, width_pt=None, width_cm=25, height_cm=12)
    draw_model_comparison(model_points=None, sim_name_to_file_name=sim_name_to_file_name,
                          hdf_folder=None,
                          start_year=start_date.year, end_year=end_date.year, stations=stations,
                          stfl_name="STFL",
                          drainage_area_reldiff_min=0.1,
                          plot_upstream_area_averaged=False,
                          sim_name_to_color=sim_name_to_color)

def compare_obs_and_model_lake_levels_at_points():
    from crcm5.analyse_hdf import compare_lake_levels_with_obs

    plot_utils.apply_plot_params(font_size=14, width_pt=None, width_cm=25, height_cm=18)

    start_date = datetime(1980, 1, 1)
    end_date = datetime(2010, 12, 31)

    compare_lake_levels_with_obs.main(start_date=start_date, end_date=end_date)
    configure()

def validate_forcing_temp_and_precip_with_cru():
    """
    Compare DFS temperature (K) and precip, with the ones from the CRU data set
    """
    from .forcing import compare_forcing_clim_with_CRU
    plot_utils.apply_plot_params(font_size=18, width_pt=None, width_cm=13, height_cm=25)
    compare_forcing_clim_with_CRU.main()
    compare_forcing_clim_with_CRU.main(dfs_var_name="precip",
                                       cru_var_name="pre",
                                       cru_file="data/cru_data/CRUTS3.1/cru_ts_3_10.1901.2009.pre.dat.nc")