From bf25fd4d7d87eeb10d13ba0d81dd794c98fe96e3 Mon Sep 17 00:00:00 2001 From: Ruud van der Ent Date: Thu, 7 Jul 2016 16:25:32 +0200 Subject: [PATCH] see readme --- Con_E_Recyc_Output.py | 4 +- Con_E_Recyc_Plots.py | 13 ++-- Con_P_Recyc_Plots.py | 129 ++++++++++++++++++++++------------------ Hor_Fluxes_Output.py | 134 ++++++++++++++++++++++++++++++++++++++++++ readme.txt | 4 ++ 5 files changed, 221 insertions(+), 63 deletions(-) create mode 100644 Hor_Fluxes_Output.py diff --git a/Con_E_Recyc_Output.py b/Con_E_Recyc_Output.py index 6b0931c..6a1ed47 100644 --- a/Con_E_Recyc_Output.py +++ b/Con_E_Recyc_Output.py @@ -5,8 +5,8 @@ @author: Ent00002 """ #delayed runs -#import time -#time.sleep(55004) +import time +time.sleep(44000) #%% Import libraries diff --git a/Con_E_Recyc_Plots.py b/Con_E_Recyc_Plots.py index 100d148..078dfdf 100644 --- a/Con_E_Recyc_Plots.py +++ b/Con_E_Recyc_Plots.py @@ -205,8 +205,8 @@ def data_path(years,timetracking): map.drawparallels(np.arange(-90,90,30),linewidth=0.1) x_shift,y_shift = np.meshgrid(longitude,latitude) -x = x_shift - 0.75 -y = y_shift - 0.75 +x = x_shift +y = y_shift # contour data over the map. clevs = np.linspace(0,1.00,11) @@ -227,7 +227,7 @@ def data_path(years,timetracking): x_shift,y_shift = np.meshgrid(longitude,latitude) x = x_shift - 0.75 -y = y_shift - 0.75 +y = y_shift + 0.75 lol = map.pcolormesh(x,y,eps_c_landtotal, latlon=True, cmap=cmap,vmin=0,vmax=1) map.colorbar(lol,location='right',pad="5%", label = '(-)') @@ -243,7 +243,12 @@ def data_path(years,timetracking): plt.imshow(eps_c_landtotal) plt.colorbar() -#%% Timetracking results +#%% Timetracking results +# these results are very crude and just meant for a quick check. For publications these have been calculated more detailed + +####################################### +# possible check to do the same with daily output or timestep output +######################################## # for some reasons the sizes of the variables in the "variable explorer" of Spyder get wrongly displayed after running code below, but calling them from the memory works fine if timetracking == 1: diff --git a/Con_P_Recyc_Plots.py b/Con_P_Recyc_Plots.py index 3d51683..2ee8f2b 100644 --- a/Con_P_Recyc_Plots.py +++ b/Con_P_Recyc_Plots.py @@ -76,46 +76,46 @@ def data_path(years,timetracking): P_total = np.squeeze(np.sum(np.sum(P_per_year_per_month,axis = 0), axis = 0)) # loading horizontal fluxes -#loading_HF2S = sio.loadmat(datapath[1]) -#Fa_E_down_per_year_per_month = loading_HF2S['Fa_E_down_per_year_per_month'] -#Fa_E_top_per_year_per_month = loading_HF2S['Fa_E_top_per_year_per_month'] -#Fa_N_down_per_year_per_month = loading_HF2S['Fa_N_down_per_year_per_month'] -#Fa_N_top_per_year_per_month = loading_HF2S['Fa_N_top_per_year_per_month'] -#Fa_Vert_per_year_per_month = loading_HF2S['Fa_Vert_per_year_per_month'] -# -#Fa_E_down_total = np.squeeze(np.sum(np.sum(Fa_E_down_per_year_per_month,axis = 0),axis = 0)) -#Fa_E_top_total = np.squeeze(np.sum(np.sum(Fa_E_top_per_year_per_month,axis = 0),axis = 0)) -#Fa_N_down_total = np.squeeze(np.sum(np.sum(Fa_N_down_per_year_per_month,axis = 0),axis = 0)) -#Fa_N_top_total = np.squeeze(np.sum(np.sum(Fa_N_top_per_year_per_month,axis = 0),axis = 0)) -#Fa_E_total = Fa_E_down_total + Fa_E_top_total -#Fa_N_total = Fa_N_down_total + Fa_N_top_total -# -## annual average fluxes per unit width (m3/m/a) -#Fa_N_total_m2a = np.zeros(np.shape(Fa_N_total)) -#Fa_N_down_total_m2a = np.zeros(np.shape(Fa_N_down_total)) -#Fa_N_top_total_m2a = np.zeros(np.shape(Fa_N_top_total)) -# -#for i in range(len(latitude)): -# Fa_N_total_m2a[i,:] = (Fa_N_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) -# Fa_N_down_total_m2a[i,:] = (Fa_N_down_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) -# Fa_N_top_total_m2a[i,:] = (Fa_N_top_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) -# -#Fa_E_total_m2a = (Fa_E_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) -#Fa_E_down_total_m2a = (Fa_E_down_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) -#Fa_E_top_total_m2a = (Fa_E_top_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) -# -## fluxes for plotting -#FluxDispMatSmall = np.zeros([5,5]) -#FluxDispMatSmall[3,2] = 1 -#FluxDispMat = np.tile(FluxDispMatSmall,[np.floor(len(latitude)/5+1),np.floor(len(longitude)/5)+1]) -#FluxDispMat = FluxDispMat[:len(latitude),:len(longitude)] -# -#Fa_E_total_m2a_part = Fa_E_total_m2a * FluxDispMat -#Fa_N_total_m2a_part = Fa_N_total_m2a * FluxDispMat -#Fa_E_down_total_m2a_part = Fa_E_down_total_m2a * FluxDispMat -#Fa_N_down_total_m2a_part = Fa_N_down_total_m2a * FluxDispMat -#Fa_E_top_total_m2a_part = Fa_E_top_total_m2a * FluxDispMat -#Fa_N_top_total_m2a_part = Fa_N_top_total_m2a * FluxDispMat +loading_HF2S = sio.loadmat(datapath[1]) +Fa_E_down_per_year_per_month = loading_HF2S['Fa_E_down_per_year_per_month'] +Fa_E_top_per_year_per_month = loading_HF2S['Fa_E_top_per_year_per_month'] +Fa_N_down_per_year_per_month = loading_HF2S['Fa_N_down_per_year_per_month'] +Fa_N_top_per_year_per_month = loading_HF2S['Fa_N_top_per_year_per_month'] +Fa_Vert_per_year_per_month = loading_HF2S['Fa_Vert_per_year_per_month'] + +Fa_E_down_total = np.squeeze(np.sum(np.sum(Fa_E_down_per_year_per_month,axis = 0),axis = 0)) +Fa_E_top_total = np.squeeze(np.sum(np.sum(Fa_E_top_per_year_per_month,axis = 0),axis = 0)) +Fa_N_down_total = np.squeeze(np.sum(np.sum(Fa_N_down_per_year_per_month,axis = 0),axis = 0)) +Fa_N_top_total = np.squeeze(np.sum(np.sum(Fa_N_top_per_year_per_month,axis = 0),axis = 0)) +Fa_E_total = Fa_E_down_total + Fa_E_top_total +Fa_N_total = Fa_N_down_total + Fa_N_top_total + +# annual average fluxes per unit width (m3/m/a) +Fa_N_total_m2a = np.zeros(np.shape(Fa_N_total)) +Fa_N_down_total_m2a = np.zeros(np.shape(Fa_N_down_total)) +Fa_N_top_total_m2a = np.zeros(np.shape(Fa_N_top_total)) + +for i in range(len(latitude)): + Fa_N_total_m2a[i,:] = (Fa_N_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) + Fa_N_down_total_m2a[i,:] = (Fa_N_down_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) + Fa_N_top_total_m2a[i,:] = (Fa_N_top_total[i,:] / (0.5*(L_N_gridcell[i] + L_S_gridcell[i])))/len(years) + +Fa_E_total_m2a = (Fa_E_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) +Fa_E_down_total_m2a = (Fa_E_down_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) +Fa_E_top_total_m2a = (Fa_E_top_total / L_EW_gridcell) / len(years) # Annual average flux per unit width (m3/m/a) + +# fluxes for plotting +FluxDispMatSmall = np.zeros([5,5]) +FluxDispMatSmall[3,2] = 1 +FluxDispMat = np.tile(FluxDispMatSmall,[np.floor(len(latitude)/5+1),np.floor(len(longitude)/5)+1]) +FluxDispMat = FluxDispMat[:len(latitude),:len(longitude)] + +Fa_E_total_m2a_part = Fa_E_total_m2a * FluxDispMat +Fa_N_total_m2a_part = Fa_N_total_m2a * FluxDispMat +Fa_E_down_total_m2a_part = Fa_E_down_total_m2a * FluxDispMat +Fa_N_down_total_m2a_part = Fa_N_down_total_m2a * FluxDispMat +Fa_E_top_total_m2a_part = Fa_E_top_total_m2a * FluxDispMat +Fa_N_top_total_m2a_part = Fa_N_top_total_m2a * FluxDispMat # recycling metrics rho_c_total = P_c_total / P_total @@ -200,23 +200,18 @@ def data_path(years,timetracking): map.drawparallels(np.arange(-90,90,30),linewidth=0.1) x_shift,y_shift = np.meshgrid(longitude,latitude) -x = x_shift - 0.75 -y = y_shift - 0.75 +x = x_shift +y = y_shift # contour data over the map -#Fa_E_total_m2a_part[Fa_E_total_m2a_part==0] = np.nan -#Fa_N_total_m2a_part[Fa_N_total_m2a_part==0] = np.nan +Fa_E_total_m2a_part[Fa_E_total_m2a_part==0] = np.nan +Fa_N_total_m2a_part[Fa_N_total_m2a_part==0] = np.nan clevs = np.linspace(0,1.00,11) clevs2 = np.linspace(0,1,1100) -cs = map.contourf(x_shift,y_shift,rho_c_landtotal,clevs2,linewidths=1.5,latlon=True, cmap=cmap) +cs = map.contourf(x,y,rho_c_landtotal,clevs2,linewidths=1.5,latlon=True, cmap=cmap) map.colorbar(cs,location='right',pad="5%", ticks = clevs, label = '(-)') -# standard way of plotting (looks worse) -#cs = map.contourf(x,y,rho_c_landtotal,linewidths=1.5,latlon=True) -#map.colorbar(cs,location='right', label = '(-)') - - -#windplot = map.quiver(x,y,Fa_E_total_m2a_part,Fa_N_total_m2a_part,color = 'black',scale = 199999999,latlon=True,width = 0.002, headaxislength = 5, headwidth = 2) +windplot = map.quiver(x,y,Fa_E_total_m2a_part,Fa_N_total_m2a_part,color = 'black',scale = 199999999,latlon=True,width = 0.002, headaxislength = 5, headwidth = 2) #adjust scale till arrows are visible plt.title('Continental precipitation recycling ratio ' r'$\rho_c$', fontsize = 14) @@ -238,7 +233,7 @@ def data_path(years,timetracking): lol = map.pcolormesh(x,y,rho_c_landtotal, latlon=True, cmap=cmap,vmin=0,vmax=1) map.colorbar(lol,location='right',pad="5%", label = '(-)') -#windplot = map.quiver(x,y,Fa_E_total_m2a_part,Fa_N_total_m2a_part,color = 'black',scale = 199999999,latlon=True,width = 0.002, headaxislength = 5, headwidth = 2) +windplot = map.quiver(x,y,Fa_E_total_m2a_part,Fa_N_total_m2a_part,color = 'black',scale = 199999999,latlon=True,width = 0.002, headaxislength = 5, headwidth = 2) #adjust scale till arrows are visible plt.title('Continental precipitation recycling ratio ' r'$\rho_c$', fontsize = 14) @@ -250,7 +245,12 @@ def data_path(years,timetracking): plt.imshow(rho_c_landtotal) plt.colorbar -#%% Timetracking results +#%% Timetracking results +# these results are very crude and just meant for a quick check. For publications these have been calculated more detailed + +####################################### +# possible check to do the same with daily output or timestep output +######################################## if timetracking == 1: Sa_track_down_average = np.squeeze(np.mean(np.mean(Sa_track_down_per_year_per_month, axis=0), axis = 0)) #m3 @@ -264,8 +264,8 @@ def data_path(years,timetracking): Sa_time_top_average = np.squeeze(np.mean(np.mean(Sa_time_top_per_year_per_month, axis=0), axis = 0)) / 86400 #d quick assesment Sa_time_average = (( np.squeeze(np.mean(np.mean(Sa_time_down_per_year_per_month * Sa_track_down_per_year_per_month, axis=0), axis = 0)) + np.squeeze(np.mean(np.mean(Sa_time_top_per_year_per_month * Sa_track_top_per_year_per_month,axis=0), axis = 0)) - ) / Sa_track_average / 86400) #d - + ) / Sa_track_average / 86400) #d + Sa_track_down_average = np.squeeze(np.mean(np.mean(Sa_track_down_per_year_per_month, axis=0), axis = 0)) #m3 Sa_track_top_average = np.squeeze(np.mean(np.mean(Sa_track_top_per_year_per_month, axis=0), axis = 0)) #m3 @@ -273,10 +273,25 @@ def data_path(years,timetracking): / np.sum(np.sum(Sa_track_down_average[1:-1,:], axis = 0), axis = 0)) #d age_top_track_water_global_mean = (np.sum(np.sum(Sa_time_top_average[1:-1,:] * Sa_track_top_average[1:-1,:], axis = 0 ), axis = 0) / np.sum(np.sum(Sa_track_top_average[1:-1,:], axis = 0), axis = 0)) #d - age_track_water_global_mean = (( np.sum(np.sum(Sa_time_down_average[1:-1,:] * Sa_track_down_average[1:-1,:], axis = 0), axis = 0) + + # weighted with the volume of moisture in the atmosphere + age_track_water_global_mean_atmos_weight2 = (( np.sum(np.sum(Sa_time_down_average[1:-1,:] * Sa_track_down_average[1:-1,:], axis = 0), axis = 0) + np.sum(np.sum(Sa_time_top_average[1:-1,:] * Sa_track_top_average[1:-1,:], axis = 0 ), axis = 0) ) / ( np.sum(np.sum(Sa_track_down_average[1:-1,:], axis = 0), axis = 0) + np.sum(np.sum(Sa_track_top_average[1:-1,:], axis = 0), axis = 0) )) #d - + + age_track_water_global_mean_atmos_weight = ( ( np.sum(Sa_time_down_per_year_per_month[:,:,1:-1,:] * Sa_track_down_per_year_per_month[:,:,1:-1,:]) + + np.sum(Sa_time_top_per_year_per_month[:,:,1:-1,:] * Sa_track_top_per_year_per_month[:,:,1:-1,:]) ) + / ( np.sum(Sa_track_down_per_year_per_month[:,:,1:-1,:]) + np.sum(Sa_track_top_per_year_per_month[:,:,1:-1,:]) ) ) / 86400 #d + + # weighted with precipitation + age_track_water_global_mean_precip_weight = ( np.sum( P_time_per_year_per_month[:,:,1:-1,:] * P_track_per_year_per_month[:,:,1:-1,:] ) + / np.sum(P_track_per_year_per_month) / 86400 )# d + + # weigthed with area + age_track_water_global_mean_area_weight = ( np.sum(Sa_time_average[1:-1,:] * A_gridcell2D[1:-1,:]) + / np.sum(A_gridcell2D) )# d + + # spatial figure time of tracked precipitation P_c_time_total = (np.squeeze(np.sum(np.sum(P_time_per_year_per_month * P_track_per_year_per_month, axis =0), axis = 0)) / P_c_total / 86400) #d diff --git a/Hor_Fluxes_Output.py b/Hor_Fluxes_Output.py new file mode 100644 index 0000000..8349666 --- /dev/null +++ b/Hor_Fluxes_Output.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Jun 16 13:24:45 2016 + +@author: Ent00002 +""" + +#delayed runs +#import time +#time.sleep(7500) + +#%% Import libraries + +import matplotlib.pyplot as plt +import numpy as np +from netCDF4 import Dataset +get_ipython().magic(u'matplotlib inline') +import numpy.matlib +import scipy.io as sio +import calendar +import time +import datetime +from getconstants import getconstants +from timeit import default_timer as timer + +#%% BEGIN OF INPUT (FILL THIS IN) +years = np.arange(2002,2009) #fill in the years +yearpart = np.arange(0,366) # for a full (leap)year fill in 0:366 + +# Manage the extent of your dataset (FILL THIS IN) +# Define the latitude and longitude cell numbers to consider and corresponding lakes that should be considered part of the land +latnrs = np.arange(7,114) +lonnrs = np.arange(0,240) + +# the lake numbers below belong to the ERA-Interim data on 1.5 degree starting at Northern latitude 79.5 and longitude -180 +lake_mask_1 = np.array([9,9,9,12,12,21,21,22,22,23,24,25,23,23,25,25,53,54,61,23,24,23,24,25,27,22,23,24,25,26,27,28,22,25,26,27,28,23,23,12,18]) +lake_mask_2 = np.array([120+19,120+40,120+41,120+43,120+44,120+61,120+62,120+62,120+63,120+62,120+62,120+62,120+65,120+66,120+65,120+66,142-120,142-120,143-120,152-120,152-120,153-120,153-120,153-120,153-120,154-120,154-120,154-120,154-120,154-120,154-120,154-120,155-120,155-120,155-120,155-120,155-120,159-120,160-120,144-120,120+55]) +lake_mask = np.transpose(np.vstack((lake_mask_1,lake_mask_2))) #recreate the arrays of the matlab model + +daily = 1 # 1 for writing out daily data, 0 for only monthly data + +#END OF INPUT + +#%% Datapaths (FILL THIS IN) +invariant_data = 'Interim_data/full/invariants.nc'#invariants + +def data_path(y,a,years): + load_fluxes_and_storages = 'interdata/' + str(y) + '-' + str(a) + 'fluxes_storages.mat' + + save_path = 'outputdata/Hor_Fluxes_full' + str(years[0]) + '-' + str(years[-1]) + '.mat' + + save_path_daily = 'outputdata/Hor_Fluxes_daily_full' + str(y) + '.mat' + + return load_fluxes_and_storages,save_path,save_path_daily + +#%% Runtime % Results + +start1 = timer() + +# obtain the constants +latitude,longitude,lsm,g,density_water,timestep,A_gridcell,L_N_gridcell,L_S_gridcell,L_EW_gridcell,gridcell = getconstants(latnrs,lonnrs,lake_mask,Dataset,invariant_data,np) + +startyear = years[0] +Fa_E_down_per_year_per_month = np.zeros((len(years),12,len(latitude),len(longitude))) +Fa_E_top_per_year_per_month = np.zeros((len(years),12,len(latitude),len(longitude))) +Fa_N_down_per_year_per_month = np.zeros((len(years),12,len(latitude),len(longitude))) +Fa_N_top_per_year_per_month = np.zeros((len(years),12,len(latitude),len(longitude))) +Fa_Vert_per_year_per_month = np.zeros((len(years),12,len(latitude),len(longitude))) + +for i in range(len(years)): + y = years[i] + ly = int(calendar.isleap(y)) + final_time = 364+ly + + Fa_E_down_per_day = np.zeros((365+ly,len(latitude),len(longitude))) + Fa_E_top_per_day = np.zeros((365+ly,len(latitude),len(longitude))) + Fa_N_down_per_day = np.zeros((365+ly,len(latitude),len(longitude))) + Fa_N_top_per_day = np.zeros((365+ly,len(latitude),len(longitude))) + Fa_Vert_per_day = np.zeros((365+ly,len(latitude),len(longitude))) + + for j in range(len(yearpart)): + start = timer() + a = yearpart[j] + datapath = data_path(y,a,years) + + if a > final_time: # a = 365 and not a leapyear + pass + else: + # load horizontal fluxes + loading_FS = sio.loadmat(datapath[0],verify_compressed_data_integrity=False) + Fa_E_top = loading_FS['Fa_E_top'] + Fa_N_top = loading_FS['Fa_N_top'] + Fa_E_down = loading_FS['Fa_E_down'] + Fa_N_down = loading_FS['Fa_N_down'] + Fa_Vert = loading_FS['Fa_Vert'] + + # save per day + Fa_E_down_per_day[a,:,:] = np.sum(Fa_E_down, axis =0) + Fa_E_top_per_day[a,:,:] = np.sum(Fa_E_top, axis =0) + Fa_N_down_per_day[a,:,:] = np.sum(Fa_N_down, axis =0) + Fa_N_top_per_day[a,:,:] = np.sum(Fa_N_top, axis =0) + Fa_Vert_per_day[a,:,:] = np.sum(Fa_Vert, axis =0) + + # timer + end = timer() + print 'Runtime output for day ' + str(a+1) + ' in year ' + str(y) + ' is',(end - start),' seconds.' + + # save daily fluxes on disk + if daily == 1: + sio.savemat(datapath[2], + {'Fa_E_down_per_day':Fa_E_down_per_day,'Fa_E_top_per_day':Fa_E_top_per_day, + 'Fa_N_down_per_day':Fa_N_down_per_day,'Fa_N_top_per_day':Fa_N_top_per_day, + 'Fa_Vert_per_day':Fa_Vert_per_day},do_compression=True) + + for m in range(12): + first_day = int(datetime.date(y,m+1,1).strftime("%j")) + last_day = int(datetime.date(y,m+1,calendar.monthrange(y,m+1)[1]).strftime("%j")) + days = np.arange(first_day,last_day+1)-1 # -1 because Python is zero-based + + Fa_E_down_per_year_per_month[y-startyear,m,:,:] = (np.squeeze(np.sum(Fa_E_down_per_day[days,:,:], axis = 0))) + Fa_E_top_per_year_per_month[y-startyear,m,:,:] = (np.squeeze(np.sum(Fa_E_top_per_day[days,:,:], axis = 0))) + Fa_N_down_per_year_per_month[y-startyear,m,:,:] = (np.squeeze(np.sum(Fa_N_down_per_day[days,:,:], axis = 0))) + Fa_N_top_per_year_per_month[y-startyear,m,:,:] = (np.squeeze(np.sum(Fa_N_top_per_day[days,:,:], axis = 0))) + Fa_Vert_per_year_per_month[y-startyear,m,:,:] = (np.squeeze(np.sum(Fa_Vert_per_day[days,:,:], axis = 0))) + +# save monthly fluxes on disk +sio.savemat(datapath[1], + {'Fa_E_down_per_year_per_month':Fa_E_down_per_year_per_month,'Fa_E_top_per_year_per_month':Fa_E_top_per_year_per_month, + 'Fa_N_down_per_year_per_month':Fa_N_down_per_year_per_month,'Fa_N_top_per_year_per_month':Fa_N_top_per_year_per_month, + 'Fa_Vert_per_year_per_month':Fa_Vert_per_year_per_month}) + +end1 = timer() +print 'The total runtime is',(end1-start1),' seconds.' + diff --git a/readme.txt b/readme.txt index 3004508..c6fbcf5 100644 --- a/readme.txt +++ b/readme.txt @@ -5,6 +5,10 @@ I might be able to offer some kind of support if it is in my interest as well, b Released under the GNU General Public License +WAM-2layers v2.4.07 | 7-7-2016 +- added the Hor_fluxes output file +- small changes to the example plotting files + WAM-2layers v2.4.06 | 6-7-2016 - solve another bug in backward time tracking - more accurate area size calculation