diff --git a/src/lisfloodutilities/water-demand-historic/step2_domestic_demand.py b/src/lisfloodutilities/water-demand-historic/step2_domestic_demand.py index 36f1aeb..a38859a 100644 --- a/src/lisfloodutilities/water-demand-historic/step2_domestic_demand.py +++ b/src/lisfloodutilities/water-demand-historic/step2_domestic_demand.py @@ -401,25 +401,25 @@ def main(): # Plot figure #-------------------------------------------------------------------------- - # Initialize figure - f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) + # # Initialize figure + # f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) - # Subpanel 1 - ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=12,cmap=plt.get_cmap('YlGnBu')) - ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') - - # Subpanel 2 - try: - timeindex = year-1971 - ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[:,:,np.arange(timeindex*12,timeindex*12+12)],axis=2)),vmin=0,vmax=12,cmap=plt.get_cmap('YlGnBu')) - ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') - except: - pass + # # Subpanel 1 + # ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=12,cmap=plt.get_cmap('YlGnBu')) + # ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') + + # # Subpanel 2 + # try: + # timeindex = year-1971 + # ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[:,:,np.arange(timeindex*12,timeindex*12+12)],axis=2)),vmin=0,vmax=12,cmap=plt.get_cmap('YlGnBu')) + # ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') + # except: + # pass - # Save figure - f.set_size_inches(10, 10) - plt.savefig(os.path.join(config['output_folder'],'step2_domestic_demand','figures','dom_'+str(year)+'.png'),dpi=150) - plt.close() + # # Save figure + # f.set_size_inches(10, 10) + # plt.savefig(os.path.join(config['output_folder'],'step2_domestic_demand','figures','dom_'+str(year)+'.png'),dpi=150) + # plt.close() print("Time elapsed is "+str(time.time()-t0)+" sec") diff --git a/src/lisfloodutilities/water-demand-historic/step3a_industrial_demand.py b/src/lisfloodutilities/water-demand-historic/step3a_industrial_demand.py index 6015393..2455de0 100644 --- a/src/lisfloodutilities/water-demand-historic/step3a_industrial_demand.py +++ b/src/lisfloodutilities/water-demand-historic/step3a_industrial_demand.py @@ -18,7 +18,7 @@ import pandas as pd import numpy as np from netCDF4 import Dataset -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt from skimage.transform import resize from tools import * import rasterio @@ -666,39 +666,39 @@ def main(): # Plot figure #-------------------------------------------------------------------------- - # Initialize figure - f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True) + # # Initialize figure + # f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True) - # Subpanel 1 - ax1.plot(years,table_aquastat_industry_withdrawal[jj,:],'*',label='FAO industry withdrawal',color=(0.5,0.5,0.5)) - ax1.plot(years,table_gcam_withdrawal_industry[jj,:],'+',label='GCAM industry withdrawal',color=(0,0.5,0)) - ax1.plot(years,table_gcam_withdrawal_industry[jj,:]-table_gcam_withdrawal_thermoelectric[jj,:],'+',label='GCAM manufacturing withdrawal',color=(0,1,0)) - ax1.plot(years,table_withdrawal_industry[jj,:],label='Beck industry withdrawal',color=(0.5,0,0)) - ax1.plot(years,table_withdrawal_manufacturing[jj,:],label='Beck manufacturing withdrawal',color=(1,0,0)) - ax1.legend() - ax1.set_ylabel('Withdrawal (km^3/year)') - ax1.set_title(country_name) + # # Subpanel 1 + # ax1.plot(years,table_aquastat_industry_withdrawal[jj,:],'*',label='FAO industry withdrawal',color=(0.5,0.5,0.5)) + # ax1.plot(years,table_gcam_withdrawal_industry[jj,:],'+',label='GCAM industry withdrawal',color=(0,0.5,0)) + # ax1.plot(years,table_gcam_withdrawal_industry[jj,:]-table_gcam_withdrawal_thermoelectric[jj,:],'+',label='GCAM manufacturing withdrawal',color=(0,1,0)) + # ax1.plot(years,table_withdrawal_industry[jj,:],label='Beck industry withdrawal',color=(0.5,0,0)) + # ax1.plot(years,table_withdrawal_manufacturing[jj,:],label='Beck manufacturing withdrawal',color=(1,0,0)) + # ax1.legend() + # ax1.set_ylabel('Withdrawal (km^3/year)') + # ax1.set_title(country_name) - # Subpanel 2 - ax2.plot(years,table_Huang_withdrawal_industry[jj,:],label='Huang industry withdrawal',color=(0,0,0.5)) - ax2.plot(years,table_Huang_withdrawal_manufacturing[jj,:],label='Huang manufacturing withdrawal',color=(0,0,1)) - ax2.plot(years,table_withdrawal_industry[jj,:],label='Beck industry withdrawal',color=(0.5,0,0)) - ax2.plot(years,table_withdrawal_manufacturing[jj,:],label='Beck manufacturing withdrawal',color=(1,0,0)) - ax2.legend() - ax2.set_ylabel('Withdrawal (km^3/year)') - ax2.set_title(country_name) + # # Subpanel 2 + # ax2.plot(years,table_Huang_withdrawal_industry[jj,:],label='Huang industry withdrawal',color=(0,0,0.5)) + # ax2.plot(years,table_Huang_withdrawal_manufacturing[jj,:],label='Huang manufacturing withdrawal',color=(0,0,1)) + # ax2.plot(years,table_withdrawal_industry[jj,:],label='Beck industry withdrawal',color=(0.5,0,0)) + # ax2.plot(years,table_withdrawal_manufacturing[jj,:],label='Beck manufacturing withdrawal',color=(1,0,0)) + # ax2.legend() + # ax2.set_ylabel('Withdrawal (km^3/year)') + # ax2.set_title(country_name) - # Subpanel 3 - ax3.plot(years,table_worldbank_gdp[jj,:],label='GDP',color='g') - ax3.plot(years,table_worldbank_mva[jj,:],label='MVA',color='b') - ax3.legend() - ax3.set_ylabel('Constant 2010 US$') - ax3.set_title(country_name) + # # Subpanel 3 + # ax3.plot(years,table_worldbank_gdp[jj,:],label='GDP',color='g') + # ax3.plot(years,table_worldbank_mva[jj,:],label='MVA',color='b') + # ax3.legend() + # ax3.set_ylabel('Constant 2010 US$') + # ax3.set_title(country_name) - # Save figure - f.set_size_inches(10, 10) - plt.savefig(os.path.join(config['output_folder'],'step3a_industrial_demand','figures',str(jj).zfill(3)+'_'+country_name+'.png'),dpi=150) - plt.close() + # # Save figure + # f.set_size_inches(10, 10) + # plt.savefig(os.path.join(config['output_folder'],'step3a_industrial_demand','figures',str(jj).zfill(3)+'_'+country_name+'.png'),dpi=150) + # plt.close() # Save to csv pd.DataFrame(table_withdrawal_data_source,index=country_codes['name'],columns=['Withdrawal data source']).to_csv(os.path.join(config['output_folder'],'step3a_industrial_demand','tables','withdrawal_data_source.csv')) diff --git a/src/lisfloodutilities/water-demand-historic/step3b_industrial_demand.py b/src/lisfloodutilities/water-demand-historic/step3b_industrial_demand.py index c4b6676..58daab5 100644 --- a/src/lisfloodutilities/water-demand-historic/step3b_industrial_demand.py +++ b/src/lisfloodutilities/water-demand-historic/step3b_industrial_demand.py @@ -18,7 +18,7 @@ import pandas as pd import numpy as np from netCDF4 import Dataset -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt from skimage.transform import resize from tools import * import rasterio @@ -298,27 +298,27 @@ def main(): # Plot figure #-------------------------------------------------------------------------- - # Initialize figure - f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) + # # Initialize figure + # f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) - # Subpanel 1 - ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) - ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') - - # Subpanel 2 - try: - if varname=='ene': varindex = 0 - if varname=='ind': varindex = 1 - timeindex = year-1971 - ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[varindex,:,:,np.arange(timeindex*12,timeindex*12+12)],axis=0)),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) - ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') - except: - pass + # # Subpanel 1 + # ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) + # ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') + + # # Subpanel 2 + # try: + # if varname=='ene': varindex = 0 + # if varname=='ind': varindex = 1 + # timeindex = year-1971 + # ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[varindex,:,:,np.arange(timeindex*12,timeindex*12+12)],axis=0)),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) + # ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') + # except: + # pass - # Save figure - f.set_size_inches(10, 10) - plt.savefig(os.path.join(config['output_folder'],'step3b_industrial_demand','figures',varname+'_'+str(year)+'.png'),dpi=150) - plt.close() + # # Save figure + # f.set_size_inches(10, 10) + # plt.savefig(os.path.join(config['output_folder'],'step3b_industrial_demand','figures',varname+'_'+str(year)+'.png'),dpi=150) + # plt.close() #-------------------------------------------------------------------------- diff --git a/src/lisfloodutilities/water-demand-historic/step4_livestock_demand.py b/src/lisfloodutilities/water-demand-historic/step4_livestock_demand.py index 6e17d40..8099558 100644 --- a/src/lisfloodutilities/water-demand-historic/step4_livestock_demand.py +++ b/src/lisfloodutilities/water-demand-historic/step4_livestock_demand.py @@ -18,7 +18,7 @@ import pandas as pd import numpy as np from netCDF4 import Dataset -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt from skimage.transform import resize from tools import * import rasterio @@ -106,10 +106,10 @@ def main(): for jj in np.arange(country_codes.shape[0]): country_code = country_codes.iloc[jj]['country-code'] for ii in np.arange(len(years)): - sel_ir = (aquastat['m49']==country_code) & (aquastat['Variable Name']=='Irrigation water withdrawal') & (aquastat['Year']==years[ii]) + sel_ir = (aquastat['m49']==country_code) & (aquastat['Variable']=='Irrigation water withdrawal') & (aquastat['Year']==years[ii]) if np.sum(sel_ir)>0: table_aquastat_irrigation_withdrawal[jj,ii] = aquastat['Value'][sel_ir].values # km3/year - sel_ag = (aquastat['m49']==country_code) & (aquastat['Variable Name']=='Agricultural water withdrawal') & (aquastat['Year']==years[ii]) + sel_ag = (aquastat['m49']==country_code) & (aquastat['Variable']=='Agricultural water withdrawal') & (aquastat['Year']==years[ii]) if np.sum(sel_ag)>0: table_aquastat_agriculture_withdrawal[jj,ii] = aquastat['Value'][sel_ag].values # km3/year @@ -291,13 +291,13 @@ def main(): stats.spearmanr(aquastat[sel], huang[sel]) stats.spearmanr(gcam[sel], huang[sel]) - # What do the different distributions look like? - plt.figure(figsize=(8,6)) - plt.hist(aquastat[sel], bins=300, alpha=0.5, range=(0,15), label="aquastat") - plt.hist(gcam[sel], bins=300, alpha=0.5, range=(0,15), label="gcam") - plt.hist(huang[sel], bins=300, alpha=0.5, range=(0,15), label="huang") - plt.legend() - plt.show(block=False) + # # What do the different distributions look like? + # plt.figure(figsize=(8,6)) + # plt.hist(aquastat[sel], bins=300, alpha=0.5, range=(0,15), label="aquastat") + # plt.hist(gcam[sel], bins=300, alpha=0.5, range=(0,15), label="gcam") + # plt.hist(huang[sel], bins=300, alpha=0.5, range=(0,15), label="huang") + # plt.legend() + # plt.show(block=False) # Can we derive a simple correction factor to improve the GCAM+GLW estimates? Not really because mean and median give opposite results... np.mean(aquastat[sel])/np.mean(gcam[sel]) @@ -449,26 +449,26 @@ def main(): # Plot figure #-------------------------------------------------------------------------- - # Initialize figure - f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) + # # Initialize figure + # f, (ax1, ax2) = plt.subplots(2, 1, sharex=True) - # Subpanel 1 - ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) - ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') - - # Subpanel 2 - try: - varindex = 0 - timeindex = year-1971 - ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[varindex,:,:,np.arange(timeindex*12,timeindex*12+12)],axis=0)),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) - ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') - except: - pass + # # Subpanel 1 + # ax1.imshow(np.sqrt(imresize_mean(data_annual_map,(360,720))),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) + # ax1.set_title('Beck et al. (2022) withdrawal (mm/month)') + + # # Subpanel 2 + # try: + # varindex = 0 + # timeindex = year-1971 + # ax2.imshow(np.sqrt(np.sum(Huang_withdrawal[varindex,:,:,np.arange(timeindex*12,timeindex*12+12)],axis=0)),vmin=0,vmax=15,cmap=plt.get_cmap('YlGnBu')) + # ax2.set_title('Huang et al. (2018) withdrawal (mm/month)') + # except: + # pass - # Save figure - f.set_size_inches(10, 10) - plt.savefig(os.path.join(config['output_folder'],'step4_livestock_demand','figures',varname+'_'+str(year)+'.png'),dpi=150) - plt.close() + # # Save figure + # f.set_size_inches(10, 10) + # plt.savefig(os.path.join(config['output_folder'],'step4_livestock_demand','figures',varname+'_'+str(year)+'.png'),dpi=150) + # plt.close() #--------------------------------------------------------------------------