Homec4science
Diffusion solar_potential (master) Browse

solarPOT/Potential_hourly/poa_irradiance.py
b9297dda7962master

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 map_coords_to_pixel(coords, offset, pixel_res, col_names = ['lon', 'lat'], decimals = 5 ):
tmp_lon = coords['lon'] - offset[0]
tmp_lat = coords['lat'] - offset[1]
tmp_lon_rounded = round_to_interval(tmp_lon, pixel_res[0])
tmp_lat_rounded = round_to_interval(tmp_lat, pixel_res[1])
return_df = coords.copy()
return_df[col_names[0]] = (tmp_lon_rounded + offset[0]).round( decimals )
return_df[col_names[1]] = (tmp_lat_rounded + offset[1]).round( decimals )
return return_df
##################### 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 #####
'''
albedo_in = sys.argv[4] # 'map' or 'constant'
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)
################################################
# 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 )
# In case ALB is desired but not found in variables: use constant
if not ('ALB' in meteo.variables):
print('Using constant albedo')
albedo_in = 'constant'
dLat = meteo.lat.to_pandas().diff().mean()
dLon = meteo.lon.to_pandas().diff().mean()
# Project rooftops to lonlat
latlong = pyproj.Proj(proj='latlong')
swiss = pyproj.Proj(init='EPSG:21781')
swiss(600000, 200000, inverse=True)
transformer = lambda x, y: pyproj.transform(swiss, latlong, x, y)
# transform x and y into lon and lat and append to list_of_locations
rooftops['lon'], rooftops['lat'] = transformer(rooftops['XCOORD'].values, rooftops['YCOORD'].values)
# 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( lon = slice( rooftops.lon.min()-dLon, rooftops.lon.max()+dLon),
lat = slice( rooftops.lat.min()-dLat, rooftops.lat.max()+dLat) )
# link rooftops with meteo pixels
px_coord = ['lon_px', 'lat_px']
rooftops[px_coord] = map_coords_to_pixel(rooftops[['lon','lat']], [5.75, 45.75], [dLon, dLat], px_coord )[px_coord]
tt = time.time()
meteo_df = meteo.to_dataframe().reset_index().dropna()
meteo_df[px_coord] = map_coords_to_pixel(meteo_df[['lon','lat']], [5.75, 45.75], [dLon, dLat], px_coord )[px_coord]
print('Created meteo dataframe %.2fs' %(time.time()-tt))
# merge rooftops with meteo data
tt = time.time()
merged_dataset = rooftops.merge(meteo_df.drop(columns = ['x','y','lon','lat']), on = px_coord, how = 'left')
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', 'day']
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', 'day']).to_xarray()
data_out = xr.merge([data_xr, meteo.time])
print('Indexed data in %.2fs' %(time.time()-tt))
# Reduce file size
tt = time.time()
for var in data_out.variables:
if var in ['timestamp', 'day', 'time', '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