Homec4science
Diffusion solar_potential (master) Browse

solarPOT/Uncertainty/uncertainty_Gt2.py
cd293a1ab61cmaster

Last Change
Raw File

uncertainty_Gt2.py
View Options

import xarray as xr
import pandas as pd
import numpy as np
import os
import sys
import time
start = int(sys.argv[1])
end = min(int(sys.argv[2]), 9639231)
batch = int(sys.argv[3])
outfp = sys.argv[4]
print('Processing rooftops %d to %d' %(start, end))
t0 = time.time()
# Set constants
CHUNKSIZE = 10000 # only use 1eX
SELECT = range(start, end) # max value: 9639231
# SET FILEPATHS
data_path = '/scratch/walch/workspace_uncertainty/'
phys_pot_fp = os.path.join( data_path, 'tmp/tilted_irrad_interp_proj_*.csv' )
sigma_Gh_fp = os.path.join( data_path, 'SIS_dhm25_200_CH_ELM50_1000_CPU_total_unc_interp_roofs.nc' )
sigma_GB_fp = os.path.join( data_path, 'SISDIR_dhm25_200_CH_ELM50_1000_CPU_total_unc_interp_roofs.nc' )
cov_G_fp = os.path.join( data_path, 'covrr_Gh_GB_nearest_interp_roofs.nc' )
G_poa_fp = '/scratch/walch/workspace_tilted_irrad/tilted_irradiance_interp_clean.nc'
Ssh_fp = os.path.join( data_path, 'CH_partly_shading_replaced.nc' )
SVF_fp = '/work/hyenergy/output/solar_potential/geographic_potential/tilted_irradiance/SVF_CH_replaced.csv'
SVF_bias_fp = os.path.join( data_path, 'SVF_RV_mean_std.csv' )
cov_DSM_fp = os.path.join( data_path, 'covariance_DSM.csv' )
PV_pot_fp = '/scratch/walch/workspace_tilted_irrad/PV_all_interp.nc'
# # Compute uncertainty for tilted radiation
# ## Step 1: load data
# 1. load physical potential ($G_h, G_B, G_D$), which was mapped to roof surfaces from physical_potential...nc
# 2. load (total) variance on physical potential ($\sigma^2_{Gh}, \sigma^2_{GB}$)
# 3. load covariance between $G_h$ and $G_B$ ($Cov(G_{h}, G_B)$)
# 4. load tilted irradiance ($G_{Bt}, G_{Dt}, G_{Rt}$), so that $R_B, R_D, R_R$ can be extracted
# 5. load partly shading factor $S_{sh}$, its bias and variance $\sigma^2_{Sh}$
# 6. load sky view factor $SVF$, its bias and variance $\sigma^2_{SVF}$
# 7. load the covariance of $S_{sh}$ and $SVF$ ($Cov(S_{sh}, SVF$)
# 8. load the total tilted irradiance $G_t$ to get the uncertainty in percentage
# ### 1.1 Mapped physical potential
phys_pot_mapped = xr.open_mfdataset( phys_pot_fp, chunks = {'DF_UID' : CHUNKSIZE} ).sortby('DF_UID').isel(DF_UID = SELECT)
phys_pot_mapped = phys_pot_mapped[['SIS', 'SISDIR', 'DHI']]
G = phys_pot_mapped.rename({'SIS' : 'Gh', 'SISDIR' : 'GB', 'DHI' : 'GD'})
# ### 1.2 Uncertainties on physical potential
sigma_Gh = xr.open_mfdataset(sigma_Gh_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
sigma_GB = xr.open_mfdataset(sigma_GB_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
var_G = xr.merge([sigma_Gh.rename({'SIS' : 'sigma_Gh'}) , sigma_GB.rename({'SISDIR' : 'sigma_GB'}),
sigma_Gh.rename({'SIS' : 'var_Gh'})**2, sigma_GB.rename({'SISDIR' : 'var_GB'})**2])
# ### 1.3 Covariance of Gh and GB
cov_G = xr.open_mfdataset(cov_G_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
cov_G = cov_G[['corr_Gh_GB']]
# ### 1.4 Tilted irradiance
G_poa = xr.open_mfdataset(G_poa_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
# ### 1.5 Partly shaded factor with bias and uncertainty
partly_shading = xr.open_mfdataset(Ssh_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
Ssh_bias_unc = pd.read_csv( os.path.join( data_path, 'Ssh_RV_mean_std.csv' ))
# get timestamp of useful hours and convert to xarray
Ssh_bias_unc = Ssh_bias_unc[Ssh_bias_unc['std'] > 0]
Ssh_bias_unc['timestamp'] = [pd.Timestamp(year = 2001, month = int(np.floor(idx/24)+1), day = 15, hour = int(idx%24))
for idx in Ssh_bias_unc.hour_idx]
Ssh_bias_xr = Ssh_bias_unc.set_index('timestamp').to_xarray()
# debias S_sh by subtracting the mean (mean was created as 2m - 50cm, hence subtraction and not addition)
Ssh = xr.merge([(partly_shading['partly_shading'] - Ssh_bias_xr['mean'] ).rename('S_Sh'),
Ssh_bias_xr['std'].rename('var_Sh')**2])
# ### 1.6 SVF with bias and uncertainty
SVF_data = pd.read_csv(SVF_fp)
SVF_mean_bias = pd.read_csv( SVF_bias_fp, index_col = 0)
# get debiased expected value and variance
# get debiased expected value and variance
svf_mean = xr.DataArray(data = SVF_data['SVF_unbiased'] - float(SVF_mean_bias['mean']),
dims = 'DF_UID', coords = {'DF_UID' : SVF_data.DF_UID}, name = 'SVF')
svf_var = xr.DataArray(data = float(SVF_mean_bias['std']**2), name = 'var_SVF')
SVF = xr.merge([svf_mean, svf_var])
# ### 1.7 Covariance between SVF and Ssh
cov_DSM = pd.read_csv( cov_DSM_fp, index_col = 0)
cov_DSM.drop(['Csh', 'SVF'], inplace = True)
cov_DSM['timestamp'] = [pd.Timestamp(year = 2001, month = int(np.floor(idx/24)+1), day = 15, hour = int(idx%24))
for idx in cov_DSM.index.values.astype(int)]
cov_Sh_SVF = cov_DSM.set_index('timestamp')[['cov_SVF']].to_xarray().rename({'cov_SVF' : 'cov_Sh_SVF'})
# ### 1.8 Load tilted radiation data
PV_pot = xr.open_mfdataset(PV_pot_fp, chunks = {'DF_UID' : CHUNKSIZE}).sortby('DF_UID').isel(DF_UID = SELECT)
Gt = PV_pot[['tilted_irradiance']].rename({'tilted_irradiance' : 'Gt'})
print('Loaded all datasets')
# # Step 2: Data processing
# 1. Compute variance and covariances of $G_D$ from the data of $G_h$ and $G_B$
# 2. Compute tilted factors $R_{Bt}, R_{Dt}, R_{Rt}$ from the poa radiation and physical potential
# 3. Compute the variances of the direct, diffuse and direct components ($\sigma^2_{Sh-Bt}, \sigma^2_{SVF-Dt}, \sigma^2_{Rt}$)
# 4. Compute the covariances between the direct components ($Cov(G_{Sh-Bt}, G_{Rt}), Cov(G_{Sh-Bt}, G_{SVF-Dt}), Cov(G_{SVF-Dt}, G_{Rt})$)
# 5. Compute the variance of the total tilted radiation ($\sigma^2_{Gt}$)
# 6. Compare with the total of the tilted radiation $G_t$
# ### 2.1 Computation of the variance of G_D and the covariances
# compute covariance from correlation coefficient and standard deviations to be consistent
cov_G['cov_Gh_GB'] = cov_G.corr_Gh_GB * sigma_Gh.SIS * sigma_GB.SISDIR
var_G['var_GD'] = xr.ufuncs.maximum( var_G.var_Gh + var_G.var_GB - 2 * cov_G.cov_Gh_GB, 0 )
var_G['sigma_GD'] = xr.ufuncs.maximum( xr.ufuncs.sqrt(var_G.var_GD), 0 )
cov_G['cov_GB_GD'] = 0.5 * (var_G.var_Gh - var_G.var_GB - var_G.var_GD)
cov_G['cov_Gh_GD'] = 0.5 * (var_G.var_Gh - var_G.var_GB + var_G.var_GD)
cov_G['corr_GB_GD'] = cov_G.cov_GB_GD / (var_G.sigma_GD * var_G.sigma_GB)
cov_G['corr_Gh_GD'] = cov_G.cov_Gh_GD / (var_G.sigma_GD * var_G.sigma_Gh)
# ### 2.2 Computation of tilted factors R
R_B = G_poa.poa_direct / G.GB
R_D = G_poa.poa_sky_diffuse / G.GD
R_R = G_poa.poa_ground_diffuse / G.Gh
R = xr.merge([ R_B.rename('RBt'), R_D.rename('RDt'), R_R.rename('RRt')])
# ### 2.3 Computation of variances of tilted components
# $\sigma^2_{Sh-Bt} = R_{Bt}^2 (\sigma^2_{Sh} G_B^2 + \sigma^2_{GB} S_{sh}^2 + \sigma^2_{Sh} \sigma^2_{GB})$
#
# $\sigma^2_{SVF-Dt} = R_{Dt}^2 (\sigma^2_{SVF} G_D^2 + \sigma^2_{GD} SVF^2 + \sigma^2_{SVF} \sigma^2_{GD})$
#
# $\sigma^2_{Rt} = R_{Rt}^2 \sigma^2_{G_h}$
var_Sh_Bt = R.RBt**2 * (Ssh.var_Sh * G.GB**2 + var_G.var_GB * Ssh.S_Sh**2 + Ssh.var_Sh * var_G.var_GB)
var_SVF_Dt = R.RBt**2 * (SVF.var_SVF * G.GD**2 + var_G.var_GD * SVF.SVF**2 + SVF.var_SVF * var_G.var_GD)
var_Rt = R.RRt**2 * var_G.var_Gh
var_Xt = xr.merge([var_Sh_Bt.rename('var_Sh_Bt'), var_SVF_Dt.rename('var_SVF_Dt'), var_Rt.rename('var_Rt')])
# ### 2.4 Computation of covariances of tilted components
# $Cov(G_{SBt}, G_{Rt}) = R_{Bt} R_{Rt} S_{sh} Cov(G_B, G_h)$
#
# $Cov(G_{SDt}, G_{Rt}) = R_{Dt} R_{Rt} SVF Cov(G_D, G_h)$
#
# $Cov(G_{SBt}, G_{SDt}) = R_{Bt} R_{Rt} S_{sh} SVF Cov(G_B, G_D)$ in the case of **independence** between $S_{sh}$ and $SVF$
#
# $Cov(G_{SBt}, G_{SFt}) = R_{Bt} R_{Rt} (Cov(G_B, G_D) Cov(S_{sh}, SVF) + G_B G_D Cov(S_{sh}, SVF) + S_{sh} SVF Cov(G_B, G_D)$ in the case of **dependence**
cov_GSBT_GRt = Ssh.S_Sh * R.RBt * R.RRt * cov_G.cov_Gh_GB
cov_GSDT_GRt = SVF.SVF * R.RDt * R.RRt * cov_G.cov_Gh_GD
# assume INDEPENDENCE between S_Sh and SVF
cov_GSBT_GSDT_indep = Ssh.S_Sh * SVF.SVF * R.RDt * R.RRt * cov_G.cov_GB_GD
# assume DEPENDENCE between S_Sh and SVF
cov_GSBT_GSDT_dep = R.RDt * R.RRt * (cov_G.cov_GB_GD * cov_Sh_SVF.cov_Sh_SVF
+ G.GB * G.GD * cov_Sh_SVF.cov_Sh_SVF
+ Ssh.S_Sh * SVF.SVF * cov_G.cov_GB_GD)
cov_Xt = xr.merge([cov_GSBT_GRt.rename('cov_GSBT_GRt'), cov_GSDT_GRt.rename('cov_GSDT_GRt'),
cov_GSBT_GSDT_dep.rename('cov_GSBT_GSDT_dep'), cov_GSBT_GSDT_indep.rename('cov_GSBT_GSDT_indep')])
# ### 2.5 Computation of variance of tilted radiation
var_Gt_indep = xr.ufuncs.maximum( var_Xt.var_Sh_Bt + var_Xt.var_SVF_Dt + var_Xt.var_Rt
+ 2 * (cov_Xt.cov_GSBT_GRt + cov_Xt.cov_GSDT_GRt + cov_Xt.cov_GSBT_GSDT_indep ), 0)
var_Gt_dep = xr.ufuncs.maximum( var_Xt.var_Sh_Bt + var_Xt.var_SVF_Dt + var_Xt.var_Rt
+ 2 * (cov_Xt.cov_GSBT_GRt + cov_Xt.cov_GSDT_GRt + cov_Xt.cov_GSBT_GSDT_dep ), 0)
var_Gt = xr.merge([var_Gt_indep.rename('var_Gt_indep'), var_Gt_dep.rename('var_Gt_dep')])
var_Gt['sigma_Gt_dep'] = np.maximum(xr.ufuncs.sqrt(var_Gt.var_Gt_dep), 0)
# get rid of na values (all na) and interpolate the remaining NaNs
var_Gt = var_Gt.interpolate_na(dim='timestamp')
var_Gt = var_Gt.where(var_Gt > 0).fillna(0)
print(var_Gt)
tt = time.time()
fp = os.path.join( outfp, 'var_Gt_%d.nc' %(batch))
var_Gt.to_netcdf(fp)
print('Saved %s in %.3fs' %(fp, time.time()-tt))
"""
# ## Step 3: Analyse results and compare with tilted radiation
Gt['sigma_Gt_dep'] = np.maximum(xr.ufuncs.sqrt(var_Gt.var_Gt_dep), 0)
Gt['sigma_Gt_perc'] = Gt.sigma_Gt_dep / Gt.Gt
Gt['Gt_min'] = xr.ufuncs.maximum(Gt.Gt - Gt.sigma_Gt_dep, 0)
Gt['Gt_max'] = Gt.Gt + Gt.sigma_Gt_dep
# get rid of na values (all na) and interpolate the remaining NaNs
Gt = Gt.dropna(dim = 'timestamp', how = 'all')
Gt = Gt.interpolate_na(dim='timestamp')
Gt = Gt.where(Gt > 0).fillna(0)
print('Computed uncertainty; saving files in %s' %outfp)
tt = time.time()
fp = os.path.join( outfp, 'tilted_irrad_unc_%d.nc' %(batch))
Gt.to_netcdf(fp)
print('Saved %s in %.3fs' %(fp, time.time()-tt))
big_file = xr.merge([Gt, var_G, cov_G, var_Gt, var_Xt, cov_Xt])
tt = time.time()
fp = os.path.join( outfp, 'all_vars_%d.nc' %(batch))
big_file.to_netcdf(fp)
print('Saved %s in %.3fs' %(fp, time.time()-tt))
tt = time.time()
fp = os.path.join( outfp, 'var_G_%d.nc' %(batch))
var_G.to_netcdf(fp)
print('Saved %s in %.3fs' %(fp, time.time()-tt))
tt = time.time()
fp = os.path.join( outfp, 'cov_G_%d.nc' %(batch))
# cov_G.to_netcdf(fp)
# print('Saved %s in %.3fs' %(fp, time.time()-tt))
tt = time.time()
fp = os.path.join( outfp, 'var_Xt_%d.nc' %(batch))
# var_Xt.to_netcdf(fp)
# print('Saved %s in %.3fs' %(fp, time.time()-tt))
tt = time.time()
fp = os.path.join( outfp, 'cov_Xt_%d.nc' %(batch))
# cov_Xt.to_netcdf(fp)
# print('Saved %s in %.3fs' %(fp, time.time()-tt))
"""
print('Total execution time: %.3fs' %(time.time()-t0))
c4science · Help