plot_bias_vs_ctrl_veri.py

import numpy as np
import pandas as pd
import dacepy as dp
import xarray as xr
from netCDF4 import Dataset

import re
import glob
import os
import os.path
import sys

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.crs as ccrs

#exp = sys.argv[1]
exp = "ctrl"  # fc06 24h_05 ctrl

save_to='/hpc/uwork/extjroqu/bacy_plots/p01_postprocessing/plots/eva/'

## reading the data##
dr_extra = "/hpc/uhome/extjroqu/bacy_data/"
extp = 'extpardir/icon_extpar_9999_R13B07_L_20200527_tiles.nc'
dext = xr.open_dataset(dr_extra+extp)

irri_crop = dext['LU_CLASS_FRACTION'][0,:]
gridfile = Dataset(dr_extra+"griddir/icon_grid_9999_R13B07_L.nc") #3km eu domain
clon, clat  = np.rad2deg( gridfile.variables["clon"]) , np.rad2deg( gridfile.variables["clat"]) #[::-1] )
indices_above_zero = np.where(irri_crop > 0)
i_clon = clon[indices_above_zero]
i_clat = clat[indices_above_zero]
i_ext = irri_crop[indices_above_zero]

# Define a custom colormap that transitions from green to white to red
cdict = {
    'red':   ((0.0, 0.0, 0.0),
              (0.5, 1.0, 1.0),
              (1.0, 1.0, 1.0)),

    'green': ((0.0, 1.0, 1.0),
              (0.5, 1.0, 1.0),
              (1.0, 0.0, 0.0)),

    'blue':  ((0.0, 0.0, 0.0),
              (0.5, 1.0, 1.0),
              (1.0, 0.0, 0.0))
}

# Create the colormap using the dictionary
red_green_cmap = mcolors.LinearSegmentedColormap('RedGreen', cdict)

def basic_plot(forecast, dataset_fc_fg, clon, clat, irri_crop,title,m_time):
    labsize=16
    projection = ccrs.PlateCarree()
    fig, ax = plt.subplots(figsize=(8.7, 7.5), subplot_kw=dict(projection=projection))
    cmap = "Blues" # plt.cm.Greys
    #Normal plot again
    levels = np.linspace(0,1,10)
    cnf = ax.tricontourf(clon, clat, irri_crop,
                     levels=levels,
                     cmap=cmap,
                     extend='both',
                     #norm=norm,
                     zorder=0)
    extent=[min_lon_s,max_lon_s,min_lat_s,max_lat_s]
    ax.set_extent(extent)
    if m_time == "June" or m_time == "May":
        sc = ax.scatter(forecast.lon, forecast.lat , s=16, c=( abs(forecast['result'].values) - abs(dataset_fc_fg['result'].values)) / abs(dataset_fc_fg['result'].values),
                cmap=red_green_cmap, vmin=-1, vmax=1)
    else:
        sc = ax.scatter(forecast.lon, forecast.lat , s=1, c=( (forecast.values) - (dataset_fc_fg.values)) / (dataset_fc_fg.values), cmap=plt.cm.bwr, vmin=-1, vmax=1)
    gl = ax.gridlines(draw_labels=True,
                 linewidth=0.5,
                 color='dimgray',
                 alpha=0.4,
                 zorder=2)
    gl.xlabel_style = {'size': labsize, 'color': 'black'}
    gl.ylabel_style = {'size': labsize, 'color': 'black'}
    gl.top_labels = False
    gl.right_labels = False
    ax.coastlines(linewidth=0.4, zorder=2)
    plt.colorbar(sc,ax=ax,shrink=0.5,pad=0.03)
    plt.tight_layout()
    plt.savefig(os.path.join(save_to,"rd_"+exp+"_"+country+"_"+title+"_"+m_time+".png"),dpi=300, bbox_inches='tight')
    plt.close()

def only_plot(lon,lat,data,dtype):
    labsize=16
    projection = ccrs.PlateCarree()
    fig, ax = plt.subplots(figsize=(8.7, 7.5), subplot_kw=dict(projection=projection))
    cmap = "Blues" # plt.cm.Greys
    #Normal plot again
    levels = np.linspace(0,1,10)
    cnf = ax.tricontourf(clon, clat, irri_crop,
                     levels=levels,
                     cmap=cmap,
                     extend='both',
                     #norm=norm,
                     zorder=0)
    ax.set_extent([min_lon_s,max_lon_s,min_lat_s,max_lat_s], crs=projection)
    sc = ax.scatter(lon['lon'],lat['lat'], s=16, c=(data['result']), cmap=red_green_cmap, vmin=-1, vmax=1)
    gl = ax.gridlines(draw_labels=True,
                 linewidth=0.5,
                 color='dimgray',
                 alpha=0.4,
                 zorder=2)
    gl.xlabel_style = {'size': labsize, 'color': 'black'}
    gl.ylabel_style = {'size': labsize, 'color': 'black'}
    gl.top_labels = False
    gl.right_labels = False
    ax.coastlines(linewidth=0.4, zorder=2)
    plt.colorbar(sc,ax=ax, label='Bias',shrink=0.5,pad=0.03)
    plt.tight_layout()
    plt.savefig(os.path.join(save_to,"bias_"+exp+"_"+country+"_"+str(data.month.values.item())+"_"+dtype+".png"),dpi=300, bbox_inches='tight')
    plt.close()

country_list = ["Spain","Italy","Greece"]

for i in country_list:
    country = i

    if country == "Spain":
        from_path = '/hpc/uwork/extjroqu/reanalyses/evaluation/veri_'+exp+'_nwp0_spain/'
        from_path_c = '/hpc/uwork/extjroqu/reanalyses/evaluation/veri_ctrl_nwp0_spain/'
    elif country == "Italy":
        from_path = '/hpc/uwork/extjroqu/reanalyses/evaluation/veri_'+exp+'_nwp0_italy/'
        from_path_c = '/hpc/uwork/extjroqu/reanalyses/evaluation/veri_ctrl_nwp0_italy/'
    else:
        from_path ='/hpc/uwork/extjroqu/reanalyses/evaluation/veri_'+exp+'_nwp0_greece/'
        from_path_c = '/hpc/uwork/extjroqu/reanalyses/evaluation/veri_ctrl_nwp0_greece/'

    # Reading data
    fc_d = xr.open_dataset(from_path+"test_fc_mon_v39_"+country+"_june.nc")
    fg_d =  xr.open_dataset(from_path_c+"test_fc_mon_v39_"+country+"_june.nc")

    fc_d_all = xr.open_dataset(from_path+"veri_all/test_fc_mon_v39_"+country+"_june_all.nc")
    fg_d_all =  xr.open_dataset(from_path_c+"veri_all/test_fc_mon_v39_"+country+"_june_all.nc")

    if country == 'Spain':
        min_lon_s = -10 #-10
        max_lon_s = 2 #4.5
        min_lat_s = 35.5
        max_lat_s = 44
    elif country == 'Italy':
        min_lon_s = 6.5 #-10
        max_lon_s = 16.5 #4.5
        min_lat_s = 37.5
        max_lat_s = 47.5
    else: #country = 'Greece'
        min_lon_s = 19 #-10
        max_lon_s = 29 #4.5
        min_lat_s = 34
        max_lat_s = 44

##

    basic_plot(fc_d,fg_d,clon, clat, irri_crop,"FC_irr","June") # irri
    basic_plot(fc_d_all,fg_d_all,clon, clat, irri_crop,"FC_all","June") # all

    only_plot(fc_d,fc_d,fc_d,"FC")

    print("DONE with "+exp+" and "+country)