Homec4science
Diffusion solar_potential (master) Browse

solarPOT/Potential_hourly/1_compute_physical_potential.py
b9297dda7962master

Last Change
Raw File

1_compute_physical_potential.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import os
import sys
import pvlib.solarposition as solar_pos
from pvlib.irradiance import get_total_irradiance, get_extra_radiation, get_ground_diffuse
from pvlib.atmosphere import get_relative_airmass
import pyproj
# Function definitions
def cosd(angle_in_degrees):
return np.cos(np.radians(angle_in_degrees))
def SISDIR_to_DNI( SISDIR , zenith_in_degrees ):
normal_incidence = SISDIR/cosd(zenith_in_degrees)
return( np.maximum( normal_incidence, 0 ) )
# # INPUT (later argvs)
YEAR = int(sys.argv[1])
OUTPATH = '/work/hyenergy/raw/METEO/MeteoSwiss_satellite/'
fp = '/work/hyenergy/raw/METEO/MeteoSwiss_satellite/raw_data' # '/Volumes/Seagate Backup Plus Drive/DOKUMENTE/EPFL/data/Meteo_data/Meteoswiss_satellite/raw_data/'
# global horizontal irradiance
Gh_file = lambda month,day : os.path.join(fp, 'global/msg.SIS.H_ch02.lonlat_%d%02d%02d000000.nc' %(YEAR,month,day))
# direct beam irradiance
Gb_file = lambda month,day : os.path.join(fp, 'direct/msg.SISDIR.H_ch02.lonlat_%d%02d%02d000000.nc' %(YEAR,month,day))
# Surface albedo
ALB_file = lambda month,day : os.path.join(fp, 'albedo/msg.ALB.H_ch02.lonlat_%d%02d%02d000000.nc' %(YEAR,month,day))
### LOAD AND APPEND DATA
G_all = []
for month in range(1,13):
G_month = []
for day in range(1,32):
try:
Gh = xr.open_dataset(Gh_file(month, day))
Gb = xr.open_dataset(Gb_file(month, day))
ALB = xr.open_dataset(ALB_file(month, day))
G = xr.merge([Gh,Gb, ALB]).rename({'time' : 'timestamp'})
G['time'] = G.timestamp
G.coords['day' ] = day
# Change timestamp to the 15th of each month (including the 30-minute time shift for interp)
G.coords['timestamp'] = [pd.to_datetime('2001-%d-15 %02d:30' %(month, t.dt.hour)) for t in G.timestamp]
G_month.append(G)
print('Successfully appended %d.%d.%d' %(day, month, YEAR))
except:
print('Skipping %d.%d.%d' %(day, month, YEAR))
if len(G_month) > 0:
G_all.append( xr.concat(G_month, dim = 'day') )
print('Successfully added month %d' %(month))
G_all = xr.concat(G_all, dim = 'timestamp')
# ### Interpolate data to represent mean of timestamp
timestamps_interp = (G_all.timestamp - pd.Timedelta(0.5, unit = 'h')).values
G_all_interp = G_all.interp(timestamp = timestamps_interp)
# ### Drop zeros to reduce file size
G_all_nonzero = G_all_interp.where(G_all_interp.SIS > 0, drop = True)
# ### Merge with time
time = G_all.time
time.coords['timestamp'] = timestamps_interp
G_all_interp_w_time = xr.merge([G_all_nonzero, time.sel(timestamp = G_all_nonzero.timestamp)])
# ### Save result
# Reduce file size
for var in G_all_interp_w_time.variables:
print(var)
if var in ['timestamp', 'day', 'lat', 'lon', 'time']: continue
G_all_interp_w_time[var] = G_all_interp_w_time[var].astype('float32')
print(G_all_interp_w_time)
G_all_interp_w_time.to_netcdf(os.path.join(OUTPATH, '%d_SIS_SISDIR_ALB_interp.nc' %YEAR))
print('Saved %s' %(os.path.join(OUTPATH, '%d_SIS_SISDIR_ALB_interp.nc' %YEAR)))
# Free up memory
del G_all_interp_w_time
del G_all_nonzero
del G_all_interp
del G_all
#### COMPUTE SUN POSITION
# Approximate sun position from mean Swiss coordinates
G_all_nonzero = xr.open_dataset(os.path.join(OUTPATH, '%d_SIS_SISDIR_ALB_interp.nc' %YEAR))
times = G_all_nonzero.time.to_dataframe().dropna().reset_index()
# get sun position (one for Switzerland) - can later be replaced by very exact sun positions
LAT = 46.818188
LON = 8.227512
solpos = solar_pos.get_solarposition(times.time.values, LAT, LON)
# derive relative airmass and extraradiation
solpos['rel_airmass'] = get_relative_airmass(solpos['apparent_zenith'])
solpos['DNI_extra' ] = get_extra_radiation(solpos.index)
# Merge data
solpos_merged = solpos.merge(times, left_index = True, right_on = 'time').set_index(['timestamp', 'day'])
meteo = xr.merge([G_all_nonzero,
solpos_merged.to_xarray()[['apparent_zenith', 'azimuth', 'rel_airmass', 'DNI_extra']]])
# Clean data
meteo = meteo.where(meteo.apparent_zenith <= 90, drop = True)
# ## Prepare weather data:
meteo['DHI'] = np.maximum((meteo['SIS'] - meteo['SISDIR']), 0)
meteo['DNI'] = SISDIR_to_DNI(meteo['SISDIR'],meteo['apparent_zenith'])
# ## Geographic transformation to attach x and y values to lon and lat
list_of_locations = G_all_nonzero.isel(timestamp = 0, day = 0).to_dataframe().reset_index()[['lat','lon']]
latlong = pyproj.Proj(proj='latlong')
swiss = pyproj.Proj(init='EPSG:21781')
swiss(600000, 200000, inverse=True)
transformer = lambda lon, lat: pyproj.transform(latlong, swiss, lon, lat)
# transform x and y into lon and lat and append to list_of_locations
list_of_locations['x'], list_of_locations['y'] = transformer(list_of_locations['lon'].values,
list_of_locations['lat'].values)
# select only spatial dimensions and convert to an xarray dataframe
locations_array = list_of_locations.set_index(['lon','lat']).to_xarray()
meteo = xr.merge([meteo, locations_array])
# ## Save output
# Reduce file size
for var in meteo.variables:
print(var)
if var in ['timestamp', 'day', 'lat', 'lon', 'time']: continue
meteo[var] = meteo[var].astype('float32')
print(meteo)
meteo.to_netcdf(os.path.join(OUTPATH, '%d_physical_potential_interp.nc' %YEAR))
print('Saved %s' %(os.path.join(OUTPATH, '%d_physical_potential_interp.nc' %YEAR)))
c4science · Help