"""
##############################################################################################
Project : CRESCENDO/AEROCOM
Filename : check.py
Author : Ramiro Checa-Garcia
email : rcheca@lsce.ipsl.fr
Purpose : Specific type of checks for key variables
Revision History ----------------------------------------------------------
Date Author Ref Revision
2018-Oct R.Checa First version.
2018-Nov R.Checa Added the global mean at level
TODO LIST:
##############################################################################################
"""
import xarray as xr
import pandas as pd
import numpy as np
[docs]def myprint(mystr, finfo=None):
print(mystr)
if finfo!=None:
finfo.write(mystr+'\n')
return
[docs]def total_load(ncname, varname, dicinfo, area_fname='data/area_grid.nc', area_varid='area', finfo=None):
"""
"""
data = xr.open_dataset(ncname)
areadata = xr.open_dataset(area_fname)[area_varid]
'''
try:
print(area_fname, area_varid)
areadata = xr.open_dataset(area_fname)[area_varid]
except:
import iris
variable_constraint = iris.Constraint(cube_func=(lambda c: c.var_name == varname))
cube = iris.load(fn_process, constraints=variable_constraint)[0]
try:
cube.coord('latitude').guess_bounds()
cube.coord('longitude').guess_bounds()
except ValueError:
pass
cube_area = iris.analysis.cartography.area_weights(cube)
areadata=cube_area[0].data
'''
load_field = data[varname]*areadata
load_month = load_field.sum(['lat','lon'])
load_out = load_month.mean(dim='time')/1.e9
myprint(' ... [checking values] [Tg]: '+ str(load_out.values), finfo=finfo)
return varname, 'Tg', load_out.values
def _total_tendency(tendency_data, varname, area):
"""
"""
tendency_field = tendency_data[varname]*area
tendency_month = tendency_field.sum(['lat','lon'])
seconds_month = np.array([k * 86400. for k in [31, 28.25, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]])
tendency_total= np.array([tendency_month[imonth]*seconds_month[imonth] for imonth in range(12)])
acc_tot = np.round(tendency_total.sum()/1.e9,3)
acc_mon = np.round(tendency_total/1.e9,1)
return acc_tot, acc_mon
[docs]def tendency_mass(ncname, varname, area_fname, area_varid, finfo=None):
data = xr.open_dataset(ncname)
area = xr.open_dataset(area_fname)[area_varid]
new = data.groupby('time.month')
new1= new.mean(dim='time')
new1.coords['month']=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
val_tendency, monthly = _total_tendency(new1, varname, area)
myprint(' ... [checking values] [Tg yr-1]: '+ str(val_tendency), finfo=finfo)
myprint(' ... [checking values] monthly: '+ ', '.join([str(y) for y in monthly]), finfo=finfo)
return varname, 'Tg yr-1', val_tendency
[docs]def global_mean(ncname, varname, dicinfo, finfo=None, level=None):
dg_to_rd = np.pi/180.0001
data = xr.open_dataset(ncname) # time, lat, lon
if level==None:
mymethod='global_mean'
datavar = data[varname].values
vlevel=''
else:
mymethod='global_mean_atlev'
#print(data[varname])
datavar = data[varname].isel(pres=level).values
myprint(' evaluating at the pressure ... '+str(data['pres'].values[level])+
' '+str(data['pres'].attrs['units']), finfo=finfo)
vlevel = ' '+str(data['pres'].values[level])+' '+str(data['pres'].attrs['units'])
lat = data.coords['lat'].values
lon = data.coords['lon'].values
time= data.coords['time'].values
lat_nh = np.where(lat>0)
lat_sh = np.where(lat<0)
lataxis=0
lonaxis=1
ts_global=[]
ts_nhemis=[]
ts_shemis=[]
for itime in range(0,len(time)):
data2Dgb = np.squeeze(datavar[itime,: ,:])
data2Dnh = np.squeeze(datavar[itime,lat_nh,:])
data2Dsh = np.squeeze(datavar[itime,lat_sh,:])
data2Davggb = np.average(np.average(data2Dgb, axis=lonaxis),
axis=lataxis, weights=np.cos(lat*dg_to_rd))
data2Davgnh = np.average(np.average(data2Dnh, axis=lonaxis),
axis=lataxis, weights=np.cos(lat[lat_nh]*dg_to_rd))
data2Davgsh = np.average(np.average(data2Dsh, axis=lonaxis),
axis=lataxis, weights=np.cos(lat[lat_sh]*dg_to_rd))
ts_global.append(data2Davggb)
ts_nhemis.append(data2Davgnh)
ts_shemis.append(data2Davgsh)
a_gl = np.array(ts_global) #; a_gl_avg = str(np.round(a_gl.mean(),4))
a_nh = np.array(ts_nhemis) #; a_nh_avg = str(np.round(a_nh.mean(),4))
a_sh = np.array(ts_shemis) # ; a_sh_avg = str(np.round(a_sh.mean(),4))
#a_gl_min = "{0:0.3f}".format(a_gl.min())
#a_gl_max = "{0:0.3f}".format(a_gl.max())
#a_gl_std = "{0:0.3f}".format(a_gl.std())
a_gl1 = "> 90S-90N {0:>8.5g} | NH {1:>8.5g} | SH {2:>8.5g} ".format(a_gl.mean(), a_nh.mean(), a_sh.mean())
a_gl2 = "| min {0:>8.5g} | max {1:>8.5g} | std {2:>8.5g} ".format(a_gl.min(), a_gl.max(), a_gl.std())
myprint(' ... [checking values]: '+
'['+str(data[varname].attrs['units']).rjust(20)+']'+a_gl1+a_gl2, finfo=finfo)
if float(dicinfo['method'][mymethod]['factor'])!=1:
na_gl = a_gl*float(dicinfo['method'][mymethod]['factor'])
na_nh = a_nh*float(dicinfo['method'][mymethod]['factor'])
na_sh = a_sh*float(dicinfo['method'][mymethod]['factor'])
na_gl1 = "> 90S-90N {0:>8.5g} | NH {1:>8.5g} | SH {2:>8.5g} ".format(na_gl.mean(), na_nh.mean(), na_sh.mean())
na_gl2 = "| min {0:>8.5g} | max {1:>8.5g} | std {2:>8.5g} ".format(na_gl.min(), na_gl.max(), na_gl.std() )
myprint(' ... [checking values]: '+
'['+str(dicinfo['method'][mymethod]['units']).rjust(20)+']'
+na_gl1+na_gl2, finfo=finfo)
else:
na_gl= a_gl
return varname, dicinfo['method'][mymethod]['units'], na_gl.mean(), vlevel
[docs]def global_mean_atlev(ncname, varname, dicinfo, level=0, finfo=None):
return global_mean(ncname, varname, dicinfo, finfo=finfo, level=level)