Page MenuHomec4science
Files F111273368

solar_position.py
No OneTemporary
Actions

File Metadata

Created
Wed, Apr 30, 12:29

solar_position.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import os
import util
import pvlib.solarposition as solar_pos
import pyproj
import itertools
import tqdm
import sys
########## INPUT/argv ##################
reference_file = sys.argv[1]
target_file = sys.argv[2]
###########.CONSTANTS #############
ALWAYS_NONZERO_HOUR = 10
DAILY_DIMENSION_NAME = 'date'
### pre-processing ##
ref = xr.open_dataset(reference_file)
# create list of all locations (for projection to lat and lon)
list_of_locations = ref.isel(hour = ALWAYS_NONZERO_HOUR, date = 0).to_dataframe().dropna().reset_index()
# set up projection from swiss coordinate system CH03
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
list_of_locations['lon'], list_of_locations['lat'] = transformer(list_of_locations['x'].values, list_of_locations['y'].values)
# select only spatial dimensions and convert to an xarray dataframe
locations_array = util.to_xarray(list_of_locations.loc[:,['x','y','lon','lat']], ['x','y'])
### Set up TIME SERIES
# extract only time series information
time_series = ref.isel(x = int(len(ref.x.values)/2), y = int(len(ref.y.values)/2))
time_series = util.to_pandas(time_series, coords = [DAILY_DIMENSION_NAME, 'hour'])
time_series['timestamp'] = util.to_timestamp(time_series[DAILY_DIMENSION_NAME], time_series['hour'])
#### COMPUTE SUN POSITIONS ######
len_x = len(locations_array.x)
len_y = len(locations_array.y)
len_t = len(time_series)
# create dataset of zeroes
zenith = xr.DataArray(data = np.zeros((len_x, len_y, len_t)),
coords = [locations_array.x, locations_array.y, time_series.timestamp],
dims = ['x','y', 'timestamp'],
name = 'apparent_zenith')
azimuth = xr.DataArray(data = np.zeros((len_x, len_y, len_t)),
coords = [locations_array.x, locations_array.y, time_series.timestamp],
dims = ['x','y', 'timestamp'],
name = 'azimuth')
solar_position = xr.merge([zenith, azimuth])
# lambda which applies the function that gets the sun position for a given location
solar_position_at_pixel = lambda ix,iy : solar_pos.get_solarposition(
time_series.timestamp,
latitude = float( locations_array.lat[ ix,iy ]),
longitude = float( locations_array.lon[ ix,iy ]))
# function which sets the correct positions in the solar_position dataset to the actual sun positions
def assign_values(ix, iy):
timeseries = solar_position_at_pixel(ix,iy)
solar_position.apparent_zenith[ix, iy, :] = timeseries.apparent_zenith
solar_position.azimuth[ix, iy, :] = timeseries.azimuth
# for-loop which loops over all locations
tt = util.Timer()
for ix,iy in tqdm.tqdm(itertools.product(range(len_x),range(len_y))):
if locations_array.lon[ix, iy] is not None:
assign_values(ix,iy)
tt.stop()
# save the result:
solar_position.to_netcdf(target_file)

Event Timeline

c4science · Help