Homec4science
Diffusion solar_potential (master) Browse

solarPOT/EW_scenarios_MMH/poa_irradiance.py
29f993e4425amaster

Last Change
Raw File

poa_irradiance.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import os
import sys
import time
import pvlib.irradiance as pvlib
import pyproj
def round_to_interval(x, round_int):
return np.round(x / round_int) * round_int
def merge_ALB_on_meteo(albedo_df, meteo_df):
# function takes two dataframes of albedo and meteo and matches based on month, x_idx. y_idx
meteo_df['month'] = meteo_df['timestamp'].dt.month
return meteo_df.merge(albedo_df, on = ['x', 'y', 'month']).drop('month', axis = 1)
def map_coords_to_pixel(coords, offset, pixel_res, decimals = 1 ):
# offset should be either scalar or of type np.array()
tmp_coords = coords - offset
tmp_coords_rounded = round_to_interval(tmp_coords, pixel_res)
return (tmp_coords_rounded + offset).round( decimals )
##################### INPUTS ####################
meteo_file = sys.argv[1]
print('\nMeteo data loaded from %s' %meteo_file)
rooftop_file = sys.argv[2]
print('Rooftop data loaded from %s' %rooftop_file)
savefile = sys.argv[3]
print('Data will be saved to %s' %savefile)
'''
#### TMP #####
fp = '/Volumes/Seagate Backup Plus Drive/DOKUMENTE/EPFL/data/'
meteo_file = os.path.join(fp, 'Meteo_data/Meteoswiss_satellite/raw_data/2010_physical_potential_interp.nc')
print('\nMeteo data loaded from %s' %meteo_file)
rooftop_file = os.path.join(fp, 'Results/PV_w_uncertainty/CH_ROOFS_all_replaced.csv')
print('Rooftop data loaded from %s' %rooftop_file)
#savefile = 'hourly_poa_script.nc'
#print('Data will be saved to %s' %savefile)
#### END - TMP #####
'''
if len(sys.argv) > 4:
albedo_file = sys.argv[4]
albedo_in = 'map'
print('Albedo data loaded from %s' %albedo_file)
if len(sys.argv) > 6:
start_index = int(sys.argv[5])
n_rows = int(sys.argv[6])
get_slice = True
print('Loading slice: %d rows, starting from ID %d\n' %(n_rows, start_index))
else:
get_slice = False
start_index = 1
n_rows = None
##################### CONSTANTS ##########################
ROOFTOP_FEATURES = ['XCOORD', 'YCOORD', 'DF_UID', 'FLAECHE', 'panel_tilt', 'panel_direction']
UNIQUE_FEATURE = 'DF_UID'
DIFFUSE_MODEL = 'perez'
albedo_input = {
'constant' : lambda ds : 0.2 ,
'map' : lambda ds : ds['ALB']
}
ZENITH_THRESH = 2 # Threshold for considering 0 direct radiation (in deg. from horizon)
PIXEL_RES = 200
OFFSET = 87.5
################################################
# Create dataframe of rooftops
field_names = pd.read_csv(rooftop_file, nrows = 0).columns
rooftops_all = pd.read_csv(rooftop_file, skiprows = start_index, nrows = n_rows, header = None, names = field_names)
rooftops = rooftops_all.loc[:, ROOFTOP_FEATURES]
if len(rooftops) == 0:
sys.exit('No rooftops found.')
# Open datasets
meteo = xr.open_dataset( meteo_file )
# if only a part of the data is queried, use only the subset of the meteo data in that range, for efficiency reasons
if get_slice:
meteo = meteo.sel( x = slice( rooftops.XCOORD.min()-PIXEL_RES, rooftops.XCOORD.max()+PIXEL_RES), y = slice( rooftops.YCOORD.min()-PIXEL_RES, rooftops.YCOORD.max()+PIXEL_RES) )
if albedo_file is not None:
albedo = xr.open_dataset( albedo_file )
if get_slice:
albedo = albedo.sel( x = slice( rooftops.XCOORD.min()-PIXEL_RES, rooftops.XCOORD.max()+PIXEL_RES), y = slice( rooftops.YCOORD.min()-PIXEL_RES, rooftops.YCOORD.max()+PIXEL_RES) )
# link rooftops with meteo pixels
px_coord = ['x', 'y']
rooftops[px_coord] = map_coords_to_pixel( rooftops[['XCOORD','YCOORD']], OFFSET, PIXEL_RES ).round(1)
tt = time.time()
meteo_df = meteo.to_dataframe().reset_index().dropna()
meteo_df[px_coord] = meteo_df[px_coord].round(1)
print('Created meteo dataframe %.2fs' %(time.time()-tt))
# if an albedo was given, merge a column of ALB to meteo_df
if albedo_file is not None:
albedo_df = albedo.ALB.to_dataframe().reset_index().dropna()
albedo_df[px_coord] = albedo_df[px_coord].round(1)
meteo_df = merge_ALB_on_meteo( albedo_df, meteo_df )
# merge rooftops with meteo data
tt = time.time()
merged_dataset = rooftops.merge(meteo_df, on = px_coord, how = 'left').drop(columns = px_coord)
print('Merged in %.2fs' %(time.time()-tt))
tmp_fp = os.path.join( os.path.split( savefile )[0], 'TMP_merged_dataset.csv')
# merged_dataset.to_csv( tmp_fp )
print('Saved temporary merged dataset in %s' %tmp_fp)
# exit()
################################################
# Compute tilted irradiance for entire dataframe
tt = time.time()
irradiance_components = pvlib.get_total_irradiance(
surface_tilt = merged_dataset['panel_tilt'],
surface_azimuth = merged_dataset['panel_direction'],
solar_zenith = merged_dataset['apparent_zenith'],
solar_azimuth = merged_dataset['azimuth'],
dni = merged_dataset['DNI'],
ghi = merged_dataset['SIS'],
dhi = merged_dataset['DHI'],
dni_extra = merged_dataset['DNI_extra'],
airmass = merged_dataset['rel_airmass'],
albedo = albedo_input[albedo_in](merged_dataset),
model = DIFFUSE_MODEL
)
print('\nComputed irradiance in %.2fs' %(time.time()-tt))
# add UNIQUE_FEATURE and timestamp to the tilted irradiance and drop data that can be calculated (global and diffuse)
tt = time.time()
tilted_irradiance_features = [ UNIQUE_FEATURE, 'timestamp' ]
tilted_irradiance = pd.concat([ merged_dataset.loc[: ,tilted_irradiance_features], irradiance_components], axis = 1)
tilted_irradiance.drop(['poa_global' , 'poa_diffuse'], axis = 1, inplace = True)
print('\nAdded columns for unique identification to dataframe in %.2fs' %(time.time()-tt))
# Clean data
tilted_irradiance.loc[merged_dataset.apparent_zenith >= (90-ZENITH_THRESH), 'poa_direct'] = 0
### Convert pandas dataframe back to xarray to save in a more efficient format
tt = time.time()
data_xr = tilted_irradiance.set_index(['DF_UID', 'timestamp' ]).to_xarray()
data_out = xr.merge([data_xr])
print('Indexed data in %.2fs' %(time.time()-tt))
# Reduce file size
tt = time.time()
for var in data_out.variables:
if var in ['timestamp', 'DF_UID']: continue
data_out[var] = data_out[var].astype('float32')
print('Converted data to float32 in %.2fs' %(time.time()-tt))
print(data_out)
# save output
tt = time.time()
data_out.to_netcdf( savefile )
print('\nSaved file as %s in %.2fs' %(savefile, time.time()-tt))
c4science · Help