Homec4science
Diffusion solar_potential (aw_solar) Browse

solarPOT/Potential_hourly/merge_and_postprocess.py
29f993e4425aaw_solar

Last Change
Raw File

merge_and_postprocess.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import sys
import os
import time
# # Define file paths
COMM_FP = '/work/hyenergy/output/solar_potential/Sonnendach_Gemeinden_Dachflaechen.csv' # '/Users/alinawalch/Documents/EPFL/data/EMPA_data/Sonnendach_Gemeinden_Dachflaechen.csv'
ROOF_FP = '/work/hyenergy/output/solar_potential/geographic_potential/CH_ROOFS_all_replaced.csv' # '/Volumes/Seagate Backup Plus Drive/DOKUMENTE/EPFL/data/Results/PV_w_uncertainty/CH_ROOFS_all_replaced.csv'
AREA_FP = '/work/hyenergy/output/solar_potential/geographic_potential/available_area/PRED_available_area_unc_filled.csv' # '/Volumes/Seagate Backup Plus Drive/DOKUMENTE/EPFL/data/Results/PV_w_uncertainty/PRED_available_area_unc_filled.csv'
PV_PATH = sys.argv[1]
OUTPATH = sys.argv[2]
## CONSTANTS
GT_IN = 'tilted_irrad_1*.nc'
PV_IN = 'pv_potential_1*.nc'
DF_OUT = '2010_annual_roofs.csv'
HR_OUT = '2010_hourly_comms.nc'
PV_OUT = '2010_hourly_pv_kW.csv'
GT_OUT = '2010_hourly_gt_kW_m2.csv'
FLAT_TILT = 10
MAX_ASPECT = 90
# Print file paths
print('\nCommune mapping file path: %s' %COMM_FP)
print('Available area file path: %s' %AREA_FP)
print('Roof surfaces file path: %s' %ROOF_FP)
print('PV and irradiance input path: %s' %PV_PATH)
print('Output path: %s\n\n' %OUTPATH)
# # Load commune map
tt = time.time()
comm_map = pd.read_csv(COMM_FP, encoding = 'cp1250', usecols = ['DF_UID', 'BFS_NUMMER', 'NAME'] )
# # Load roof data
roofs = pd.read_csv(ROOF_FP, usecols = ['DF_UID', 'NEIGUNG', 'AUSRICHTUNG'])
av_area = pd.read_csv(AREA_FP, usecols = ['DF_UID', 'available_area', 'available_area_std'])
roof_info = roofs.merge(av_area, on = 'DF_UID').merge(comm_map, on = 'DF_UID', how = 'left')
print('Loaded roof info in %.2fs' %(time.time()-tt))
# # Get suitable roofs
tt = time.time()
roof_info['suitable'] = ( ( (roof_info.NEIGUNG <= FLAT_TILT) | (abs(roof_info.AUSRICHTUNG) <= MAX_ASPECT) )
& (roof_info.available_area >= 8) ).astype(int)
# # Load prediction
tilted_irrad = xr.open_mfdataset( os.path.join(PV_PATH, GT_IN), combine = 'by_coords' )
pv_potential = xr.open_mfdataset( os.path.join(PV_PATH, PV_IN), combine = 'by_coords' )
roof_info = roof_info.merge( tilted_irrad.yearly_kWh_m2.to_series().rename('Gt_kWh_m2'), on = 'DF_UID' )
roof_info = roof_info.merge( pv_potential.yearly_PV_kWh.to_series().rename('EPV_kWh'), on = 'DF_UID' )
print('Processed annual PV data in %.2fs\n' %(time.time()-tt))
# ### Save rooftop data
tt = time.time()
roof_save = roof_info.drop(columns = ['NEIGUNG', 'AUSRICHTUNG', 'NAME'])
print(roof_save.head())
roof_save.to_csv( os.path.join(OUTPATH, DF_OUT), index = False )
print('\nSaved annual PV data in %.2fs' %(time.time()-tt))
# # Group suitable potential by communes
tt = time.time()
time_xr = pv_potential.time.isel( DF_UID = 0 ).drop('DF_UID')
roof_sel = roof_info.loc[roof_info.DF_UID.isin(pv_potential.DF_UID.values), ['DF_UID', 'suitable', 'BFS_NUMMER']]
roof_xr = roof_sel.set_index('DF_UID').to_xarray()
hourly_suitable = xr.merge([pv_potential.pv_potential, tilted_irrad.tilted_irradiance]) * roof_xr.suitable
hourly_suitable.coords['BFS_NUMMER'] = roof_xr.BFS_NUMMER.astype('int')
pv_communes_hourly = hourly_suitable.pv_potential.groupby('BFS_NUMMER').sum().astype('float32') / 1000
gt_communes_hourly = ( hourly_suitable.tilted_irradiance.where( hourly_suitable.tilted_irradiance > 0 )
.groupby('BFS_NUMMER').mean().astype('float32') )
communes_output = xr.merge([pv_communes_hourly, gt_communes_hourly, time_xr])
print('\nCreated hourly commune potential in %.2fs\n' %(time.time()-tt))
print(communes_output)
# ### Save intermediate files
tt = time.time()
communes_output.to_netcdf(os.path.join(OUTPATH, HR_OUT))
print('\nComputed hourly commune potential and saved output in %.2fs' %(time.time()-tt))
# ## Create output csv files
tt = time.time()
communes_xr = xr.open_dataset(os.path.join(OUTPATH, HR_OUT))
comm_map_name = comm_map.groupby('BFS_NUMMER').first().NAME
comm_df = communes_xr.to_dataframe().dropna().reset_index()
comm_df['NAME'] = comm_df.BFS_NUMMER.map( comm_map_name )
pv_comm_df = comm_df.pivot(values='pv_potential', index='time', columns='BFS_NUMMER')
gt_comm_df = comm_df.pivot(values='tilted_irradiance', index='time', columns='BFS_NUMMER')
# ### Save data
pv_comm_df.to_csv( os.path.join(OUTPATH, PV_OUT) )
gt_comm_df.to_csv( os.path.join(OUTPATH, GT_OUT) )
print('Saved hourly commune potential and irradiance in %.2fs\n' %(time.time()-tt))
print(pv_comm_df.head())
print(gt_comm_df.head())
c4science · Help