plot_validation_results.py

#!/usr/bin/env python

'''Verifies reanalysis data against AEMET station data, both brought into a comparable and harmonized format using 1. csv2nc.py and 2. get_neighbour.py prior to calling 3. validate.py (i.e. this script).
Author: Swen Brands, brandssf@ifca.unican.es
'''

#load packages
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import os
import pandas as pd
import xskillscore as xs
from math import radians, cos, sin, asin, sqrt #needed to calculate haversine distance
home = os.getenv('HOME')
exec(open('functions_radiation.py').read())
exec(open(home+'/datos/tareas/proyectos/pticlima/seasonal/python/functions_seasonal.py').read())

#set input parameters
model_dataset = ['era5_land','era5'] #list of model or reanalysis datasets: era5 or era5_land
ref_dataset = ['aemet','aemet'] # #list of model or reference observational dataset paired with to <model_dataset>: aemet
rundir = home+'/datos/tareas/proyectos/pticlima/pyPTIclima/pySolar' #script directory, you should be there or point to this directory when running these scripts via python
dir_netcdf = home+'/datos/tareas/proyectos/pticlima/radiation/results/validation' #path to outupt netcdf files produced with this script, containing an xarray dataset with all verification results
dir_figs = home+'/datos/tareas/proyectos/pticlima/radiation/results/validation' #path to output figures file generated by this script
corr_outlier = 'yes' #load the outlier-correted validation results; yes or no

scores = ['pearson_r','spearman_r','bias','relbias','mae','rmse'] #the scores to be plotted

precision = 'float32' #precision of the variable in the output netCDF files
dpival = 300 #resultion of the output figure in dpi
figformat = 'pdf' #format of the output figures: pdf, png, etc.
colormap_ascend = 'Spectral_r' #ascendig colormap used for plotting: Spectral_r 
colormap_div = 'seismic' #diverging (zero-centered) colormap used for plotting: seismic

##EXECUTE ##############################################################
#create output directories ir they do not exist.
if os.path.isdir(dir_figs) != True:
    os.makedirs(dir_figs)

minvals = np.empty((len(model_dataset),len(scores)))
maxvals = np.empty((len(model_dataset),len(scores)))
#get global minimum and maximum values
for dd in range(len(model_dataset)):
    colormaps = [] #the colormap is the same for each dataset and only depends on the score so this is a list with a length equalling the length of <scores>
    #load netcdf files containing the verification results
    if model_dataset[dd] == 'era5_land':
        filename_results = 'verification_results_day_'+model_dataset[dd]+'_vs_'+ref_dataset[dd]+'_corroutlier_'+corr_outlier+'_20171201_20221230.nc'
    elif model_dataset[dd] == 'era5':
        filename_results = 'verification_results_day_'+model_dataset[dd]+'_vs_'+ref_dataset[dd]+'_corroutlier_'+corr_outlier+'_20171201_20221231.nc'
    else:
        raise Exception('ERROR: unknown entry for <model_dataset> !')
    nc_results = xr.open_dataset(dir_netcdf+'/'+filename_results)
    
    #extract the min and max values for each score
    for sc in np.arange(len(scores)):
        if scores[sc] in ('bias','relbias'):
            maxvals[dd,sc] = np.abs(nc_results[scores[sc]]).max().values
            minvals[dd,sc] = np.abs(nc_results[scores[sc]]).max().values*-1
            colormaps.append(colormap_div)
        elif scores[sc] in ('pearson_r','spearman_r'):
            maxvals[dd,sc] = 1
            minvals[dd,sc] = 0
            colormaps.append(colormap_ascend)
        elif scores[sc] in ('mae','mape','rmse'):
            maxvals[dd,sc] = nc_results[scores[sc]].max().values
            minvals[dd,sc] = nc_results[scores[sc]].min().values
            colormaps.append(colormap_ascend)
        else:
            raise Exception('ERROR: unknown value for <scores[sc]> !')

#close file
nc_results.close()
#get global min and max values covering the range of all considered model / reanalysis datasets
minvals = np.min(minvals,axis=0)
maxvals = np.max(maxvals,axis=0)

#then plot the results with this min and max values
for dd in range(len(model_dataset)):
    #load netcdf files containing the verification results
    if model_dataset[dd] == 'era5_land':
        filename_results = 'verification_results_day_'+model_dataset[dd]+'_vs_'+ref_dataset[dd]+'_corroutlier_'+corr_outlier+'_20171201_20221230.nc'
    elif model_dataset[dd] == 'era5':
        filename_results = 'verification_results_day_'+model_dataset[dd]+'_vs_'+ref_dataset[dd]+'_corroutlier_'+corr_outlier+'_20171201_20221231.nc'
    else:
        raise Exception('ERROR: unknown entry for <model_dataset> !')
    nc_results = xr.open_dataset(dir_netcdf+'/'+filename_results)

    ##plot matrices of verification results for the distinct score (x-axis = seasons, y-axis = stations and save to <figformat>
    for sc in np.arange(len(scores)):
        savename = dir_figs+'/'+scores[sc]+'_day_'+model_dataset[dd]+'_vs_aemet_corroutlier_'+corr_outlier+'.'+figformat
        plot_pcolormesh(nc_results,scores[sc],minvals[sc],maxvals[sc],savename,colormaps[sc],dpival)

##close input nc files and produced xr dataset
nc_results.close()
print('INFO: plot_validation_results.py has been run successfully !')