From 691d353cde175663cc248509fab9bbeac5b66bc4 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Fri, 19 Jul 2019 16:05:59 -0700 Subject: [PATCH 01/29] first pass at horizon code in pvlib --- pvlib/horizon.py | 387 ++++++++++++++++++++++++++++++++++++++++++++ pvlib/irradiance.py | 32 +++- pvlib/location.py | 6 +- pvlib/modelchain.py | 25 ++- 4 files changed, 444 insertions(+), 6 deletions(-) create mode 100644 pvlib/horizon.py diff --git a/pvlib/horizon.py b/pvlib/horizon.py new file mode 100644 index 0000000000..5767c6c5e9 --- /dev/null +++ b/pvlib/horizon.py @@ -0,0 +1,387 @@ +import random +import pytz +import time +import itertools + +import numpy as np +import scipy as sp +from scipy.interpolate import RegularGridInterpolator +import pandas as pd + +import matplotlib.pyplot as plt + +import googlemaps + + +def latitude_to_geocentric(phi): + a = 6378.137 + b = 6356.752 + return np.arctan(b**2/a**2*np.tan(phi)) + +def latitude_to_geodetic(phi): + a = 6378.137 + b = 6356.752 + return np.arctan(a**2/b**2*np.tan(phi)) + +def xyz_from_lle(point): + lat = point[0] + lon = point[1] + elev = point[2] + + a = 6378137.0 + b = 6356752.0 + + # convert to radians + phi = lat*np.pi/180.0 + theta = lon*np.pi/180.0 + + # compute radius of earth at each point + r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 + r = r / (a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2) + r = np.sqrt(r) + + h = r + elev + alpha = latitude_to_geocentric(phi) + beta = theta + x = h * np.cos(alpha) * np.cos(beta) + y = h * np.cos(alpha) * np.sin(beta) + z = h * np.sin(alpha) + v = np.array((x, y, z)) + return v + +def lle_from_xyz(point): + a = 6378137.0 + b = 6356752.0 + + x = point[0] + y = point[1] + z = point[2] + + # get corresponding point on earth's surface + t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) + point_s = t * point + x_s = point_s[0] + y_s = point_s[1] + z_s = point_s[2] + + elev = np.linalg.norm(point-point_s) + r = np.linalg.norm(point_s) + + alpha = np.arcsin(z_s / r) + phi = latitude_to_geodetic(alpha) + lat = phi*180.0/np.pi + + lon = np.arctan2(y, x)*180/np.pi + return (lat, lon, elev) + +def pol2cart(rho, phi): + x = rho * np.cos(phi) + y = rho * np.sin(phi) + return(x, y) + +def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): + ''' + input lat long + output grid of lat,long tuples + ''' + grid = np.ndarray((grid_size + 1, grid_size + 1, 2)) + + # fill out grid + for i in range(grid_size + 1): + for j in range(grid_size + 1): + grid[i,j,0] = lat + (i - grid_size / 2) * grid_step + grid[i,j,1] = lon + (j - grid_size / 2) * grid_step + + return grid + +def get_grid_elevations(in_grid, api_key): + ''' + input grid of lat,lon tuples + output grid of LLE tuples + ''' + in_shape = in_grid.shape + lats = in_grid.T[0].flatten() + longs = in_grid.T[1].flatten() + locations = zip(lats, longs) + gmaps = googlemaps.Client(key=api_key) + + out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) + + # Get elevation data from gmaps + elevations = [] + responses = [] + + while len(locations) > 512: + locations_to_request = locations[:512] + locations = locations[512:] + responses += gmaps.elevation(locations=locations_to_request) + responses += gmaps.elevation(locations=locations) + for entry in responses: + elevations.append(entry["elevation"]) + + for i in range(in_shape[0]): + for j in range(in_shape[1]): + lat = in_grid[i,j,0] + lon = in_grid[i,j,1] + elevation = elevations[i + j * in_shape[1]] + + out_grid[i,j,0] = lat + out_grid[i,j,1] = lon + out_grid[i,j,2] = elevation + return out_grid + +def dip_calc(pt1, pt2): + ''' + input: two LLE tuples + output: distance, dip angle, azimuth + ''' + a = 6378137.0 + b = 6356752.0 + + lat1 = pt1[0] + lon1 = pt1[1] + elev1 = pt1[2] + lat2 = pt2[0] + lon2 = pt2[1] + elev2 = pt2[2] + + # convert to radians + phi1 = lat1*np.pi/180.0 + theta1 = lon1*np.pi/180.0 + phi2 = lat2*np.pi/180.0 + theta2 = lon2*np.pi/180.0 + + v1 = xyz_from_lle((lat1, lon1, elev1)) + v2 = xyz_from_lle((lat2, lon2, elev2)) + + x1 = v1[0] + y1 = v1[1] + z1 = v1[2] + x2 = v2[0] + y2 = v2[1] + z2 = v2[2] + + delta = np.subtract(v1,v2) + distance = np.linalg.norm(delta) + + normal = np.array((2*x1/a**2, 2*y1/a**2, 2*z1/b**2)) + beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta)/np.linalg.norm(normal)) + dip_angle = beta - np.pi/2 + dip_angle_deg = dip_angle*180.0/np.pi + + # might wanna double check this formula (haversine?) + azimuth = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) + azimuth_deg = azimuth*180.0/np.pi + + + return (azimuth_deg, dip_angle_deg) + +def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): + ''' + input grid of lat,lon,elevation tuples + output list of azimuth, distance, dip angle + use grid points in first pass and then try randomly sampling triangle method + ''' + + grid_shape = grid.shape + grid_center_i = (grid_shape[0] - 1) / 2 + grid_center_j = (grid_shape[1] - 1) / 2 + site_lat = grid[grid_center_i, grid_center_j, 0] + site_lon = grid[grid_center_i, grid_center_j, 1] + site_elev = grid[grid_center_i, grid_center_j, 2] + site = (site_lat, site_lon, site_elev) + + horizon_points = [] + start = time.time() + if sampling_method == "grid": + samples = sample_using_grid(grid) + elif sampling_method == "triangles": + samples = sample_using_triangles(grid, sampling_param) + elif sampling_method == "interpolator": + samples = sample_using_interpolator(grid, sampling_param) + post_sampling = time.time() + print("Sampling took " + str(post_sampling-start) + " sec") + + dip_calc_lambda = lambda pt: dip_calc(site, pt) + horizon_points = np.array(list(map(dip_calc_lambda, samples))) + + post_calcs = time.time() + print("Dip calcs on samples took " + str(post_calcs-post_sampling) + " sec") + return horizon_points + +def sample_using_grid(grid): + # use every grid point as a sample + # returns a list of LLE tuples + grid_shape = grid.shape + grid_center_i = (grid_shape[0] - 1) / 2 + grid_center_j = (grid_shape[1] - 1) / 2 + + all_samples = [] + for i in range(grid_shape[0]): + for j in range(grid_shape[1]): + # make sure the site is excluded + if i != grid_center_i or j != grid_center_j: + lat = grid[i, j, 0] + lon = grid[i, j, 1] + elev = grid[i, j, 2] + all_samples.append((lat, lon, elev)) + return all_samples + +def sample_using_triangles(grid, samples_per_triangle=10): + # uniformly sample all triangles between neighboring grid points + # returns a list of LLE tuples + + all_samples = [] + for i in range(grid.shape[0]): + for j in range(grid.shape[1]): + center = (grid[i,j,0], grid[i,j,1], grid[i,j,2]) + if i != 0 and j != 0: + left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) + top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) + all_samples += uniformly_sample_triangle(center, top, left , samples_per_triangle) + + if i != 0 and j != grid.shape[1] - 1: + right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) + top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) + all_samples += uniformly_sample_triangle(center, top, right , samples_per_triangle) + + if i != grid.shape[0] - 1 and j != 0: + left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) + bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) + all_samples += uniformly_sample_triangle(center, bottom, left , samples_per_triangle) + + if i != grid.shape[0] - 1 and j != grid.shape[1] - 1: + right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) + bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) + all_samples += uniformly_sample_triangle(center, bottom, right , samples_per_triangle) + return all_samples + +def sample_using_interpolator(grid, num_samples): + x = grid.T[0][0] + y = grid.T[1].T[0] + + x_base = x[0] + x_range = x[-1] - x[0] + y_base = y[0] + y_range = y[-1] - y[0] + + grid_shape = grid.shape + grid_center_i = (grid_shape[0] - 1) / 2 + grid_center_j = (grid_shape[1] - 1) / 2 + site_lat = grid[grid_center_i, grid_center_j, 0] + site_lon = grid[grid_center_i, grid_center_j, 1] + + elevs = grid.T[2].T + interpolator = RegularGridInterpolator((x,y), elevs) + + + r = np.linspace(0, x_range/2, num_samples[0]) + theta = np.linspace(0, 2 * np.pi, num_samples[1]) + polar_pts = np.array(list(itertools.product(r, theta))) + + pts = np.array([pol2cart(e[0], e[1]) for e in polar_pts]) + pts += np.array((site_lat, site_lon)) + + interpolated_elevs = interpolator(pts).reshape(num_samples[0]*num_samples[1], 1) + samples = np.concatenate((pts, interpolated_elevs), axis=1) + return samples + + +def uniformly_sample_triangle(p1, p2, p3, num_samples): + # returns a list of LLE tuples + c1 = xyz_from_lle(p1) + c2 = xyz_from_lle(p2) + c3 = xyz_from_lle(p3) + + points = [] + for i in range(num_samples): + r1 = np.random.rand() + r2 = np.random.rand() + sqrt_r1 = np.sqrt(r1) + + # use uniform sampling from http://www.sherrytowers.com/randomly_sample_points_in_triangle.pdf + random_pt = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 + points.append(lle_from_xyz(random_pt)) + return points + +def round_to_nearest(x, base): + return base * round(float(x) / base) + +def filter_points(horizon_points, bucket_size=1): + wedges = {} + for pair in horizon_points: + azimuth = pair[0] + dip = pair[1] + azimuth_wedge = round_to_nearest(azimuth, bucket_size) + + if azimuth_wedge in wedges: + wedges[azimuth_wedge] = max(dip, wedges[azimuth_wedge]) + else: + wedges[azimuth_wedge] = dip + + filtered_points = [] + for key in wedges.keys(): + filtered_points.append((key, wedges[key])) + + sorted_points = sorted(filtered_points, key=lambda tup: tup[0]) + return sorted_points + +def visualize(horizon_profile, pvsyst_scale=False): + + azimuths = [] + dips = [] + for pair in horizon_profile: + azimuth = pair[0] + azimuths.append(azimuth) + dips.append(pair[1]) + plt.figure(figsize=(10,6)) + if pvsyst_scale: + plt.ylim(0, 90) + plt.plot(azimuths, dips, "-") + plt.show + + +def polar_plot(horizon_profile): + + azimuths = [] + dips = [] + for pair in horizon_profile: + azimuth = pair[0] + azimuths.append(np.radians(azimuth)) + dips.append(pair[1] + 5) + plt.figure(figsize=(10,6)) + sp = plt.subplot(1, 1, 1, projection='polar') + sp.set_theta_zero_location('N') + sp.set_theta_direction(-1) + plt.plot(azimuths, dips, "o") + plt.show + +def horizon_table(horizon_points): + for pair in horizon_points: + azimuth = pair[0] + print(str(azimuth) + ": " + str(pair[1])) + +def invert_for_pvsyst(horizon_points, hemisphere="north"): + # look at that northern hemisphere bias right there + # not even sorry. + assert hemisphere == "north" or hemisphere == "south" + + inverted_points = [] + for pair in horizon_points: + azimuth = pair[0] + if hemisphere == "north": + azimuth -= 180 + if azimuth < -180: + azimuth += 360 + elif hemisphere == "south": + azimuth = -azimuth + inverted_points.append((azimuth, pair[1])) + sorted_points = sorted(inverted_points, key=lambda tup: tup[0]) + return sorted_points + +def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): + grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) + elev_grid = get_grid_elevations(grid, GMAPS_API_KEY) + horizon_points = calculate_horizon_points(elev_grid, sampling_method="interpolator", sampling_param=(1000,1000)) + filtered_points = filter_points(horizon_points, bucket_size=1) + return filtered_points diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 000cf3216c..553b5fa3cf 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -365,7 +365,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth, poa_sky_diffuse = get_sky_diffuse( surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=dni_extra, airmass=airmass, model=model, - model_perez=model_perez) + model_perez=model_perez, **kwargs) poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo, surface_type) @@ -383,7 +383,8 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=None, airmass=None, model='isotropic', - model_perez='allsitescomposite1990'): + model_perez='allsitescomposite1990', + **kwargs): r""" Determine in-plane sky diffuse irradiance component using the specified sky diffuse irradiance model. @@ -444,6 +445,9 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, sky = perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith, solar_azimuth, airmass, model=model_perez) + elif model == 'horizon_adjusted': + assert("horizon_profile" in kwargs) + sky = horizon_adjusted(dhi, horizon_profile, surface_tilt, surface_azimuth) else: raise ValueError('invalid model selection {}'.format(model)) @@ -1205,6 +1209,30 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, return sky_diffuse +def calculate_dtf(dip_angles, surface_tilt, surface_azimuth): + tilt_rad = np.radians(surface_tilt) + plane_az_rad = np.radians(surface_azimuth) + a = np.sin(tilt_rad) * np.cos(plane_az_rad) + b = np.sin(tilt_rad) * np.sin(plane_az_rad) + c = np.cos(tilt_rad) + dtf = 0 + for pair in dip_angles: + az = np.radians(pair[0]) + dip = np.radians(pair[1]) + first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) + second_term = .5 * c * np.cos(dip)**2 + dtf += (first_term + second_term) / 180 + return dtf + +def horizon_adjusted(dhi, horizon_profile, surface_tilt, surface_azimuth): + + dtf = calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) + sky_diffuse = dhi * dtf + + return sky_diffuse + + + def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0): """ Calculate the clearsky index. diff --git a/pvlib/location.py b/pvlib/location.py index 68b2062907..ca0753fb83 100644 --- a/pvlib/location.py +++ b/pvlib/location.py @@ -58,7 +58,7 @@ class Location(object): """ def __init__(self, latitude, longitude, tz='UTC', altitude=0, - name=None, **kwargs): + name=None, horizon_profile=None, **kwargs): self.latitude = latitude self.longitude = longitude @@ -76,7 +76,7 @@ def __init__(self, latitude, longitude, tz='UTC', altitude=0, raise TypeError('Invalid tz specification') self.altitude = altitude - + self.horizon_profile = horizon_profile self.name = name def __repr__(self): @@ -317,3 +317,5 @@ def get_sun_rise_set_transit(self, times, method='pyephem', **kwargs): 'one of pyephem, spa, geometric' .format(method)) return result + + \ No newline at end of file diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index 1455a57ad3..a0a6994efc 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -715,10 +715,30 @@ def no_extra_losses(self): self.losses = 1 return self + def adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith): + # assumes horizon_profile is in 1 deg. azimuth increments + dip_angles = {} + for pair in horizon_profile: + dip_angles[pair[0]] = pair[1] + adjusted_irradiance = poa_direct + for i in range(len(poa_direct)): + rounded_solar_azimuth = round(solar_azimuth[i]) + horizon_dip = dip_angles[rounded_solar_azimuth] + if (solar_zenith > (90 - horizon_dip)): + adjusted_irradiance[i] = 0 + + return adjusted_irradiance + + def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) + poa_direct = self.total_irrad['poa_direct'] + if self.location.horizon_profile != None: + solar_zenith = self.solar_position['apparent_zenith'] + solar_azimuth = self.solar_position['azimuth'] + poa_direct = adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith) self.effective_irradiance = self.spectral_modifier * ( - self.total_irrad['poa_direct']*self.aoi_modifier + + poa_direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) return self @@ -889,7 +909,8 @@ def prepare_inputs(self, times=None, weather=None): self.weather['ghi'], self.weather['dhi'], airmass=self.airmass['airmass_relative'], - model=self.transposition_model) + model=self.transposition_model, + horizon_profile=self.location.horizon_profile) if self.weather.get('wind_speed') is None: self.weather['wind_speed'] = 0 From 99426ba167ee991766ee410767ed8c8f0b7784aa Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 23 Jul 2019 15:18:17 -0700 Subject: [PATCH 02/29] added most of horizon shading code --- pvlib/__init__.py | 1 + pvlib/horizon.py | 8 ++++++ pvlib/irradiance.py | 60 +++++++++++++++++++++++++++++++++++++++------ pvlib/modelchain.py | 22 ++++++----------- setup.py | 4 +-- 5 files changed, 70 insertions(+), 25 deletions(-) diff --git a/pvlib/__init__.py b/pvlib/__init__.py index 238d481615..2f96da15b2 100644 --- a/pvlib/__init__.py +++ b/pvlib/__init__.py @@ -12,6 +12,7 @@ from pvlib import spa from pvlib import modelchain from pvlib import singlediode +from pvlib import horizon # for backwards compatibility for pvlib.tmy module from pvlib import tmy diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 5767c6c5e9..4fe91d3eed 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -385,3 +385,11 @@ def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): horizon_points = calculate_horizon_points(elev_grid, sampling_method="interpolator", sampling_param=(1000,1000)) filtered_points = filter_points(horizon_points, bucket_size=1) return filtered_points + + +def fake_horizon_profile(max_dip): + fake_profile = [] + for i in range(-180, 181): + fake_profile.append((i, random.random()*max_dip)) + + return fake_profile \ No newline at end of file diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 553b5fa3cf..9bb680fdef 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -447,7 +447,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model=model_perez) elif model == 'horizon_adjusted': assert("horizon_profile" in kwargs) - sky = horizon_adjusted(dhi, horizon_profile, surface_tilt, surface_azimuth) + sky = horizon_adjusted(dhi, kwargs["horizon_profile"], surface_tilt, surface_azimuth) else: raise ValueError('invalid model selection {}'.format(model)) @@ -684,7 +684,7 @@ def isotropic(surface_tilt, dhi): ''' sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5 - + # print((1 + tools.cosd(surface_tilt)) * 0.5) return sky_diffuse @@ -1208,25 +1208,52 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, else: return sky_diffuse +def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): + plane_az = np.radians(surface_azimuth) + tilt = np.radians(surface_tilt) + az = np.radians(direction) + + x = np.cos(az) + y = np.sin(az) + z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + np.sin(plane_az) * np.sin(az)) + mask = z <=0 + numer = x*x + y*y + denom = x*x + y*y + z*z + dip = np.degrees(np.arccos(np.sqrt(numer/denom))) + mask = np.logical_or(np.isnan(dip), z <=0) + + if isinstance(dip, pd.Series): + dip[np.isnan(dip)] = 0 + dip[mask] = 0 + else: + dip = np.where(mask, 0, dip) + return dip + -def calculate_dtf(dip_angles, surface_tilt, surface_azimuth): +def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): tilt_rad = np.radians(surface_tilt) plane_az_rad = np.radians(surface_azimuth) a = np.sin(tilt_rad) * np.cos(plane_az_rad) b = np.sin(tilt_rad) * np.sin(plane_az_rad) c = np.cos(tilt_rad) - dtf = 0 - for pair in dip_angles: - az = np.radians(pair[0]) - dip = np.radians(pair[1]) + + dtf = np.zeros(len(surface_tilt)) + for pair in horizon_profile: + az = np.full(len(surface_tilt), np.radians(pair[0])) + horizon_dip = np.full(len(surface_tilt), np.radians(pair[1])) + temp = np.radians(collection_plane_dip_angle(surface_tilt, + surface_azimuth, + pair[0])) + dip = np.maximum(horizon_dip, temp) + #print(dip.shape) first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) second_term = .5 * c * np.cos(dip)**2 dtf += (first_term + second_term) / 180 return dtf def horizon_adjusted(dhi, horizon_profile, surface_tilt, surface_azimuth): - dtf = calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) + sky_diffuse = dhi * dtf return sky_diffuse @@ -2913,3 +2940,20 @@ def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1, (zenith < zenith_threshold_for_zero_dni) & (dni > max_dni)] = max_dni return dni + +def adjust_direct_for_horizon(poa_direct, horizon_profile, + solar_azimuth, solar_zenith): + # assumes horizon_profile is in 1 deg. azimuth increments + dip_angles = {} + for pair in horizon_profile: + dip_angles[pair[0]] = pair[1] + adjusted_irradiance = poa_direct + for i in range(len(poa_direct)): + rounded_solar_azimuth = round(solar_azimuth[i]) + if rounded_solar_azimuth > 180: + rounded_solar_azimuth -= 360 + horizon_dip = dip_angles[rounded_solar_azimuth] + if (solar_zenith[i] > (90 - horizon_dip)): + adjusted_irradiance[i] = 0 + + return adjusted_irradiance \ No newline at end of file diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index a0a6994efc..1690d0253c 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -715,28 +715,20 @@ def no_extra_losses(self): self.losses = 1 return self - def adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith): - # assumes horizon_profile is in 1 deg. azimuth increments - dip_angles = {} - for pair in horizon_profile: - dip_angles[pair[0]] = pair[1] - adjusted_irradiance = poa_direct - for i in range(len(poa_direct)): - rounded_solar_azimuth = round(solar_azimuth[i]) - horizon_dip = dip_angles[rounded_solar_azimuth] - if (solar_zenith > (90 - horizon_dip)): - adjusted_irradiance[i] = 0 - - return adjusted_irradiance - def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) poa_direct = self.total_irrad['poa_direct'] + if self.location.horizon_profile != None: + horizon_profile = self.location.horizon_profile solar_zenith = self.solar_position['apparent_zenith'] solar_azimuth = self.solar_position['azimuth'] - poa_direct = adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith) + poa_direct = pvlib.irradiance.adjust_direct_for_horizon(poa_direct, + horizon_profile, + solar_azimuth, + solar_zenith) + self.effective_irradiance = self.spectral_modifier * ( poa_direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) diff --git a/setup.py b/setup.py index f2e320ccbb..e884be3301 100755 --- a/setup.py +++ b/setup.py @@ -67,8 +67,8 @@ setuptools_kwargs = { 'zip_safe': False, 'scripts': [], - 'include_package_data': True, - 'python_requires': '~=3.5' + 'include_package_data': True +# 'python_requires': '~=3.5' } # set up pvlib packages to be installed and extensions to be compiled From 39302d60d0009c38c50f1f4f926c433fd91d12aa Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 25 Jul 2019 11:08:49 -0700 Subject: [PATCH 03/29] wrapped up most code and added docstrings to irradiance --- pvlib/irradiance.py | 138 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 112 insertions(+), 26 deletions(-) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 9bb680fdef..d34dc7ba41 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -447,7 +447,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model=model_perez) elif model == 'horizon_adjusted': assert("horizon_profile" in kwargs) - sky = horizon_adjusted(dhi, kwargs["horizon_profile"], surface_tilt, surface_azimuth) + sky = horizon_adjusted(surface_tilt, surface_azimuth, dhi, kwargs["horizon_profile"]) else: raise ValueError('invalid model selection {}'.format(model)) @@ -1209,19 +1209,48 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, return sky_diffuse def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): + """ + Determine the dip angle created by the surface of a tilted plane + in a given direction. The dip angle is limited to be non-negative. + + Parameters + ---------- + surface_tilt : numeric + Surface tilt angles in decimal degrees. surface_tilt must be >=0 + and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. surface_azimuth must + be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + direction : numeric + The direction along which the dip angle is to be calculated in + decimal degrees. The convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + Returns + -------- + + dip_angle : numeric + The dip angle in direction created by the tilted plane. + + """ plane_az = np.radians(surface_azimuth) tilt = np.radians(surface_tilt) az = np.radians(direction) x = np.cos(az) y = np.sin(az) - z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + np.sin(plane_az) * np.sin(az)) + z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + \ + np.sin(plane_az) * np.sin(az)) mask = z <=0 numer = x*x + y*y denom = x*x + y*y + z*z dip = np.degrees(np.arccos(np.sqrt(numer/denom))) mask = np.logical_or(np.isnan(dip), z <=0) - + if isinstance(dip, pd.Series): dip[np.isnan(dip)] = 0 dip[mask] = 0 @@ -1231,29 +1260,64 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): + """ + Calculate the diffuse tilt factor that is adjusted with the horizon + profile. + """ tilt_rad = np.radians(surface_tilt) plane_az_rad = np.radians(surface_azimuth) a = np.sin(tilt_rad) * np.cos(plane_az_rad) b = np.sin(tilt_rad) * np.sin(plane_az_rad) c = np.cos(tilt_rad) - dtf = np.zeros(len(surface_tilt)) + # this gets either an int or an array of zeros + dtf = 0.0 * surface_tilt for pair in horizon_profile: - az = np.full(len(surface_tilt), np.radians(pair[0])) - horizon_dip = np.full(len(surface_tilt), np.radians(pair[1])) + az = np.radians(pair[0]) + horizon_dip = np.radians(pair[1]) temp = np.radians(collection_plane_dip_angle(surface_tilt, - surface_azimuth, - pair[0])) + surface_azimuth, + pair[0])) dip = np.maximum(horizon_dip, temp) - #print(dip.shape) + first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) second_term = .5 * c * np.cos(dip)**2 - dtf += (first_term + second_term) / 180 + dtf += 2 * (first_term + second_term) / len(horizon_profile) return dtf -def horizon_adjusted(dhi, horizon_profile, surface_tilt, surface_azimuth): +def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): + ''' + Determine diffuse irradiance from the sky on a tilted surface using + an isotropic model that is adjusted with a horizon profile. + + Parameters + ---------- + surface_tilt : numeric + Surface tilt angles in decimal degrees. surface_tilt must be >=0 + and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. surface_azimuth must + be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + dhi : numeric + Diffuse horizontal irradiance in W/m^2. DHI must be >=0. + + horizon_profile : list + A list of (azimuth, dip_angle) tuples that defines the horizon + profile + + + Returns + -------- + + sky_diffuse : numeric + The sky diffuse component of the solar radiation on a tilted + surface. + ''' dtf = calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) - sky_diffuse = dhi * dtf return sky_diffuse @@ -2943,17 +3007,39 @@ def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1, def adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith): - # assumes horizon_profile is in 1 deg. azimuth increments - dip_angles = {} - for pair in horizon_profile: - dip_angles[pair[0]] = pair[1] - adjusted_irradiance = poa_direct - for i in range(len(poa_direct)): - rounded_solar_azimuth = round(solar_azimuth[i]) - if rounded_solar_azimuth > 180: - rounded_solar_azimuth -= 360 - horizon_dip = dip_angles[rounded_solar_azimuth] - if (solar_zenith[i] > (90 - horizon_dip)): - adjusted_irradiance[i] = 0 - - return adjusted_irradiance \ No newline at end of file + ''' + Adjusts modeled DNI to account for the presence of a horizon. At + each time step check to see if the sun is behind the horizon. If so, + set the DNI to zero. + + Parameters + ---------- + poa_direct : numeric + The modeled direct normal irradiance in the plane of array. + solar_zenith : numeric + Solar zenith angle. + solar_azimuth : numeric + Solar azimuth angle. + horizon_profile : list + A list of (azimuth, dip_angle) tuples that defines the horizon + profile given. Every azimuth from -180 to 180, in 1 degree + increments, must be in this list. + + Returns + ------- + adjusted_irradiance : Series + POA direct normal irradiance after adjusting for the horizon. + ''' + dip_angles = {} + for pair in horizon_profile: + dip_angles[pair[0]] = pair[1] + adjusted_irradiance = poa_direct + for i in range(len(poa_direct)): + rounded_solar_azimuth = round(solar_azimuth[i]) + if rounded_solar_azimuth > 180: + rounded_solar_azimuth -= 360 + horizon_dip = dip_angles[rounded_solar_azimuth] + if (solar_zenith[i] > (90 - horizon_dip)): + adjusted_irradiance[i] = 0 + + return adjusted_irradiance \ No newline at end of file From 2daa5c6533754d9f4256be059612ba7eea9c7f15 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Mon, 29 Jul 2019 11:14:07 -0700 Subject: [PATCH 04/29] added more documentation --- pvlib/horizon.py | 282 ++++++++++++++++++++++++++++++++++++-------- pvlib/irradiance.py | 2 +- 2 files changed, 237 insertions(+), 47 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 4fe91d3eed..6b5fef8e48 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,3 +1,10 @@ +""" +The ``horizon`` module contains functions for horizon profile modeling. +There are various geometric utilities that are useful in horizon calculations +as well as a method that uses the googlemaps elevation API to create a +horizon profile. +""" + import random import pytz import time @@ -14,16 +21,29 @@ def latitude_to_geocentric(phi): + """ + Converts a geodetic (common) latitude to a geocentric latitude. + [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf + """ a = 6378.137 b = 6356.752 return np.arctan(b**2/a**2*np.tan(phi)) def latitude_to_geodetic(phi): + """ + Converts a geocentric latitude to a geodetic (common) latitude. + [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf + """ a = 6378.137 b = 6356.752 return np.arctan(a**2/b**2*np.tan(phi)) def xyz_from_lle(point): + """ + Converts a (lat, lon, elev) tuple into a (x, y, z) tuple. + The center of the earth is the origin in the xyz-space. + The input latitude is assumed to be a common latitude (geodetic). + """ lat = point[0] lon = point[1] elev = point[2] @@ -50,6 +70,11 @@ def xyz_from_lle(point): return v def lle_from_xyz(point): + """ + Converts a (x, y, z) tuple into a (lat, lon, elev) tuple. + The center of the earth is the origin in the xyz-space. + The output latitude is assumed to be a common latitude (geodetic). + """ a = 6378137.0 b = 6356752.0 @@ -75,14 +100,18 @@ def lle_from_xyz(point): return (lat, lon, elev) def pol2cart(rho, phi): + """ + Converts polar coordiantes to cartesian coordinates in 2-d space. + """ x = rho * np.cos(phi) y = rho * np.sin(phi) - return(x, y) + return (x, y) def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): ''' - input lat long - output grid of lat,long tuples + Creates a grid around a location (lat/lon pair) with a specified + grid size and step. The returned grid will be a 3-d matrix with + shape: grid_size+1 x grid_size+1 x 2. ''' grid = np.ndarray((grid_size + 1, grid_size + 1, 2)) @@ -94,16 +123,18 @@ def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): return grid -def get_grid_elevations(in_grid, api_key): - ''' - input grid of lat,lon tuples - output grid of LLE tuples - ''' +def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): + """ + Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). + Queries the googlemaps elevation API to get elevation data at each lat/lon + point. Outputs the original grid with the elevation data appended along + the third axis so the shape is grid_size+1 x grid_size+1 x 3. + """ in_shape = in_grid.shape lats = in_grid.T[0].flatten() longs = in_grid.T[1].flatten() locations = zip(lats, longs) - gmaps = googlemaps.Client(key=api_key) + gmaps = googlemaps.Client(key=GMAPS_API_KEY) out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) @@ -134,6 +165,30 @@ def dip_calc(pt1, pt2): ''' input: two LLE tuples output: distance, dip angle, azimuth + + Calculates the dip angle from pt1 to pt2 where dip angle is defined as + the angle that the line connecting pt1 to pt2 makes with the plane normal + to the Earth's surface at pt2. Also computes the azimuth defined as degrees + East of North the bearing of pt2 from pt1. This uses the Haversine formula. + + Parameters + ---------- + pt1 : tuple + (lat, lon, elev) tuple that corresponds to the origin from which + the dip angle is to be calculated. The observer point. + + pt1 : tuple + (lat, lon, elev) tuple that corresponds to the origin from which + the dip angle is to be calculated. The observee point. + + Returns + ------- + bearing_deg: numeric + The bearing from pt1 to pt2 in degrees East of North + + dip_angle: + The dip angle that pt2 makes with the horizontal as observed at pt1. + ''' a = 6378137.0 b = 6356752.0 @@ -170,18 +225,37 @@ def dip_calc(pt1, pt2): dip_angle_deg = dip_angle*180.0/np.pi # might wanna double check this formula (haversine?) - azimuth = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) - azimuth_deg = azimuth*180.0/np.pi + bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) + bearing_deg = bearing*180.0/np.pi - return (azimuth_deg, dip_angle_deg) + return (bearing_deg, dip_angle_deg) def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): - ''' - input grid of lat,lon,elevation tuples - output list of azimuth, distance, dip angle - use grid points in first pass and then try randomly sampling triangle method - ''' + """ + Calculates a horizon profile from a grid of (lat, lon, elev) tuples. + The "site" is assumed to be at the center of the grid. + + Parameters + ---------- + grid : ndarray + Assumes + + sampling_method : string + A string that specifies the sampling method used to generate the + horizon profile. + + sampling_param : variable + A parameter that is passed into the function specified by + sampling_method. + + Returns + ------- + horizon_points: list + List of (azimuth, dip_angle) tuples that define the horizon at the + point at the center of the grid. + + """ grid_shape = grid.shape grid_center_i = (grid_shape[0] - 1) / 2 @@ -192,26 +266,26 @@ def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): site = (site_lat, site_lon, site_elev) horizon_points = [] - start = time.time() + if sampling_method == "grid": samples = sample_using_grid(grid) elif sampling_method == "triangles": samples = sample_using_triangles(grid, sampling_param) elif sampling_method == "interpolator": samples = sample_using_interpolator(grid, sampling_param) - post_sampling = time.time() - print("Sampling took " + str(post_sampling-start) + " sec") - + dip_calc_lambda = lambda pt: dip_calc(site, pt) horizon_points = np.array(list(map(dip_calc_lambda, samples))) - post_calcs = time.time() - print("Dip calcs on samples took " + str(post_calcs-post_sampling) + " sec") return horizon_points def sample_using_grid(grid): - # use every grid point as a sample - # returns a list of LLE tuples + """ + Calculates the dip angle from the site (center of the grid) + to every point on the grid and uses the results as the + horizon profile. + + """ grid_shape = grid.shape grid_center_i = (grid_shape[0] - 1) / 2 grid_center_j = (grid_shape[1] - 1) / 2 @@ -228,8 +302,26 @@ def sample_using_grid(grid): return all_samples def sample_using_triangles(grid, samples_per_triangle=10): - # uniformly sample all triangles between neighboring grid points - # returns a list of LLE tuples + """ + Creates triangles using nearest neighbors for every grid point and randomly + samples each of these triangles to find dip angles for the horizon profile. + + Parameters + ---------- + grid : ndarray + Grid that contains lat lon and elev information. + + samples_per_triangle : numeric + The number of random samples to be uniformly taken from the surface + of each triangle. + + Returns + ------- + all_samples: list + List of (azimuth, dip_angle) tuples that were sampled from the grid + + [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf + """ all_samples = [] for i in range(grid.shape[0]): @@ -238,25 +330,61 @@ def sample_using_triangles(grid, samples_per_triangle=10): if i != 0 and j != 0: left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) - all_samples += uniformly_sample_triangle(center, top, left , samples_per_triangle) + all_samples += uniformly_sample_triangle(center, + top, + left, + samples_per_triangle) if i != 0 and j != grid.shape[1] - 1: right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) - all_samples += uniformly_sample_triangle(center, top, right , samples_per_triangle) + all_samples += uniformly_sample_triangle(center, + top, + right, + samples_per_triangle) if i != grid.shape[0] - 1 and j != 0: left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) - all_samples += uniformly_sample_triangle(center, bottom, left , samples_per_triangle) + all_samples += uniformly_sample_triangle(center, + bottom, + left, + samples_per_triangle) if i != grid.shape[0] - 1 and j != grid.shape[1] - 1: right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) - all_samples += uniformly_sample_triangle(center, bottom, right , samples_per_triangle) + all_samples += uniformly_sample_triangle(center, + bottom, + right, + samples_per_triangle) return all_samples def sample_using_interpolator(grid, num_samples): + """ + Creates a "grid" using polar coordinates and uses the scipy's grid + interpolator to estimate elevation values at each point on the polar grid + from the input (rectangular) grid that has true elevation values. Dip + calculations are done at each point on the polar grid and the results + are returned. + + Parameters + ---------- + grid : ndarray + Grid that contains lat lon and elev information. + + num_samples : tuple + A tuple containing two integers. The first is the desired number of + points along the radial axis of the polar grid. The second is the + desired number of points along the angular axis of the polar grid. + + + Returns + ------- + all_samples: list + List of (azimuth, dip_angle) tuples taken from the polar grid + + """ x = grid.T[0][0] y = grid.T[1].T[0] @@ -288,7 +416,31 @@ def sample_using_interpolator(grid, num_samples): def uniformly_sample_triangle(p1, p2, p3, num_samples): - # returns a list of LLE tuples + """ + Randomly sample the surface of a triangle defined by three (lat, lon, elev) + points uniformly [1]. + + Parameters + ---------- + pt1 : tuple + A (lat, lon, elev) tuple that defines one vertex of the triangle. + pt2 : tuple + A (lat, lon, elev) tuple that defines another vertex of the triangle. + pt3 : tuple + A (lat, lon, elev) tuple that defines the last vertex of the triangle. + + num_samples : tuple + The number of random samples to be uniformly taken from the surface + of the triangle. + + Returns + ------- + points: list + List of (lat, lon, elev) tuples that lie on the surface of the + triangle. + + [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf + """ c1 = xyz_from_lle(p1) c2 = xyz_from_lle(p2) c3 = xyz_from_lle(p3) @@ -299,20 +451,41 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): r2 = np.random.rand() sqrt_r1 = np.sqrt(r1) - # use uniform sampling from http://www.sherrytowers.com/randomly_sample_points_in_triangle.pdf random_pt = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 points.append(lle_from_xyz(random_pt)) return points def round_to_nearest(x, base): + """ + Helper function to round x to nearest base. + """ return base * round(float(x) / base) -def filter_points(horizon_points, bucket_size=1): +def filter_points(horizon_points, bin_size=1): + """ + Bins the horizon_points by azimuth values. The azimuth value of each + point in horizon_points is rounded to the nearest bin and then the + max value in each bin is returned. + + Parameters + ---------- + horizon_points : list + List of (azimuth, dip_angle) tuples that define the horizon. + + bin_size : int + The bin size of azimuth values. + + Returns + ------- + sorted_points: list + List of (azimuth, dip_angle) values that correspond to the greatest + dip_angle in each azimuth bin. + """ wedges = {} for pair in horizon_points: azimuth = pair[0] dip = pair[1] - azimuth_wedge = round_to_nearest(azimuth, bucket_size) + azimuth_wedge = round_to_nearest(azimuth, bin_size) if azimuth_wedge in wedges: wedges[azimuth_wedge] = max(dip, wedges[azimuth_wedge]) @@ -327,7 +500,9 @@ def filter_points(horizon_points, bucket_size=1): return sorted_points def visualize(horizon_profile, pvsyst_scale=False): - + """ + Plots a horizon profile with azimuth on the x-axis and dip angle on the y. + """ azimuths = [] dips = [] for pair in horizon_profile: @@ -342,7 +517,11 @@ def visualize(horizon_profile, pvsyst_scale=False): def polar_plot(horizon_profile): - + """ + Plots a horizon profile on a polar plot with dip angle as the raidus and + azimuth as the theta value. An offset of 5 is added to the dip_angle to + make the plot more readable with low dip angles. + """ azimuths = [] dips = [] for pair in horizon_profile: @@ -356,12 +535,13 @@ def polar_plot(horizon_profile): plt.plot(azimuths, dips, "o") plt.show -def horizon_table(horizon_points): - for pair in horizon_points: - azimuth = pair[0] - print(str(azimuth) + ": " + str(pair[1])) def invert_for_pvsyst(horizon_points, hemisphere="north"): + """ + Modify the azimuth values in horizon_points to match PVSyst's azimuth + convention (which is dependent on hemisphere) + """ + # look at that northern hemisphere bias right there # not even sorry. assert hemisphere == "north" or hemisphere == "south" @@ -380,16 +560,26 @@ def invert_for_pvsyst(horizon_points, hemisphere="north"): return sorted_points def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): + """ + Uses the functions defined in this modules to generate a complete horizon + profile for a location (specified by lat/lon). An API key for the + googlemaps elevation API is needeed. + """ grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) - elev_grid = get_grid_elevations(grid, GMAPS_API_KEY) - horizon_points = calculate_horizon_points(elev_grid, sampling_method="interpolator", sampling_param=(1000,1000)) - filtered_points = filter_points(horizon_points, bucket_size=1) + elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) + horizon_points = calculate_horizon_points(elev_grid, + sampling_method="interpolator", + sampling_param=(1000,1000)) + filtered_points = filter_points(horizon_points, bin_size=1) return filtered_points def fake_horizon_profile(max_dip): + """ + Creates a bogus horizon profile by randomly generating dip_angles at + integral azimuth values. Used for testing purposes. + """ fake_profile = [] for i in range(-180, 181): - fake_profile.append((i, random.random()*max_dip)) - + fake_profile.append((i, random.random() * max_dip)) return fake_profile \ No newline at end of file diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index d34dc7ba41..78b280d62e 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -684,7 +684,7 @@ def isotropic(surface_tilt, dhi): ''' sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5 - # print((1 + tools.cosd(surface_tilt)) * 0.5) + return sky_diffuse From 6cb73992660fb8b6b34a24826118dab3d1380fa0 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Mon, 29 Jul 2019 11:16:21 -0700 Subject: [PATCH 05/29] reverted setup --- setup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index e884be3301..f2e320ccbb 100755 --- a/setup.py +++ b/setup.py @@ -67,8 +67,8 @@ setuptools_kwargs = { 'zip_safe': False, 'scripts': [], - 'include_package_data': True -# 'python_requires': '~=3.5' + 'include_package_data': True, + 'python_requires': '~=3.5' } # set up pvlib packages to be installed and extensions to be compiled From df724e5a0a05f91d3187daa85b91a9b8ea082c82 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Mon, 29 Jul 2019 12:56:41 -0700 Subject: [PATCH 06/29] linted --- pvlib/horizon.py | 220 ++++++++++++++++++++++---------------------- pvlib/irradiance.py | 28 +++--- pvlib/modelchain.py | 19 ++-- 3 files changed, 137 insertions(+), 130 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 6b5fef8e48..7a4fb4630f 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,19 +1,15 @@ """ The ``horizon`` module contains functions for horizon profile modeling. There are various geometric utilities that are useful in horizon calculations -as well as a method that uses the googlemaps elevation API to create a +as well as a method that uses the googlemaps elevation API to create a horizon profile. """ import random -import pytz -import time import itertools import numpy as np -import scipy as sp from scipy.interpolate import RegularGridInterpolator -import pandas as pd import matplotlib.pyplot as plt @@ -29,6 +25,7 @@ def latitude_to_geocentric(phi): b = 6356.752 return np.arctan(b**2/a**2*np.tan(phi)) + def latitude_to_geodetic(phi): """ Converts a geocentric latitude to a geodetic (common) latitude. @@ -38,6 +35,7 @@ def latitude_to_geodetic(phi): b = 6356.752 return np.arctan(a**2/b**2*np.tan(phi)) + def xyz_from_lle(point): """ Converts a (lat, lon, elev) tuple into a (x, y, z) tuple. @@ -47,19 +45,19 @@ def xyz_from_lle(point): lat = point[0] lon = point[1] elev = point[2] - + a = 6378137.0 b = 6356752.0 - + # convert to radians phi = lat*np.pi/180.0 theta = lon*np.pi/180.0 - + # compute radius of earth at each point r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 r = r / (a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2) r = np.sqrt(r) - + h = r + elev alpha = latitude_to_geocentric(phi) beta = theta @@ -69,6 +67,7 @@ def xyz_from_lle(point): v = np.array((x, y, z)) return v + def lle_from_xyz(point): """ Converts a (x, y, z) tuple into a (lat, lon, elev) tuple. @@ -77,28 +76,27 @@ def lle_from_xyz(point): """ a = 6378137.0 b = 6356752.0 - + x = point[0] y = point[1] z = point[2] - + # get corresponding point on earth's surface t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) point_s = t * point - x_s = point_s[0] - y_s = point_s[1] z_s = point_s[2] - + elev = np.linalg.norm(point-point_s) r = np.linalg.norm(point_s) - + alpha = np.arcsin(z_s / r) phi = latitude_to_geodetic(alpha) lat = phi*180.0/np.pi - + lon = np.arctan2(y, x)*180/np.pi return (lat, lon, elev) + def pol2cart(rho, phi): """ Converts polar coordiantes to cartesian coordinates in 2-d space. @@ -107,6 +105,7 @@ def pol2cart(rho, phi): y = rho * np.sin(phi) return (x, y) + def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): ''' Creates a grid around a location (lat/lon pair) with a specified @@ -114,18 +113,19 @@ def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): shape: grid_size+1 x grid_size+1 x 2. ''' grid = np.ndarray((grid_size + 1, grid_size + 1, 2)) - + # fill out grid for i in range(grid_size + 1): for j in range(grid_size + 1): - grid[i,j,0] = lat + (i - grid_size / 2) * grid_step - grid[i,j,1] = lon + (j - grid_size / 2) * grid_step - + grid[i, j, 0] = lat + (i - grid_size / 2) * grid_step + grid[i, j, 1] = lon + (j - grid_size / 2) * grid_step + return grid - + + def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): """ - Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). + Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). Queries the googlemaps elevation API to get elevation data at each lat/lon point. Outputs the original grid with the elevation data appended along the third axis so the shape is grid_size+1 x grid_size+1 x 3. @@ -135,13 +135,13 @@ def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): longs = in_grid.T[1].flatten() locations = zip(lats, longs) gmaps = googlemaps.Client(key=GMAPS_API_KEY) - + out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) - + # Get elevation data from gmaps elevations = [] responses = [] - + while len(locations) > 512: locations_to_request = locations[:512] locations = locations[512:] @@ -149,27 +149,28 @@ def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): responses += gmaps.elevation(locations=locations) for entry in responses: elevations.append(entry["elevation"]) - + for i in range(in_shape[0]): for j in range(in_shape[1]): - lat = in_grid[i,j,0] - lon = in_grid[i,j,1] + lat = in_grid[i, j, 0] + lon = in_grid[i, j, 1] elevation = elevations[i + j * in_shape[1]] - - out_grid[i,j,0] = lat - out_grid[i,j,1] = lon - out_grid[i,j,2] = elevation + + out_grid[i, j, 0] = lat + out_grid[i, j, 1] = lon + out_grid[i, j, 2] = elevation return out_grid - + + def dip_calc(pt1, pt2): ''' input: two LLE tuples output: distance, dip angle, azimuth - Calculates the dip angle from pt1 to pt2 where dip angle is defined as + Calculates the dip angle from pt1 to pt2 where dip angle is defined as the angle that the line connecting pt1 to pt2 makes with the plane normal to the Earth's surface at pt2. Also computes the azimuth defined as degrees - East of North the bearing of pt2 from pt1. This uses the Haversine formula. + East of North the bearing of pt2 from pt1. This uses the Haversine formula. Parameters ---------- @@ -186,51 +187,50 @@ def dip_calc(pt1, pt2): bearing_deg: numeric The bearing from pt1 to pt2 in degrees East of North - dip_angle: + dip_angle: The dip angle that pt2 makes with the horizontal as observed at pt1. - + ''' a = 6378137.0 b = 6356752.0 - + lat1 = pt1[0] lon1 = pt1[1] elev1 = pt1[2] lat2 = pt2[0] lon2 = pt2[1] elev2 = pt2[2] - + # convert to radians phi1 = lat1*np.pi/180.0 theta1 = lon1*np.pi/180.0 phi2 = lat2*np.pi/180.0 theta2 = lon2*np.pi/180.0 - + v1 = xyz_from_lle((lat1, lon1, elev1)) v2 = xyz_from_lle((lat2, lon2, elev2)) - + x1 = v1[0] y1 = v1[1] z1 = v1[2] - x2 = v2[0] - y2 = v2[1] - z2 = v2[2] - - delta = np.subtract(v1,v2) - distance = np.linalg.norm(delta) - + + delta = np.subtract(v1, v2) + normal = np.array((2*x1/a**2, 2*y1/a**2, 2*z1/b**2)) - beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta)/np.linalg.norm(normal)) + beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta) / + np.linalg.norm(normal)) dip_angle = beta - np.pi/2 dip_angle_deg = dip_angle*180.0/np.pi # might wanna double check this formula (haversine?) - bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) + bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), + np.cos(phi1) * np.sin(phi2) - + np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) bearing_deg = bearing*180.0/np.pi - - + return (bearing_deg, dip_angle_deg) - + + def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): """ Calculates a horizon profile from a grid of (lat, lon, elev) tuples. @@ -253,10 +253,10 @@ def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): ------- horizon_points: list List of (azimuth, dip_angle) tuples that define the horizon at the - point at the center of the grid. + point at the center of the grid. """ - + grid_shape = grid.shape grid_center_i = (grid_shape[0] - 1) / 2 grid_center_j = (grid_shape[1] - 1) / 2 @@ -264,7 +264,7 @@ def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): site_lon = grid[grid_center_i, grid_center_j, 1] site_elev = grid[grid_center_i, grid_center_j, 2] site = (site_lat, site_lon, site_elev) - + horizon_points = [] if sampling_method == "grid": @@ -274,11 +274,12 @@ def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): elif sampling_method == "interpolator": samples = sample_using_interpolator(grid, sampling_param) - dip_calc_lambda = lambda pt: dip_calc(site, pt) - horizon_points = np.array(list(map(dip_calc_lambda, samples))) + horizon_points = np.array(list(map(lambda pt: dip_calc(site, pt), + samples))) return horizon_points + def sample_using_grid(grid): """ Calculates the dip angle from the site (center of the grid) @@ -289,7 +290,7 @@ def sample_using_grid(grid): grid_shape = grid.shape grid_center_i = (grid_shape[0] - 1) / 2 grid_center_j = (grid_shape[1] - 1) / 2 - + all_samples = [] for i in range(grid_shape[0]): for j in range(grid_shape[1]): @@ -301,6 +302,7 @@ def sample_using_grid(grid): all_samples.append((lat, lon, elev)) return all_samples + def sample_using_triangles(grid, samples_per_triangle=10): """ Creates triangles using nearest neighbors for every grid point and randomly @@ -309,7 +311,7 @@ def sample_using_triangles(grid, samples_per_triangle=10): Parameters ---------- grid : ndarray - Grid that contains lat lon and elev information. + Grid that contains lat lon and elev information. samples_per_triangle : numeric The number of random samples to be uniformly taken from the surface @@ -322,44 +324,45 @@ def sample_using_triangles(grid, samples_per_triangle=10): [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf """ - + all_samples = [] for i in range(grid.shape[0]): for j in range(grid.shape[1]): - center = (grid[i,j,0], grid[i,j,1], grid[i,j,2]) + center = (grid[i, j, 0], grid[i, j, 1], grid[i, j, 2]) if i != 0 and j != 0: - left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) - top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) + left = (grid[i, j-1, 0], grid[i, j-1, 1], grid[i, j-1, 2]) + top = (grid[i-1, j, 0], grid[i-1, j, 1], grid[i-1, j, 2]) all_samples += uniformly_sample_triangle(center, top, left, samples_per_triangle) - + if i != 0 and j != grid.shape[1] - 1: - right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) - top = (grid[i-1,j,0], grid[i-1,j,1], grid[i-1,j,2]) + right = (grid[i, j+1, 0], grid[i, j+1, 1], grid[i, j+1, 2]) + top = (grid[i-1, j, 0], grid[i-1, j, 1], grid[i-1, j, 2]) all_samples += uniformly_sample_triangle(center, top, right, samples_per_triangle) - + if i != grid.shape[0] - 1 and j != 0: - left = (grid[i,j-1,0], grid[i,j-1,1], grid[i,j-1,2]) - bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) + left = (grid[i, j-1, 0], grid[i, j-1, 1], grid[i, j-1, 2]) + bottom = (grid[i+1, j, 0], grid[i+1, j, 1], grid[i+1, j, 2]) all_samples += uniformly_sample_triangle(center, bottom, left, samples_per_triangle) - + if i != grid.shape[0] - 1 and j != grid.shape[1] - 1: - right = (grid[i,j+1,0], grid[i,j+1,1], grid[i,j+1,2]) - bottom = (grid[i+1,j,0], grid[i+1,j,1], grid[i+1,j,2]) + right = (grid[i, j+1, 0], grid[i, j+1, 1], grid[i, j+1, 2]) + bottom = (grid[i+1, j, 0], grid[i+1, j, 1], grid[i+1, j, 2]) all_samples += uniformly_sample_triangle(center, bottom, right, samples_per_triangle) return all_samples + def sample_using_interpolator(grid, num_samples): """ Creates a "grid" using polar coordinates and uses the scipy's grid @@ -371,7 +374,7 @@ def sample_using_interpolator(grid, num_samples): Parameters ---------- grid : ndarray - Grid that contains lat lon and elev information. + Grid that contains lat lon and elev information. num_samples : tuple A tuple containing two integers. The first is the desired number of @@ -387,12 +390,9 @@ def sample_using_interpolator(grid, num_samples): """ x = grid.T[0][0] y = grid.T[1].T[0] - - x_base = x[0] + x_range = x[-1] - x[0] - y_base = y[0] - y_range = y[-1] - y[0] - + grid_shape = grid.shape grid_center_i = (grid_shape[0] - 1) / 2 grid_center_j = (grid_shape[1] - 1) / 2 @@ -400,20 +400,20 @@ def sample_using_interpolator(grid, num_samples): site_lon = grid[grid_center_i, grid_center_j, 1] elevs = grid.T[2].T - interpolator = RegularGridInterpolator((x,y), elevs) - + interpolator = RegularGridInterpolator((x, y), elevs) r = np.linspace(0, x_range/2, num_samples[0]) theta = np.linspace(0, 2 * np.pi, num_samples[1]) polar_pts = np.array(list(itertools.product(r, theta))) - + pts = np.array([pol2cart(e[0], e[1]) for e in polar_pts]) pts += np.array((site_lat, site_lon)) - - interpolated_elevs = interpolator(pts).reshape(num_samples[0]*num_samples[1], 1) + total_num_samples = num_samples[0]*num_samples[1] + + interpolated_elevs = interpolator(pts).reshape(total_num_samples, 1) samples = np.concatenate((pts, interpolated_elevs), axis=1) return samples - + def uniformly_sample_triangle(p1, p2, p3, num_samples): """ @@ -427,7 +427,7 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): pt2 : tuple A (lat, lon, elev) tuple that defines another vertex of the triangle. pt3 : tuple - A (lat, lon, elev) tuple that defines the last vertex of the triangle. + A (lat, lon, elev) tuple that defines the last vertex of the triangle. num_samples : tuple The number of random samples to be uniformly taken from the surface @@ -436,7 +436,7 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): Returns ------- points: list - List of (lat, lon, elev) tuples that lie on the surface of the + List of (lat, lon, elev) tuples that lie on the surface of the triangle. [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf @@ -444,28 +444,30 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): c1 = xyz_from_lle(p1) c2 = xyz_from_lle(p2) c3 = xyz_from_lle(p3) - + points = [] for i in range(num_samples): r1 = np.random.rand() r2 = np.random.rand() sqrt_r1 = np.sqrt(r1) - + random_pt = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 points.append(lle_from_xyz(random_pt)) return points + def round_to_nearest(x, base): """ Helper function to round x to nearest base. """ return base * round(float(x) / base) + def filter_points(horizon_points, bin_size=1): """ Bins the horizon_points by azimuth values. The azimuth value of each point in horizon_points is rounded to the nearest bin and then the - max value in each bin is returned. + max value in each bin is returned. Parameters ---------- @@ -473,7 +475,7 @@ def filter_points(horizon_points, bin_size=1): List of (azimuth, dip_angle) tuples that define the horizon. bin_size : int - The bin size of azimuth values. + The bin size of azimuth values. Returns ------- @@ -486,19 +488,20 @@ def filter_points(horizon_points, bin_size=1): azimuth = pair[0] dip = pair[1] azimuth_wedge = round_to_nearest(azimuth, bin_size) - + if azimuth_wedge in wedges: wedges[azimuth_wedge] = max(dip, wedges[azimuth_wedge]) else: wedges[azimuth_wedge] = dip - + filtered_points = [] for key in wedges.keys(): filtered_points.append((key, wedges[key])) - + sorted_points = sorted(filtered_points, key=lambda tup: tup[0]) return sorted_points - + + def visualize(horizon_profile, pvsyst_scale=False): """ Plots a horizon profile with azimuth on the x-axis and dip angle on the y. @@ -509,33 +512,33 @@ def visualize(horizon_profile, pvsyst_scale=False): azimuth = pair[0] azimuths.append(azimuth) dips.append(pair[1]) - plt.figure(figsize=(10,6)) + plt.figure(figsize=(10, 6)) if pvsyst_scale: plt.ylim(0, 90) plt.plot(azimuths, dips, "-") plt.show - + def polar_plot(horizon_profile): """ Plots a horizon profile on a polar plot with dip angle as the raidus and azimuth as the theta value. An offset of 5 is added to the dip_angle to make the plot more readable with low dip angles. - """ + """ azimuths = [] dips = [] for pair in horizon_profile: azimuth = pair[0] azimuths.append(np.radians(azimuth)) dips.append(pair[1] + 5) - plt.figure(figsize=(10,6)) + plt.figure(figsize=(10, 6)) sp = plt.subplot(1, 1, 1, projection='polar') sp.set_theta_zero_location('N') sp.set_theta_direction(-1) plt.plot(azimuths, dips, "o") plt.show - + def invert_for_pvsyst(horizon_points, hemisphere="north"): """ Modify the azimuth values in horizon_points to match PVSyst's azimuth @@ -545,7 +548,7 @@ def invert_for_pvsyst(horizon_points, hemisphere="north"): # look at that northern hemisphere bias right there # not even sorry. assert hemisphere == "north" or hemisphere == "south" - + inverted_points = [] for pair in horizon_points: azimuth = pair[0] @@ -559,19 +562,20 @@ def invert_for_pvsyst(horizon_points, hemisphere="north"): sorted_points = sorted(inverted_points, key=lambda tup: tup[0]) return sorted_points + def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): """ Uses the functions defined in this modules to generate a complete horizon profile for a location (specified by lat/lon). An API key for the googlemaps elevation API is needeed. """ - grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) - elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) - horizon_points = calculate_horizon_points(elev_grid, + grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) + elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) + horizon_points = calculate_horizon_points(elev_grid, sampling_method="interpolator", - sampling_param=(1000,1000)) - filtered_points = filter_points(horizon_points, bin_size=1) - return filtered_points + sampling_param=(1000, 1000)) + filtered_points = filter_points(horizon_points, bin_size=1) + return filtered_points def fake_horizon_profile(max_dip): @@ -582,4 +586,4 @@ def fake_horizon_profile(max_dip): fake_profile = [] for i in range(-180, 181): fake_profile.append((i, random.random() * max_dip)) - return fake_profile \ No newline at end of file + return fake_profile diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 78b280d62e..0ad7191fd1 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -447,7 +447,8 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model=model_perez) elif model == 'horizon_adjusted': assert("horizon_profile" in kwargs) - sky = horizon_adjusted(surface_tilt, surface_azimuth, dhi, kwargs["horizon_profile"]) + sky = horizon_adjusted(surface_tilt, surface_azimuth, dhi, + kwargs["horizon_profile"]) else: raise ValueError('invalid model selection {}'.format(model)) @@ -1208,6 +1209,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, else: return sky_diffuse + def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): """ Determine the dip angle created by the surface of a tilted plane @@ -1243,13 +1245,13 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): x = np.cos(az) y = np.sin(az) - z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + \ + z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + np.sin(plane_az) * np.sin(az)) - mask = z <=0 + mask = z <= 0 numer = x*x + y*y denom = x*x + y*y + z*z dip = np.degrees(np.arccos(np.sqrt(numer/denom))) - mask = np.logical_or(np.isnan(dip), z <=0) + mask = np.logical_or(np.isnan(dip), z <= 0) if isinstance(dip, pd.Series): dip[np.isnan(dip)] = 0 @@ -1262,29 +1264,31 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): """ Calculate the diffuse tilt factor that is adjusted with the horizon - profile. + profile. """ tilt_rad = np.radians(surface_tilt) plane_az_rad = np.radians(surface_azimuth) a = np.sin(tilt_rad) * np.cos(plane_az_rad) b = np.sin(tilt_rad) * np.sin(plane_az_rad) c = np.cos(tilt_rad) - + # this gets either an int or an array of zeros dtf = 0.0 * surface_tilt for pair in horizon_profile: az = np.radians(pair[0]) horizon_dip = np.radians(pair[1]) temp = np.radians(collection_plane_dip_angle(surface_tilt, - surface_azimuth, - pair[0])) + surface_azimuth, + pair[0])) dip = np.maximum(horizon_dip, temp) - first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) + first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * \ + (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) second_term = .5 * c * np.cos(dip)**2 dtf += 2 * (first_term + second_term) / len(horizon_profile) return dtf + def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): ''' Determine diffuse irradiance from the sky on a tilted surface using @@ -1323,7 +1327,6 @@ def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): return sky_diffuse - def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0): """ Calculate the clearsky index. @@ -3005,6 +3008,7 @@ def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1, (dni > max_dni)] = max_dni return dni + def adjust_direct_for_horizon(poa_direct, horizon_profile, solar_azimuth, solar_zenith): ''' @@ -3015,7 +3019,7 @@ def adjust_direct_for_horizon(poa_direct, horizon_profile, Parameters ---------- poa_direct : numeric - The modeled direct normal irradiance in the plane of array. + The modeled direct normal irradiance in the plane of array. solar_zenith : numeric Solar zenith angle. solar_azimuth : numeric @@ -3042,4 +3046,4 @@ def adjust_direct_for_horizon(poa_direct, horizon_profile, if (solar_zenith[i] > (90 - horizon_dip)): adjusted_irradiance[i] = 0 - return adjusted_irradiance \ No newline at end of file + return adjusted_irradiance diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index 1690d0253c..ff91cadab0 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -715,22 +715,21 @@ def no_extra_losses(self): self.losses = 1 return self - def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) poa_direct = self.total_irrad['poa_direct'] - if self.location.horizon_profile != None: - horizon_profile = self.location.horizon_profile + if self.location.horizon_profile is not None: + horizon = self.location.horizon_profile solar_zenith = self.solar_position['apparent_zenith'] - solar_azimuth = self.solar_position['azimuth'] - poa_direct = pvlib.irradiance.adjust_direct_for_horizon(poa_direct, - horizon_profile, - solar_azimuth, - solar_zenith) - + solar_az = self.solar_position['azimuth'] + adjusted = pvlib.irradiance.adjust_direct_for_horizon(poa_direct, + horizon, + solar_az, + solar_zenith) + self.effective_irradiance = self.spectral_modifier * ( - poa_direct*self.aoi_modifier + + adjusted*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) return self From 808af03626f7c70c8b3d9e10d02bd43963a2177c Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Mon, 29 Jul 2019 13:40:41 -0700 Subject: [PATCH 07/29] more lints + moved import of gmaps --- pvlib/horizon.py | 14 ++++++++------ pvlib/irradiance.py | 4 ++-- pvlib/location.py | 2 -- pvlib/modelchain.py | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 7a4fb4630f..15af03cbeb 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -13,8 +13,6 @@ import matplotlib.pyplot as plt -import googlemaps - def latitude_to_geocentric(phi): """ @@ -130,6 +128,9 @@ def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): point. Outputs the original grid with the elevation data appended along the third axis so the shape is grid_size+1 x grid_size+1 x 3. """ + + import googlemaps + in_shape = in_grid.shape lats = in_grid.T[0].flatten() longs = in_grid.T[1].flatten() @@ -217,15 +218,15 @@ def dip_calc(pt1, pt2): delta = np.subtract(v1, v2) normal = np.array((2*x1/a**2, 2*y1/a**2, 2*z1/b**2)) - beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta) / - np.linalg.norm(normal)) + beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta) + / np.linalg.norm(normal)) dip_angle = beta - np.pi/2 dip_angle_deg = dip_angle*180.0/np.pi # might wanna double check this formula (haversine?) bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), - np.cos(phi1) * np.sin(phi2) - - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) + np.cos(phi1) * np.sin(phi2) + - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) bearing_deg = bearing*180.0/np.pi return (bearing_deg, dip_angle_deg) @@ -569,6 +570,7 @@ def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): profile for a location (specified by lat/lon). An API key for the googlemaps elevation API is needeed. """ + grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) horizon_points = calculate_horizon_points(elev_grid, diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 0ad7191fd1..210a1fc4fe 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -1245,8 +1245,8 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): x = np.cos(az) y = np.sin(az) - z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + - np.sin(plane_az) * np.sin(az)) + z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + + np.sin(plane_az) * np.sin(az)) mask = z <= 0 numer = x*x + y*y denom = x*x + y*y + z*z diff --git a/pvlib/location.py b/pvlib/location.py index ca0753fb83..57cd77d727 100644 --- a/pvlib/location.py +++ b/pvlib/location.py @@ -317,5 +317,3 @@ def get_sun_rise_set_transit(self, times, method='pyephem', **kwargs): 'one of pyephem, spa, geometric' .format(method)) return result - - \ No newline at end of file diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index ff91cadab0..432541dbe5 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -729,8 +729,8 @@ def effective_irradiance_model(self): solar_zenith) self.effective_irradiance = self.spectral_modifier * ( - adjusted*self.aoi_modifier + - fd*self.total_irrad['poa_diffuse']) + adjusted*self.aoi_modifier + + fd*self.total_irrad['poa_diffuse']) return self def complete_irradiance(self, times=None, weather=None): From 757376358e148785596ed375ca33e535fe0b7ba7 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 31 Jul 2019 15:00:13 -0700 Subject: [PATCH 08/29] added scipy to setup.py and fixed bug in modelchain --- pvlib/horizon.py | 4 ++-- pvlib/modelchain.py | 12 ++++++------ setup.py | 3 ++- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 15af03cbeb..c408c9e0dd 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -225,8 +225,8 @@ def dip_calc(pt1, pt2): # might wanna double check this formula (haversine?) bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), - np.cos(phi1) * np.sin(phi2) - - np.sin(phi1) * np.cos(phi2) * np.cos(theta2-theta1)) + (np.cos(phi1) * np.sin(phi2) + - np.sin(phi1) * np.cos(phi2)*np.cos(theta2-theta1))) bearing_deg = bearing*180.0/np.pi return (bearing_deg, dip_angle_deg) diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index 432541dbe5..4e87c77fef 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -717,19 +717,19 @@ def no_extra_losses(self): def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) - poa_direct = self.total_irrad['poa_direct'] + direct = self.total_irrad['poa_direct'] if self.location.horizon_profile is not None: horizon = self.location.horizon_profile solar_zenith = self.solar_position['apparent_zenith'] solar_az = self.solar_position['azimuth'] - adjusted = pvlib.irradiance.adjust_direct_for_horizon(poa_direct, - horizon, - solar_az, - solar_zenith) + direct = pvlib.irradiance.adjust_direct_for_horizon(direct, + horizon, + solar_az, + solar_zenith) self.effective_irradiance = self.spectral_modifier * ( - adjusted*self.aoi_modifier + direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) return self diff --git a/setup.py b/setup.py index f2e320ccbb..ac55de8c55 100755 --- a/setup.py +++ b/setup.py @@ -40,7 +40,8 @@ INSTALL_REQUIRES = ['numpy >= 1.10.4', 'pandas >= 0.18.1', 'pytz', - 'requests'] + 'requests', + 'scipy'] TESTS_REQUIRE = ['nose', 'pytest', 'pytest-cov', 'pytest-mock', 'pytest-timeout'] EXTRAS_REQUIRE = { From 34218996f9738d9585d306e58bec3ef36d70eebd Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 1 Aug 2019 11:03:16 -0700 Subject: [PATCH 09/29] moved some horizon functions to tools and fixed naming. Also wrote 2 tests for new irradiance functions. --- pvlib/horizon.py | 114 +++------------------------------- pvlib/irradiance.py | 2 +- pvlib/modelchain.py | 3 +- pvlib/test/test_irradiance.py | 48 ++++++++++++++ pvlib/tools.py | 97 +++++++++++++++++++++++++++++ setup.py | 2 +- 6 files changed, 157 insertions(+), 109 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index c408c9e0dd..06dcccfe8e 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -11,99 +11,10 @@ import numpy as np from scipy.interpolate import RegularGridInterpolator +from pvlib import tools import matplotlib.pyplot as plt -def latitude_to_geocentric(phi): - """ - Converts a geodetic (common) latitude to a geocentric latitude. - [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf - """ - a = 6378.137 - b = 6356.752 - return np.arctan(b**2/a**2*np.tan(phi)) - - -def latitude_to_geodetic(phi): - """ - Converts a geocentric latitude to a geodetic (common) latitude. - [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf - """ - a = 6378.137 - b = 6356.752 - return np.arctan(a**2/b**2*np.tan(phi)) - - -def xyz_from_lle(point): - """ - Converts a (lat, lon, elev) tuple into a (x, y, z) tuple. - The center of the earth is the origin in the xyz-space. - The input latitude is assumed to be a common latitude (geodetic). - """ - lat = point[0] - lon = point[1] - elev = point[2] - - a = 6378137.0 - b = 6356752.0 - - # convert to radians - phi = lat*np.pi/180.0 - theta = lon*np.pi/180.0 - - # compute radius of earth at each point - r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 - r = r / (a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2) - r = np.sqrt(r) - - h = r + elev - alpha = latitude_to_geocentric(phi) - beta = theta - x = h * np.cos(alpha) * np.cos(beta) - y = h * np.cos(alpha) * np.sin(beta) - z = h * np.sin(alpha) - v = np.array((x, y, z)) - return v - - -def lle_from_xyz(point): - """ - Converts a (x, y, z) tuple into a (lat, lon, elev) tuple. - The center of the earth is the origin in the xyz-space. - The output latitude is assumed to be a common latitude (geodetic). - """ - a = 6378137.0 - b = 6356752.0 - - x = point[0] - y = point[1] - z = point[2] - - # get corresponding point on earth's surface - t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) - point_s = t * point - z_s = point_s[2] - - elev = np.linalg.norm(point-point_s) - r = np.linalg.norm(point_s) - - alpha = np.arcsin(z_s / r) - phi = latitude_to_geodetic(alpha) - lat = phi*180.0/np.pi - - lon = np.arctan2(y, x)*180/np.pi - return (lat, lon, elev) - - -def pol2cart(rho, phi): - """ - Converts polar coordiantes to cartesian coordinates in 2-d space. - """ - x = rho * np.cos(phi) - y = rho * np.sin(phi) - return (x, y) - - def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): ''' Creates a grid around a location (lat/lon pair) with a specified @@ -208,8 +119,8 @@ def dip_calc(pt1, pt2): phi2 = lat2*np.pi/180.0 theta2 = lon2*np.pi/180.0 - v1 = xyz_from_lle((lat1, lon1, elev1)) - v2 = xyz_from_lle((lat2, lon2, elev2)) + v1 = tools.lle_to_xyz((lat1, lon1, elev1)) + v2 = tools.lle_to_xyz((lat2, lon2, elev2)) x1 = v1[0] y1 = v1[1] @@ -407,7 +318,7 @@ def sample_using_interpolator(grid, num_samples): theta = np.linspace(0, 2 * np.pi, num_samples[1]) polar_pts = np.array(list(itertools.product(r, theta))) - pts = np.array([pol2cart(e[0], e[1]) for e in polar_pts]) + pts = np.array([tools.polar_to_cart(e[0], e[1]) for e in polar_pts]) pts += np.array((site_lat, site_lon)) total_num_samples = num_samples[0]*num_samples[1] @@ -442,9 +353,9 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf """ - c1 = xyz_from_lle(p1) - c2 = xyz_from_lle(p2) - c3 = xyz_from_lle(p3) + c1 = tools.lle_to_xyz(p1) + c2 = tools.lle_to_xyz(p2) + c3 = tools.lle_to_xyz(p3) points = [] for i in range(num_samples): @@ -453,17 +364,10 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): sqrt_r1 = np.sqrt(r1) random_pt = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 - points.append(lle_from_xyz(random_pt)) + points.append(tools.xyz_to_lle(random_pt)) return points -def round_to_nearest(x, base): - """ - Helper function to round x to nearest base. - """ - return base * round(float(x) / base) - - def filter_points(horizon_points, bin_size=1): """ Bins the horizon_points by azimuth values. The azimuth value of each @@ -488,7 +392,7 @@ def filter_points(horizon_points, bin_size=1): for pair in horizon_points: azimuth = pair[0] dip = pair[1] - azimuth_wedge = round_to_nearest(azimuth, bin_size) + azimuth_wedge = tools.round_to_nearest(azimuth, bin_size) if azimuth_wedge in wedges: wedges[azimuth_wedge] = max(dip, wedges[azimuth_wedge]) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 210a1fc4fe..15c5cb17b2 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -1273,7 +1273,7 @@ def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): c = np.cos(tilt_rad) # this gets either an int or an array of zeros - dtf = 0.0 * surface_tilt + dtf = np.multiply(0.0, surface_tilt) for pair in horizon_profile: az = np.radians(pair[0]) horizon_dip = np.radians(pair[1]) diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index 4e87c77fef..cc4d8c79cc 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -729,8 +729,7 @@ def effective_irradiance_model(self): solar_zenith) self.effective_irradiance = self.spectral_modifier * ( - direct*self.aoi_modifier - + fd*self.total_irrad['poa_diffuse']) + direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) return self def complete_irradiance(self, times=None, weather=None): diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 27598e785d..2cf82ec446 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -940,3 +940,51 @@ def test_clearness_index_zenith_independent(airmass_kt): airmass) expected = pd.Series([np.nan, 0.553744437562], index=times) assert_series_equal(out, expected) + + +def test_adjust_direct_for_horizon(irrad_data, ephem_data): + zero_horizon = [] + max_horizon = [] + for i in range(-180, 181): + zero_horizon.append((i, 0.0)) + max_horizon.append((i, 90.0)) + + zero_adjusted = irradiance.adjust_direct_for_horizon(irrad_data["dni"], + zero_horizon, + ephem_data["azimuth"], + ephem_data["zenith"]) + zero_expected = irrad_data["dni"] + assert_allclose(zero_adjusted.values, zero_expected.values) + + max_adjusted = irradiance.adjust_direct_for_horizon(irrad_data["dni"], + max_horizon, + ephem_data["azimuth"], + ephem_data["zenith"]) + + max_expected = pd.Series([0, 0, 0, 0]) + assert_allclose(max_adjusted.values, max_expected.values, atol=1e-7) + + +def test_horizon_adjusted(): + zero_horizon = [] + max_horizon = [] + for i in range(-180, 181): + zero_horizon.append((i, 0.0)) + max_horizon.append((i, 90.0)) + + surface_tilts = [0, 5, 20, 38, 89] + surface_azimuths = [0, 90, 180, 235, 355] + + adjusted = irradiance.horizon_adjusted(surface_tilts, + surface_azimuths, + 1.0, + zero_horizon) + expected = irradiance.isotropic(surface_tilts, np.ones(5)) + assert_allclose(adjusted, expected, atol=2e-3) + + adjusted = irradiance.horizon_adjusted(surface_tilts, + surface_azimuths, + 1.0, + max_horizon) + expected = np.zeros(5) + assert_allclose(adjusted, expected, atol=1e-7) diff --git a/pvlib/tools.py b/pvlib/tools.py index cbd59a8541..be8d52ff21 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -425,3 +425,100 @@ def _golden_sect_DataFrame(params, VL, VH, func): raise Exception("EXCEPTION:iterations exceeded maximum (50)") return func(df, 'V1'), df['V1'] + + +def latitude_to_geocentric(phi): + """ + Converts a geodetic (common) latitude to a geocentric latitude. + [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf + """ + a = 6378.137 + b = 6356.752 + return np.arctan(b**2/a**2*np.tan(phi)) + + +def latitude_to_geodetic(phi): + """ + Converts a geocentric latitude to a geodetic (common) latitude. + [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf + """ + a = 6378.137 + b = 6356.752 + return np.arctan(a**2/b**2*np.tan(phi)) + + +def lle_to_xyz(point): + """ + Converts a (lat, lon, elev) tuple into a (x, y, z) tuple. + The center of the earth is the origin in the xyz-space. + The input latitude is assumed to be a common latitude (geodetic). + """ + lat = point[0] + lon = point[1] + elev = point[2] + + a = 6378137.0 + b = 6356752.0 + + # convert to radians + phi = lat*np.pi/180.0 + theta = lon*np.pi/180.0 + + # compute radius of earth at each point + r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 + r = r / (a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2) + r = np.sqrt(r) + + h = r + elev + alpha = latitude_to_geocentric(phi) + beta = theta + x = h * np.cos(alpha) * np.cos(beta) + y = h * np.cos(alpha) * np.sin(beta) + z = h * np.sin(alpha) + v = np.array((x, y, z)) + return v + + +def xyz_to_lle(point): + """ + Converts a (x, y, z) tuple into a (lat, lon, elev) tuple. + The center of the earth is the origin in the xyz-space. + The output latitude is assumed to be a common latitude (geodetic). + """ + a = 6378137.0 + b = 6356752.0 + + x = point[0] + y = point[1] + z = point[2] + + # get corresponding point on earth's surface + t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) + point_s = t * point + z_s = point_s[2] + + elev = np.linalg.norm(point-point_s) + r = np.linalg.norm(point_s) + + alpha = np.arcsin(z_s / r) + phi = latitude_to_geodetic(alpha) + lat = phi*180.0/np.pi + + lon = np.arctan2(y, x)*180/np.pi + return (lat, lon, elev) + + +def polar_to_cart(rho, phi): + """ + Converts polar coordiantes to cartesian coordinates in 2-d space. + """ + x = rho * np.cos(phi) + y = rho * np.sin(phi) + return (x, y) + + +def round_to_nearest(x, base): + """ + Helper function to round x to nearest base. + """ + return base * round(float(x) / base) diff --git a/setup.py b/setup.py index ac55de8c55..b15cbaf76e 100755 --- a/setup.py +++ b/setup.py @@ -37,7 +37,7 @@ MAINTAINER_EMAIL = 'holmgren@email.arizona.edu' URL = 'https://github.com/pvlib/pvlib-python' -INSTALL_REQUIRES = ['numpy >= 1.10.4', +INSTALL_REQUIRES = ['numpy >= 1.13.3', 'pandas >= 0.18.1', 'pytz', 'requests', From b600e1666cd5546b84f874d56e0b6affe92a1cdc Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 1 Aug 2019 17:00:55 -0700 Subject: [PATCH 10/29] added test_horizon.py and code restructuring --- pvlib/horizon.py | 94 ++++++++++++++++++++++- pvlib/irradiance.py | 83 +------------------- pvlib/test/test_horizon.py | 137 ++++++++++++++++++++++++++++++++++ pvlib/test/test_irradiance.py | 4 +- 4 files changed, 233 insertions(+), 85 deletions(-) create mode 100644 pvlib/test/test_horizon.py diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 06dcccfe8e..88ab9a7057 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -139,6 +139,8 @@ def dip_calc(pt1, pt2): (np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2)*np.cos(theta2-theta1))) bearing_deg = bearing*180.0/np.pi + if bearing_deg < 0: + bearing_deg += 360 return (bearing_deg, dip_angle_deg) @@ -306,8 +308,8 @@ def sample_using_interpolator(grid, num_samples): x_range = x[-1] - x[0] grid_shape = grid.shape - grid_center_i = (grid_shape[0] - 1) / 2 - grid_center_j = (grid_shape[1] - 1) / 2 + grid_center_i = (grid_shape[0] - 1) // 2 + grid_center_j = (grid_shape[1] - 1) // 2 site_lat = grid[grid_center_i, grid_center_j, 0] site_lon = grid[grid_center_i, grid_center_j, 1] @@ -493,3 +495,91 @@ def fake_horizon_profile(max_dip): for i in range(-180, 181): fake_profile.append((i, random.random() * max_dip)) return fake_profile + + +def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): + """ + Determine the dip angle created by the surface of a tilted plane + in a given direction. The dip angle is limited to be non-negative. + + Parameters + ---------- + surface_tilt : numeric + Surface tilt angles in decimal degrees. surface_tilt must be >=0 + and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. surface_azimuth must + be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + direction : numeric + The direction along which the dip angle is to be calculated in + decimal degrees. The convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + Returns + -------- + + dip_angle : numeric + The dip angle in direction created by the tilted plane. + + """ + tilt = np.radians(surface_tilt) + bearing = np.radians(direction - surface_azimuth - 180.0) + + declination = np.degrees(np.arctan(1.0/np.tan(tilt)/np.cos(bearing))) + mask = (declination <= 0) + dip = 90.0 - declination + dip[mask] = 0.0 + + return dip + + # if temp < 0: + # temp = 90.0 + + # x = np.cos(az) + # y = np.sin(az) + # z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) + # + np.sin(plane_az) * np.sin(az)) + # mask = z <= 0 + # numer = x*x + y*y + # denom = x*x + y*y + z*z + # dip = np.degrees(np.arccos(np.sqrt(numer/denom))) + # mask = np.logical_or(np.isnan(dip), z <= 0) + + # if isinstance(dip, pd.Series): + # dip[np.isnan(dip)] = 0 + # dip[mask] = 0 + # else: + # dip = np.where(mask, 0, dip) + # return dip + + +def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): + """ + Calculate the diffuse tilt factor that is adjusted with the horizon + profile. + """ + tilt_rad = np.radians(surface_tilt) + plane_az_rad = np.radians(surface_azimuth) + a = np.sin(tilt_rad) * np.cos(plane_az_rad) + b = np.sin(tilt_rad) * np.sin(plane_az_rad) + c = np.cos(tilt_rad) + + # this gets either an int or an array of zeros + dtf = np.multiply(0.0, surface_tilt) + for pair in horizon_profile: + az = np.radians(pair[0]) + horizon_dip = np.radians(pair[1]) + temp = np.radians(collection_plane_dip_angle(surface_tilt, + surface_azimuth, + pair[0])) + dip = np.maximum(horizon_dip, temp) + + first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * \ + (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) + second_term = .5 * c * np.cos(dip)**2 + dtf += 2 * (first_term + second_term) / len(horizon_profile) + return dtf diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 15c5cb17b2..b0d84ce9d3 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -13,7 +13,7 @@ import numpy as np import pandas as pd -from pvlib import atmosphere, solarposition, tools +from pvlib import atmosphere, solarposition, tools, horizon from pvlib._deprecation import deprecated # see References section of grounddiffuse function @@ -1210,85 +1210,6 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, return sky_diffuse -def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): - """ - Determine the dip angle created by the surface of a tilted plane - in a given direction. The dip angle is limited to be non-negative. - - Parameters - ---------- - surface_tilt : numeric - Surface tilt angles in decimal degrees. surface_tilt must be >=0 - and <=180. The tilt angle is defined as degrees from horizontal - (e.g. surface facing up = 0, surface facing horizon = 90) - - surface_azimuth : numeric - Surface azimuth angles in decimal degrees. surface_azimuth must - be >=0 and <=360. The azimuth convention is defined as degrees - east of north (e.g. North = 0, South=180 East = 90, West = 270). - - direction : numeric - The direction along which the dip angle is to be calculated in - decimal degrees. The convention is defined as degrees - east of north (e.g. North = 0, South=180 East = 90, West = 270). - - Returns - -------- - - dip_angle : numeric - The dip angle in direction created by the tilted plane. - - """ - plane_az = np.radians(surface_azimuth) - tilt = np.radians(surface_tilt) - az = np.radians(direction) - - x = np.cos(az) - y = np.sin(az) - z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) - + np.sin(plane_az) * np.sin(az)) - mask = z <= 0 - numer = x*x + y*y - denom = x*x + y*y + z*z - dip = np.degrees(np.arccos(np.sqrt(numer/denom))) - mask = np.logical_or(np.isnan(dip), z <= 0) - - if isinstance(dip, pd.Series): - dip[np.isnan(dip)] = 0 - dip[mask] = 0 - else: - dip = np.where(mask, 0, dip) - return dip - - -def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): - """ - Calculate the diffuse tilt factor that is adjusted with the horizon - profile. - """ - tilt_rad = np.radians(surface_tilt) - plane_az_rad = np.radians(surface_azimuth) - a = np.sin(tilt_rad) * np.cos(plane_az_rad) - b = np.sin(tilt_rad) * np.sin(plane_az_rad) - c = np.cos(tilt_rad) - - # this gets either an int or an array of zeros - dtf = np.multiply(0.0, surface_tilt) - for pair in horizon_profile: - az = np.radians(pair[0]) - horizon_dip = np.radians(pair[1]) - temp = np.radians(collection_plane_dip_angle(surface_tilt, - surface_azimuth, - pair[0])) - dip = np.maximum(horizon_dip, temp) - - first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * \ - (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) - second_term = .5 * c * np.cos(dip)**2 - dtf += 2 * (first_term + second_term) / len(horizon_profile) - return dtf - - def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): ''' Determine diffuse irradiance from the sky on a tilted surface using @@ -1321,7 +1242,7 @@ def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): The sky diffuse component of the solar radiation on a tilted surface. ''' - dtf = calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) + dtf = horizon.calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) sky_diffuse = dhi * dtf return sky_diffuse diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py new file mode 100644 index 0000000000..9f2a86448c --- /dev/null +++ b/pvlib/test/test_horizon.py @@ -0,0 +1,137 @@ +import numpy as np +from numpy.testing import assert_allclose +from pvlib import horizon, tools + + +def test_grid_lat_lon(): + grid_size = 100 + grid_step = .2 + lat = 24.21 + lon = -35.52 + grid = horizon.grid_lat_lon(lat, lon, + grid_size=grid_size, + grid_step=grid_step) + + assert(grid.shape[0] == 101) + assert(grid[50][50][0] == lat) + assert(grid[49][51][1] == lon + grid_step) + + +def test_dip_calc(): + pt1 = np.array((71.23, -34.70, 1234)) + pt2 = np.array((71.12, -34.16, 124)) + pt3 = np.array((71.29, -35.23, 30044)) + + expected12 = (121.9895, -2.8654) + expected21 = (302.5061, 2.6593) + expected13 = (289.8132, 54.8663) + + actual12 = horizon.dip_calc(pt1, pt2) + actual13 = horizon.dip_calc(pt1, pt3) + actual21 = horizon.dip_calc(pt2, pt1) + assert_allclose(expected12, actual12, rtol=1e-3) + assert_allclose(expected13, actual13, rtol=1e-3) + assert_allclose(expected21, actual21, rtol=1e-3) + + +def test_calculate_horizon_points(): + pass + + +def test_sample_using_grid(): + test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], + [[2, 1, 1], [2, 2, 2], [2, 3, 1]], + [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) + samples = horizon.sample_using_grid(test_grid) + assert(len(samples) == 8) + + +def test_sample_using_triangles(): + test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], + [[2, 1, 1], [2, 2, 2], [2, 3, 1]], + [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) + samples = horizon.sample_using_triangles(test_grid, samples_per_triangle=2) + assert(len(samples) == 32) + + +def test_using_interpolator(): + test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], + [[2, 1, 1], [2, 2, 2], [2, 3, 1]], + [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) + samples = horizon.sample_using_interpolator(test_grid, num_samples=(5, 5)) + assert(len(samples) == 25) + + +def test_uniformly_sample_triangle(): + pt1 = np.array((71.23, -34.70, 1234)) + pt2 = np.array((69.12, -38.16, 124)) + pt3 = np.array((78.23, -36.23, 344)) + points = horizon.uniformly_sample_triangle(pt1, pt2, pt3, 5) + + p1 = tools.lle_to_xyz(pt1) + p2 = tools.lle_to_xyz(pt2) + p3 = tools.lle_to_xyz(pt3) + area = 0.5 * np.linalg.norm(np.cross(p2-p1, p3-p1)) + + for point in points: + p = tools.lle_to_xyz(point) + alpha = 0.5 * np.linalg.norm(np.cross(p2-p, p3-p)) / area + beta = 0.5 * np.linalg.norm(np.cross(p3-p, p1-p)) / area + gamma = 1 - alpha - beta + assert(0 <= alpha <= 1) + assert(0 <= beta <= 1) + assert(0 <= gamma <= 1) + + +def test_filter_points(): + bogus_horizon = [(23, 10), (23.05, 8), (22.56, 14), (55, 2)] + filtered = horizon.filter_points(bogus_horizon, bin_size=1) + assert(len(filtered) == 2) + assert(filtered[0][1] == 14) + + filtered = horizon.filter_points(bogus_horizon, bin_size=.2) + assert(len(filtered) == 3) + assert(filtered[1][1] == 10) + + +def test_collection_plane_dip_angle(): + surface_tilts = np.array([0, 5, 20, 38, 89]) + surface_azimuths = np.array([0, 90, 180, 235, 355]) + directions_easy = np.array([78, 270, 0, 145, 355]) + directions_hard = np.array([729, 220, 60, 115, 3545]) + + expected_easy = np.array([0, 5, 20, 0, 0]) + expected_hard = np.array([0, 3.21873120519, 10.3141048156, + 21.3377447931, 0]) + dips_easy = horizon.collection_plane_dip_angle(surface_tilts, + surface_azimuths, + directions_easy) + assert_allclose(dips_easy, expected_easy) + + dips_hard = horizon.collection_plane_dip_angle(surface_tilts, + surface_azimuths, + directions_hard) + assert_allclose(dips_hard, expected_hard) + + +def calculate_dtf(): + zero_horizon = [] + max_horizon = [] + for i in range(-180, 181): + zero_horizon.append((i, 0.0)) + max_horizon.append((i, 90.0)) + + surface_tilts = np.array([0, 5, 20, 38, 89]) + surface_azimuths = np.array([0, 90, 180, 235, 355]) + + adjusted = horizon.calculate_dtf(zero_horizon, + surface_tilts, + surface_azimuths) + expected = (1 + tools.cosd(surface_tilts)) * 0.5 + assert_allclose(adjusted, expected, atol=2e-3) + + adjusted = horizon.calculate_dtf(max_horizon, + surface_tilts, + surface_azimuths) + expected = np.zeros(5) + assert_allclose(adjusted, expected, atol=1e-7) diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 2cf82ec446..bc207cf6c6 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -972,8 +972,8 @@ def test_horizon_adjusted(): zero_horizon.append((i, 0.0)) max_horizon.append((i, 90.0)) - surface_tilts = [0, 5, 20, 38, 89] - surface_azimuths = [0, 90, 180, 235, 355] + surface_tilts = np.array([0, 5, 20, 38, 89]) + surface_azimuths = np.array([0, 90, 180, 235, 355]) adjusted = irradiance.horizon_adjusted(surface_tilts, surface_azimuths, From b0c9193ad2e637c1283539db45297f56aeee557f Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Fri, 2 Aug 2019 10:53:36 -0700 Subject: [PATCH 11/29] moved horizon adjustmen to isotropic. Fixed some tests --- pvlib/horizon.py | 20 ----------------- pvlib/irradiance.py | 42 +++++++++++++++++++++++++++++++---- pvlib/test/test_horizon.py | 2 +- pvlib/test/test_irradiance.py | 26 ++++++++++++---------- setup.py | 3 ++- 5 files changed, 55 insertions(+), 38 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 88ab9a7057..4361000975 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -536,26 +536,6 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): return dip - # if temp < 0: - # temp = 90.0 - - # x = np.cos(az) - # y = np.sin(az) - # z = -np.tan(tilt) * (np.cos(plane_az) * np.cos(az) - # + np.sin(plane_az) * np.sin(az)) - # mask = z <= 0 - # numer = x*x + y*y - # denom = x*x + y*y + z*z - # dip = np.degrees(np.arccos(np.sqrt(numer/denom))) - # mask = np.logical_or(np.isnan(dip), z <= 0) - - # if isinstance(dip, pd.Series): - # dip[np.isnan(dip)] = 0 - # dip[mask] = 0 - # else: - # dip = np.where(mask, 0, dip) - # return dip - def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): """ diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index b0d84ce9d3..30682ac425 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -429,7 +429,11 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model = model.lower() if model == 'isotropic': - sky = isotropic(surface_tilt, dhi) + if "horizon_profile" in kwargs: + sky = isotropic(surface_tilt, surface_azimuth, dhi, + horizon_profile=kwargs["horizon_profile"]) + else: + sky = isotropic(surface_tilt, surface_azimuth, dhi) elif model == 'klucher': sky = klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith, solar_azimuth) @@ -644,7 +648,7 @@ def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None): get_ground_diffuse) -def isotropic(surface_tilt, dhi): +def isotropic(surface_tilt, surface_azimuth, dhi, horizon_profile=None): r''' Determine diffuse irradiance from the sky on a tilted surface using the isotropic sky model. @@ -657,7 +661,8 @@ def isotropic(surface_tilt, dhi): diffuse irradiance. Thus the diffuse irradiance from the sky (ground reflected irradiance is not included in this algorithm) on a tilted surface can be found from the diffuse horizontal irradiance and the - tilt angle of the surface. + tilt angle of the surface. If a horizon profile is given as an argument + to the function then a numerical integration over the visible part of Parameters ---------- @@ -666,8 +671,17 @@ def isotropic(surface_tilt, dhi): <=180. The tilt angle is defined as degrees from horizontal (e.g. surface facing up = 0, surface facing horizon = 90) + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. surface_azimuth must + be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + dhi : numeric Diffuse horizontal irradiance in W/m^2. DHI must be >=0. + + horizon_profile : list + A list of (azimuth, dip_angle) tuples that defines the horizon + profile Returns ------- @@ -684,8 +698,28 @@ def isotropic(surface_tilt, dhi): heat collector. Trans. ASME 64, 91. ''' - sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5 + ''' + Determine diffuse irradiance from the sky on a tilted surface using + an isotropic model that is adjusted with a horizon profile. + + + Returns + -------- + + sky_diffuse : numeric + The sky diffuse component of the solar radiation on a tilted + surface. + ''' + + if horizon_profile is not None: + dtf = horizon.calculate_dtf(horizon_profile, + surface_tilt, + surface_azimuth) + else: + dtf = (1 + tools.cosd(surface_tilt)) * 0.5 + + sky_diffuse = dhi * dtf return sky_diffuse diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 9f2a86448c..34d5195c27 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -114,7 +114,7 @@ def test_collection_plane_dip_angle(): assert_allclose(dips_hard, expected_hard) -def calculate_dtf(): +def test_calculate_dtf(): zero_horizon = [] max_horizon = [] for i in range(-180, 181): diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index bc207cf6c6..b67c270ad1 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -160,12 +160,12 @@ def test_grounddiffuse_albedo_surface(irrad_data): def test_isotropic_float(): - result = irradiance.isotropic(40, 100) + result = irradiance.isotropic(40, 45, 100) assert_allclose(result, 88.30222215594891) def test_isotropic_series(irrad_data): - result = irradiance.isotropic(40, irrad_data['dhi']) + result = irradiance.isotropic(40, 45, irrad_data['dhi']) assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4) @@ -965,7 +965,7 @@ def test_adjust_direct_for_horizon(irrad_data, ephem_data): assert_allclose(max_adjusted.values, max_expected.values, atol=1e-7) -def test_horizon_adjusted(): +def test_horizon_adjusted_isotropic(): zero_horizon = [] max_horizon = [] for i in range(-180, 181): @@ -975,16 +975,18 @@ def test_horizon_adjusted(): surface_tilts = np.array([0, 5, 20, 38, 89]) surface_azimuths = np.array([0, 90, 180, 235, 355]) - adjusted = irradiance.horizon_adjusted(surface_tilts, - surface_azimuths, - 1.0, - zero_horizon) - expected = irradiance.isotropic(surface_tilts, np.ones(5)) + adjusted = irradiance.isotropic(surface_tilts, + surface_azimuths, + 1.0, + horizon_profile=zero_horizon) + expected = irradiance.isotropic(surface_tilts, + surface_azimuths, + np.ones(5)) assert_allclose(adjusted, expected, atol=2e-3) - adjusted = irradiance.horizon_adjusted(surface_tilts, - surface_azimuths, - 1.0, - max_horizon) + adjusted = irradiance.isotropic(surface_tilts, + surface_azimuths, + 1.0, + horizon_profile=max_horizon) expected = np.zeros(5) assert_allclose(adjusted, expected, atol=1e-7) diff --git a/setup.py b/setup.py index b15cbaf76e..de4125b588 100755 --- a/setup.py +++ b/setup.py @@ -41,7 +41,8 @@ 'pandas >= 0.18.1', 'pytz', 'requests', - 'scipy'] + 'scipy', + 'matplotlib'] TESTS_REQUIRE = ['nose', 'pytest', 'pytest-cov', 'pytest-mock', 'pytest-timeout'] EXTRAS_REQUIRE = { From e6beb326f3dc3d58a7b651e852526878f636c4c4 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Fri, 2 Aug 2019 11:36:10 -0700 Subject: [PATCH 12/29] removed gmaps from horizon. Deleted remnants of horizon adjusment model --- pvlib/horizon.py | 70 --------------------------------------------- pvlib/irradiance.py | 47 ++---------------------------- 2 files changed, 3 insertions(+), 114 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 4361000975..f818066c2f 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -5,7 +5,6 @@ horizon profile. """ -import random import itertools import numpy as np @@ -32,48 +31,6 @@ def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): return grid -def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): - """ - Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). - Queries the googlemaps elevation API to get elevation data at each lat/lon - point. Outputs the original grid with the elevation data appended along - the third axis so the shape is grid_size+1 x grid_size+1 x 3. - """ - - import googlemaps - - in_shape = in_grid.shape - lats = in_grid.T[0].flatten() - longs = in_grid.T[1].flatten() - locations = zip(lats, longs) - gmaps = googlemaps.Client(key=GMAPS_API_KEY) - - out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) - - # Get elevation data from gmaps - elevations = [] - responses = [] - - while len(locations) > 512: - locations_to_request = locations[:512] - locations = locations[512:] - responses += gmaps.elevation(locations=locations_to_request) - responses += gmaps.elevation(locations=locations) - for entry in responses: - elevations.append(entry["elevation"]) - - for i in range(in_shape[0]): - for j in range(in_shape[1]): - lat = in_grid[i, j, 0] - lon = in_grid[i, j, 1] - elevation = elevations[i + j * in_shape[1]] - - out_grid[i, j, 0] = lat - out_grid[i, j, 1] = lon - out_grid[i, j, 2] = elevation - return out_grid - - def dip_calc(pt1, pt2): ''' input: two LLE tuples @@ -470,33 +427,6 @@ def invert_for_pvsyst(horizon_points, hemisphere="north"): return sorted_points -def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): - """ - Uses the functions defined in this modules to generate a complete horizon - profile for a location (specified by lat/lon). An API key for the - googlemaps elevation API is needeed. - """ - - grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) - elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) - horizon_points = calculate_horizon_points(elev_grid, - sampling_method="interpolator", - sampling_param=(1000, 1000)) - filtered_points = filter_points(horizon_points, bin_size=1) - return filtered_points - - -def fake_horizon_profile(max_dip): - """ - Creates a bogus horizon profile by randomly generating dip_angles at - integral azimuth values. Used for testing purposes. - """ - fake_profile = [] - for i in range(-180, 181): - fake_profile.append((i, random.random() * max_dip)) - return fake_profile - - def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): """ Determine the dip angle created by the surface of a tilted plane diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 30682ac425..95d968a200 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -449,10 +449,6 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, sky = perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith, solar_azimuth, airmass, model=model_perez) - elif model == 'horizon_adjusted': - assert("horizon_profile" in kwargs) - sky = horizon_adjusted(surface_tilt, surface_azimuth, dhi, - kwargs["horizon_profile"]) else: raise ValueError('invalid model selection {}'.format(model)) @@ -663,6 +659,7 @@ def isotropic(surface_tilt, surface_azimuth, dhi, horizon_profile=None): surface can be found from the diffuse horizontal irradiance and the tilt angle of the surface. If a horizon profile is given as an argument to the function then a numerical integration over the visible part of + the sky dome is used as the conversion factor. Parameters ---------- @@ -678,8 +675,8 @@ def isotropic(surface_tilt, surface_azimuth, dhi, horizon_profile=None): dhi : numeric Diffuse horizontal irradiance in W/m^2. DHI must be >=0. - - horizon_profile : list + + horizon_profile : list (optional) A list of (azimuth, dip_angle) tuples that defines the horizon profile @@ -1244,44 +1241,6 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, return sky_diffuse -def horizon_adjusted(surface_tilt, surface_azimuth, dhi, horizon_profile): - ''' - Determine diffuse irradiance from the sky on a tilted surface using - an isotropic model that is adjusted with a horizon profile. - - Parameters - ---------- - surface_tilt : numeric - Surface tilt angles in decimal degrees. surface_tilt must be >=0 - and <=180. The tilt angle is defined as degrees from horizontal - (e.g. surface facing up = 0, surface facing horizon = 90) - - surface_azimuth : numeric - Surface azimuth angles in decimal degrees. surface_azimuth must - be >=0 and <=360. The azimuth convention is defined as degrees - east of north (e.g. North = 0, South=180 East = 90, West = 270). - - dhi : numeric - Diffuse horizontal irradiance in W/m^2. DHI must be >=0. - - horizon_profile : list - A list of (azimuth, dip_angle) tuples that defines the horizon - profile - - - Returns - -------- - - sky_diffuse : numeric - The sky diffuse component of the solar radiation on a tilted - surface. - ''' - dtf = horizon.calculate_dtf(horizon_profile, surface_tilt, surface_azimuth) - sky_diffuse = dhi * dtf - - return sky_diffuse - - def clearsky_index(ghi, clearsky_ghi, max_clearsky_index=2.0): """ Calculate the clearsky index. From 10baccdc37e7575b8d60a6ba80d90f24e127c3ae Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 6 Aug 2019 14:09:45 -0700 Subject: [PATCH 13/29] major code restructuring. much more numpy friendly now. still need to update docstrings --- pvlib/horizon.py | 445 +++++++++++++++++----------------- pvlib/irradiance.py | 38 ++- pvlib/test/test_horizon.py | 178 ++++++++++---- pvlib/test/test_irradiance.py | 43 +++- pvlib/tools.py | 38 +-- 5 files changed, 428 insertions(+), 314 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index f818066c2f..c808087dc1 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -4,6 +4,7 @@ as well as a method that uses the googlemaps elevation API to create a horizon profile. """ +from __future__ import division import itertools @@ -11,31 +12,48 @@ from scipy.interpolate import RegularGridInterpolator from pvlib import tools -import matplotlib.pyplot as plt -def grid_lat_lon(lat, lon, grid_size=200, grid_step=.001): +def grid_lat_lon(lat, lon, grid_radius=200, grid_step=.001): ''' - Creates a grid around a location (lat/lon pair) with a specified - grid size and step. The returned grid will be a 3-d matrix with - shape: grid_size+1 x grid_size+1 x 2. + Uses numpy's meshgrid to create grids around a location (lat/lon pair) + with a specified grid radius and step. The grid will be a square with + (2xgrid_radius)+1 points along each side. + + Parameters + ---------- + lat : numeric + The latitude of the location that is to be the center of the grid. + + lon : numeric + The longitude of the location that is to be the center of the grid. + + Returns + ------- + lat_grid: 2d-array + Latitude values at each point on the grid. Values will vary along + axis=1 and will be constant along axis=0. + + lon_grid: 2d-array + Longitude values at each point on the grid. Values will vary along + axis=0 and will be constant along axis=1. ''' - grid = np.ndarray((grid_size + 1, grid_size + 1, 2)) - # fill out grid - for i in range(grid_size + 1): - for j in range(grid_size + 1): - grid[i, j, 0] = lat + (i - grid_size / 2) * grid_step - grid[i, j, 1] = lon + (j - grid_size / 2) * grid_step + lat_start = lat - (grid_radius * grid_step) + lat_stop = lat + (grid_radius * grid_step) + lats = np.linspace(lat_start, lat_stop, 2*grid_radius + 1) + + lon_start = lon - (grid_radius * grid_step) + lon_stop = lon + (grid_radius * grid_step) + lons = np.linspace(lon_start, lon_stop, 2*grid_radius + 1) + + lat_grid, lon_grid = np.meshgrid(lats, lons) - return grid + return lat_grid, lon_grid def dip_calc(pt1, pt2): ''' - input: two LLE tuples - output: distance, dip angle, azimuth - Calculates the dip angle from pt1 to pt2 where dip angle is defined as the angle that the line connecting pt1 to pt2 makes with the plane normal to the Earth's surface at pt2. Also computes the azimuth defined as degrees @@ -43,74 +61,91 @@ def dip_calc(pt1, pt2): Parameters ---------- - pt1 : tuple - (lat, lon, elev) tuple that corresponds to the origin from which - the dip angle is to be calculated. The observer point. + pt1 : ndarray + Nx3 array that contains lat, lon, and elev values that correspond + to the origin points from which the dip angles are to be calculated. + The observer points. - pt1 : tuple - (lat, lon, elev) tuple that corresponds to the origin from which - the dip angle is to be calculated. The observee point. + pt2 : ndarray + Nx3 array that contains lat, lon, and elev values that correspond + to the target points to which the dip angles are to be calculated. + The observee points. Returns ------- - bearing_deg: numeric - The bearing from pt1 to pt2 in degrees East of North + bearing_deg: Nx1 ndarray + The bearings from pt1 to pt2 in degrees East of North. - dip_angle: - The dip angle that pt2 makes with the horizontal as observed at pt1. + dip_angle: Nx1 ndarray + The dip angles that the points in pt2 make with the horizontal + as observed from the points in pt1. ''' a = 6378137.0 b = 6356752.0 - lat1 = pt1[0] - lon1 = pt1[1] - elev1 = pt1[2] - lat2 = pt2[0] - lon2 = pt2[1] - elev2 = pt2[2] + lat1 = np.atleast_1d(pt1.T[0]) + lon1 = np.atleast_1d(pt1.T[1]) + elev1 = np.atleast_1d(pt1.T[2]) + lat2 = np.atleast_1d(pt2.T[0]) + lon2 = np.atleast_1d(pt2.T[1]) + elev2 = np.atleast_1d(pt2.T[2]) # convert to radians - phi1 = lat1*np.pi/180.0 - theta1 = lon1*np.pi/180.0 - phi2 = lat2*np.pi/180.0 - theta2 = lon2*np.pi/180.0 + phi1 = np.radians(lat1) + theta1 = np.radians(lon1) + phi2 = np.radians(lat2) + theta2 = np.radians(lon2) - v1 = tools.lle_to_xyz((lat1, lon1, elev1)) - v2 = tools.lle_to_xyz((lat2, lon2, elev2)) + v1 = tools.lle_to_xyz(np.stack([lat1, lon1, elev1], axis=1)) + v2 = tools.lle_to_xyz(np.stack([lat2, lon2, elev2], axis=1)) + x1 = np.atleast_1d(v1.T[0]) + y1 = np.atleast_1d(v1.T[1]) + z1 = np.atleast_1d(v1.T[2]) - x1 = v1[0] - y1 = v1[1] - z1 = v1[2] + delta = np.atleast_2d(np.subtract(v1, v2)) - delta = np.subtract(v1, v2) + normal = np.atleast_2d(np.stack([2*x1/a**2, 2*y1/a**2, 2*z1/b**2], axis=1)) - normal = np.array((2*x1/a**2, 2*y1/a**2, 2*z1/b**2)) - beta = np.arccos(np.dot(delta, normal)/np.linalg.norm(delta) - / np.linalg.norm(normal)) + # Take the dot product of corresponding vectors + dot = np.sum(np.multiply(delta, normal), axis=1) + beta = np.arccos(dot / np.linalg.norm(delta, axis=1) + / np.linalg.norm(normal, axis=1)) dip_angle = beta - np.pi/2 - dip_angle_deg = dip_angle*180.0/np.pi - # might wanna double check this formula (haversine?) + dip_angle_deg = np.degrees(dip_angle) + bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), (np.cos(phi1) * np.sin(phi2) - np.sin(phi1) * np.cos(phi2)*np.cos(theta2-theta1))) - bearing_deg = bearing*180.0/np.pi - if bearing_deg < 0: - bearing_deg += 360 + bearing_deg = np.degrees(bearing) + + mask = (bearing_deg < 0) + bearing_deg[mask] += 360 - return (bearing_deg, dip_angle_deg) + return bearing_deg, dip_angle_deg -def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): +def calculate_horizon_points(lat_grid, lon_grid, elev_grid, + sampling_method="grid", sampling_param=400): """ - Calculates a horizon profile from a grid of (lat, lon, elev) tuples. - The "site" is assumed to be at the center of the grid. + Calculates a horizon profile from a three grids containing lat, lon, + and elevation values. The "site" is assumed to be at the center of the + grid. Parameters ---------- - grid : ndarray - Assumes + lat_grid : ndarray + A 2d array containing latitude values that correspond to the other + two input grids. + + lon_grid : ndarray + A 2d array containing longitude values that correspond to the other + two input grids. + + elev_grid : ndarray + A 2d array containing elevation values that correspond to the other + two input grids. sampling_method : string A string that specifies the sampling method used to generate the @@ -127,54 +162,52 @@ def calculate_horizon_points(grid, sampling_method="grid", sampling_param=400): point at the center of the grid. """ + assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) - grid_shape = grid.shape - grid_center_i = (grid_shape[0] - 1) / 2 - grid_center_j = (grid_shape[1] - 1) / 2 - site_lat = grid[grid_center_i, grid_center_j, 0] - site_lon = grid[grid_center_i, grid_center_j, 1] - site_elev = grid[grid_center_i, grid_center_j, 2] - site = (site_lat, site_lon, site_elev) - - horizon_points = [] + grid_shape = lat_grid.shape + grid_center_i = (grid_shape[0] - 1) // 2 + grid_center_j = (grid_shape[1] - 1) // 2 + site_lat = lat_grid[grid_center_i, grid_center_j] + site_lon = lon_grid[grid_center_i, grid_center_j] + site_elev = elev_grid[grid_center_i, grid_center_j] + site = np.array([site_lat, site_lon, site_elev]) if sampling_method == "grid": - samples = sample_using_grid(grid) + samples = sample_using_grid(lat_grid, lon_grid, elev_grid) elif sampling_method == "triangles": - samples = sample_using_triangles(grid, sampling_param) + samples = sample_using_triangles(lat_grid, lon_grid, elev_grid, + sampling_param) elif sampling_method == "interpolator": - samples = sample_using_interpolator(grid, sampling_param) + samples = sample_using_interpolator(lat_grid, lon_grid, elev_grid, + sampling_param) - horizon_points = np.array(list(map(lambda pt: dip_calc(site, pt), - samples))) + bearing_deg, dip_angle_deg = dip_calc(site, samples) - return horizon_points + return bearing_deg, dip_angle_deg -def sample_using_grid(grid): +def sample_using_grid(lat_grid, lon_grid, elev_grid): """ Calculates the dip angle from the site (center of the grid) to every point on the grid and uses the results as the horizon profile. """ - grid_shape = grid.shape - grid_center_i = (grid_shape[0] - 1) / 2 - grid_center_j = (grid_shape[1] - 1) / 2 - - all_samples = [] - for i in range(grid_shape[0]): - for j in range(grid_shape[1]): - # make sure the site is excluded - if i != grid_center_i or j != grid_center_j: - lat = grid[i, j, 0] - lon = grid[i, j, 1] - elev = grid[i, j, 2] - all_samples.append((lat, lon, elev)) + assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) + + lats = lat_grid.flatten() + lons = lon_grid.flatten() + elevs = elev_grid.flatten() + samples = np.stack([lats, lons, elevs], axis=1) + assert(samples.shape[1] == 3) + # remove site from samples + all_samples = np.delete(samples, samples.shape[0]//2, axis=0) + assert(all_samples.shape[1] == 3) return all_samples -def sample_using_triangles(grid, samples_per_triangle=10): +def sample_using_triangles(lat_grid, lon_grid, elev_grid, + samples_per_triangle=10): """ Creates triangles using nearest neighbors for every grid point and randomly samples each of these triangles to find dip angles for the horizon profile. @@ -195,46 +228,71 @@ def sample_using_triangles(grid, samples_per_triangle=10): [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf """ + assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) - all_samples = [] - for i in range(grid.shape[0]): - for j in range(grid.shape[1]): - center = (grid[i, j, 0], grid[i, j, 1], grid[i, j, 2]) - if i != 0 and j != 0: - left = (grid[i, j-1, 0], grid[i, j-1, 1], grid[i, j-1, 2]) - top = (grid[i-1, j, 0], grid[i-1, j, 1], grid[i-1, j, 2]) - all_samples += uniformly_sample_triangle(center, - top, - left, - samples_per_triangle) - - if i != 0 and j != grid.shape[1] - 1: - right = (grid[i, j+1, 0], grid[i, j+1, 1], grid[i, j+1, 2]) - top = (grid[i-1, j, 0], grid[i-1, j, 1], grid[i-1, j, 2]) - all_samples += uniformly_sample_triangle(center, - top, - right, - samples_per_triangle) - - if i != grid.shape[0] - 1 and j != 0: - left = (grid[i, j-1, 0], grid[i, j-1, 1], grid[i, j-1, 2]) - bottom = (grid[i+1, j, 0], grid[i+1, j, 1], grid[i+1, j, 2]) - all_samples += uniformly_sample_triangle(center, - bottom, - left, - samples_per_triangle) - - if i != grid.shape[0] - 1 and j != grid.shape[1] - 1: - right = (grid[i, j+1, 0], grid[i, j+1, 1], grid[i, j+1, 2]) - bottom = (grid[i+1, j, 0], grid[i+1, j, 1], grid[i+1, j, 2]) - all_samples += uniformly_sample_triangle(center, - bottom, - right, - samples_per_triangle) - return all_samples + # start with empty array + all_samples = np.array([], dtype=np.float64).reshape(0, 3) - -def sample_using_interpolator(grid, num_samples): + for i in range(lat_grid.shape[0]): + for j in range(lat_grid.shape[1]): + center = np.array([lat_grid[i, j], + lon_grid[i, j], + elev_grid[i, j]]) + if i != 0 and j != 0: + left = np.array([lat_grid[i, j-1], + lon_grid[i, j-1], + elev_grid[i, j-1]]) + top = np.array([lat_grid[i-1, j], + lon_grid[i-1, j], + elev_grid[i-1, j]]) + samples = uniformly_sample_triangle(center, + top, + left, + samples_per_triangle) + all_samples = np.vstack([all_samples, samples]) + + if i != 0 and j != lat_grid.shape[1] - 1: + right = np.array([lat_grid[i, j+1], + lon_grid[i, j+1], + elev_grid[i, j+1]]) + top = np.array([lat_grid[i-1, j], + lon_grid[i-1, j], + elev_grid[i-1, j]]) + samples = uniformly_sample_triangle(center, + top, + right, + samples_per_triangle) + all_samples = np.vstack([all_samples, samples]) + + if i != lat_grid.shape[0] - 1 and j != 0: + left = np.array([lat_grid[i, j-1], + lon_grid[i, j-1], + elev_grid[i, j-1]]) + bottom = np.array([lat_grid[i+1, j], + lon_grid[i+1, j], + elev_grid[i+1, j]]) + samples = uniformly_sample_triangle(center, + bottom, + left, + samples_per_triangle) + all_samples = np.vstack([all_samples, samples]) + + if i != lat_grid.shape[0] - 1 and j != lat_grid.shape[1] - 1: + right = np.array([lat_grid[i, j+1], + lon_grid[i, j+1], + elev_grid[i, j+1]]) + bottom = np.array([lat_grid[i+1, j], + lon_grid[i+1, j], + elev_grid[i+1, j]]) + samples = uniformly_sample_triangle(center, + bottom, + right, + samples_per_triangle) + all_samples = np.vstack([all_samples, samples]) + return np.array(all_samples) + + +def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): """ Creates a "grid" using polar coordinates and uses the scipy's grid interpolator to estimate elevation values at each point on the polar grid @@ -259,21 +317,22 @@ def sample_using_interpolator(grid, num_samples): List of (azimuth, dip_angle) tuples taken from the polar grid """ - x = grid.T[0][0] - y = grid.T[1].T[0] + assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) - x_range = x[-1] - x[0] + lats = lat_grid[0] + lons = lon_grid.T[0] - grid_shape = grid.shape + lat_range = lats[-1] - lats[0] + + grid_shape = lat_grid.shape grid_center_i = (grid_shape[0] - 1) // 2 grid_center_j = (grid_shape[1] - 1) // 2 - site_lat = grid[grid_center_i, grid_center_j, 0] - site_lon = grid[grid_center_i, grid_center_j, 1] + site_lat = lat_grid[grid_center_i, grid_center_j] + site_lon = lon_grid[grid_center_i, grid_center_j] - elevs = grid.T[2].T - interpolator = RegularGridInterpolator((x, y), elevs) + interpolator = RegularGridInterpolator((lats, lons), elev_grid.T) - r = np.linspace(0, x_range/2, num_samples[0]) + r = np.linspace(0, lat_range//2, num_samples[0]) theta = np.linspace(0, 2 * np.pi, num_samples[1]) polar_pts = np.array(list(itertools.product(r, theta))) @@ -282,7 +341,7 @@ def sample_using_interpolator(grid, num_samples): total_num_samples = num_samples[0]*num_samples[1] interpolated_elevs = interpolator(pts).reshape(total_num_samples, 1) - samples = np.concatenate((pts, interpolated_elevs), axis=1) + samples = np.hstack((pts, interpolated_elevs)) return samples @@ -316,18 +375,17 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): c2 = tools.lle_to_xyz(p2) c3 = tools.lle_to_xyz(p3) - points = [] - for i in range(num_samples): - r1 = np.random.rand() - r2 = np.random.rand() - sqrt_r1 = np.sqrt(r1) + r1 = np.random.rand(num_samples).reshape((num_samples, 1)) + r2 = np.random.rand(num_samples).reshape((num_samples, 1)) + sqrt_r1 = np.sqrt(r1) + + random_pts = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 + random_pts = tools.xyz_to_lle(random_pts) - random_pt = (1-sqrt_r1)*c1 + sqrt_r1*(1-r2)*c2 + sqrt_r1*r2*c3 - points.append(tools.xyz_to_lle(random_pt)) - return points + return random_pts -def filter_points(horizon_points, bin_size=1): +def filter_points(horizon_azimuths, horizon_angles, bin_size=1): """ Bins the horizon_points by azimuth values. The azimuth value of each point in horizon_points is rounded to the nearest bin and then the @@ -343,14 +401,14 @@ def filter_points(horizon_points, bin_size=1): Returns ------- - sorted_points: list - List of (azimuth, dip_angle) values that correspond to the greatest - dip_angle in each azimuth bin. + """ + assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) + wedges = {} - for pair in horizon_points: - azimuth = pair[0] - dip = pair[1] + for i in range(horizon_angles.shape[0]): + azimuth = horizon_azimuths[i] + dip = horizon_angles[i] azimuth_wedge = tools.round_to_nearest(azimuth, bin_size) if azimuth_wedge in wedges: @@ -358,73 +416,15 @@ def filter_points(horizon_points, bin_size=1): else: wedges[azimuth_wedge] = dip - filtered_points = [] - for key in wedges.keys(): - filtered_points.append((key, wedges[key])) - - sorted_points = sorted(filtered_points, key=lambda tup: tup[0]) - return sorted_points - - -def visualize(horizon_profile, pvsyst_scale=False): - """ - Plots a horizon profile with azimuth on the x-axis and dip angle on the y. - """ - azimuths = [] - dips = [] - for pair in horizon_profile: - azimuth = pair[0] - azimuths.append(azimuth) - dips.append(pair[1]) - plt.figure(figsize=(10, 6)) - if pvsyst_scale: - plt.ylim(0, 90) - plt.plot(azimuths, dips, "-") - plt.show - - -def polar_plot(horizon_profile): - """ - Plots a horizon profile on a polar plot with dip angle as the raidus and - azimuth as the theta value. An offset of 5 is added to the dip_angle to - make the plot more readable with low dip angles. - """ - azimuths = [] - dips = [] - for pair in horizon_profile: - azimuth = pair[0] - azimuths.append(np.radians(azimuth)) - dips.append(pair[1] + 5) - plt.figure(figsize=(10, 6)) - sp = plt.subplot(1, 1, 1, projection='polar') - sp.set_theta_zero_location('N') - sp.set_theta_direction(-1) - plt.plot(azimuths, dips, "o") - plt.show - - -def invert_for_pvsyst(horizon_points, hemisphere="north"): - """ - Modify the azimuth values in horizon_points to match PVSyst's azimuth - convention (which is dependent on hemisphere) - """ - - # look at that northern hemisphere bias right there - # not even sorry. - assert hemisphere == "north" or hemisphere == "south" + filtered_azimuths = [] + filtered_angles = [] + for key in sorted(wedges.keys()): + filtered_azimuths.append(key) + filtered_angles.append(wedges[key]) - inverted_points = [] - for pair in horizon_points: - azimuth = pair[0] - if hemisphere == "north": - azimuth -= 180 - if azimuth < -180: - azimuth += 360 - elif hemisphere == "south": - azimuth = -azimuth - inverted_points.append((azimuth, pair[1])) - sorted_points = sorted(inverted_points, key=lambda tup: tup[0]) - return sorted_points + filtered_angles = np.array(filtered_angles) + filtered_azimuths = np.array(filtered_azimuths) + return filtered_azimuths, filtered_angles def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): @@ -467,29 +467,32 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): return dip -def calculate_dtf(horizon_profile, surface_tilt, surface_azimuth): +def calculate_dtf(horizon_azimuths, horizon_angles, + surface_tilt, surface_azimuth): """ Calculate the diffuse tilt factor that is adjusted with the horizon profile. """ + assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) tilt_rad = np.radians(surface_tilt) plane_az_rad = np.radians(surface_azimuth) a = np.sin(tilt_rad) * np.cos(plane_az_rad) b = np.sin(tilt_rad) * np.sin(plane_az_rad) c = np.cos(tilt_rad) - # this gets either an int or an array of zeros + # this gets either a float or an array of zeros dtf = np.multiply(0.0, surface_tilt) - for pair in horizon_profile: - az = np.radians(pair[0]) - horizon_dip = np.radians(pair[1]) + num_points = horizon_azimuths.shape[0] + for i in range(horizon_azimuths.shape[0]): + az = np.radians(horizon_azimuths[i]) + horizon_dip = np.radians(horizon_angles[i]) temp = np.radians(collection_plane_dip_angle(surface_tilt, surface_azimuth, - pair[0])) + horizon_azimuths[i])) dip = np.maximum(horizon_dip, temp) first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * \ (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) second_term = .5 * c * np.cos(dip)**2 - dtf += 2 * (first_term + second_term) / len(horizon_profile) + dtf += 2 * (first_term + second_term) / num_points return dtf diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 95d968a200..c58fe04b1e 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -429,9 +429,10 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model = model.lower() if model == 'isotropic': - if "horizon_profile" in kwargs: + if "horizon_azimuths" in kwargs and "horizon_angles" in kwargs: sky = isotropic(surface_tilt, surface_azimuth, dhi, - horizon_profile=kwargs["horizon_profile"]) + horizon_azimuths=kwargs["horizon_azimuths"], + horizon_angles=kwargs["horizon_angles"]) else: sky = isotropic(surface_tilt, surface_azimuth, dhi) elif model == 'klucher': @@ -644,7 +645,8 @@ def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None): get_ground_diffuse) -def isotropic(surface_tilt, surface_azimuth, dhi, horizon_profile=None): +def isotropic(surface_tilt, surface_azimuth, dhi, + horizon_azimuths=None, horizon_angles=None): r''' Determine diffuse irradiance from the sky on a tilted surface using the isotropic sky model. @@ -709,8 +711,9 @@ def isotropic(surface_tilt, surface_azimuth, dhi, horizon_profile=None): surface. ''' - if horizon_profile is not None: - dtf = horizon.calculate_dtf(horizon_profile, + if horizon_azimuths is not None and horizon_angles is not None: + dtf = horizon.calculate_dtf(horizon_azimuths, + horizon_angles, surface_tilt, surface_azimuth) else: @@ -2923,7 +2926,7 @@ def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1, return dni -def adjust_direct_for_horizon(poa_direct, horizon_profile, +def adjust_direct_for_horizon(poa_direct, horizon_angles, solar_azimuth, solar_zenith): ''' Adjusts modeled DNI to account for the presence of a horizon. At @@ -2934,30 +2937,23 @@ def adjust_direct_for_horizon(poa_direct, horizon_profile, ---------- poa_direct : numeric The modeled direct normal irradiance in the plane of array. + horizon_angles : numeric + A list or vector of 361 values where the ith element corresponds + to the dip angle of the horizon at i degrees of azimuth. solar_zenith : numeric Solar zenith angle. solar_azimuth : numeric Solar azimuth angle. - horizon_profile : list - A list of (azimuth, dip_angle) tuples that defines the horizon - profile given. Every azimuth from -180 to 180, in 1 degree - increments, must be in this list. Returns ------- adjusted_irradiance : Series POA direct normal irradiance after adjusting for the horizon. ''' - dip_angles = {} - for pair in horizon_profile: - dip_angles[pair[0]] = pair[1] - adjusted_irradiance = poa_direct - for i in range(len(poa_direct)): - rounded_solar_azimuth = round(solar_azimuth[i]) - if rounded_solar_azimuth > 180: - rounded_solar_azimuth -= 360 - horizon_dip = dip_angles[rounded_solar_azimuth] - if (solar_zenith[i] > (90 - horizon_dip)): - adjusted_irradiance[i] = 0 + adjusted_irradiance = poa_direct + rounded_solar_azimuth = np.round(solar_azimuth).astype(int) + horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] + mask = solar_zenith > horizon_zenith + adjusted_irradiance[mask] = 0 return adjusted_irradiance diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 34d5195c27..9049f00927 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -4,17 +4,19 @@ def test_grid_lat_lon(): - grid_size = 100 + grid_radius = 50 grid_step = .2 lat = 24.21 lon = -35.52 - grid = horizon.grid_lat_lon(lat, lon, - grid_size=grid_size, - grid_step=grid_step) + lat_grid, lon_grid = horizon.grid_lat_lon(lat, lon, + grid_radius=grid_radius, + grid_step=grid_step) - assert(grid.shape[0] == 101) - assert(grid[50][50][0] == lat) - assert(grid[49][51][1] == lon + grid_step) + assert(lat_grid.shape[0] == 101) + assert(lat_grid.shape == lon_grid.shape) + assert_allclose(lat_grid[50][50], lat) + assert_allclose(lon_grid[49][51], lon - grid_step) + assert_allclose(lat_grid[30][72], lat + 22*grid_step) def test_dip_calc(): @@ -22,43 +24,114 @@ def test_dip_calc(): pt2 = np.array((71.12, -34.16, 124)) pt3 = np.array((71.29, -35.23, 30044)) - expected12 = (121.9895, -2.8654) - expected21 = (302.5061, 2.6593) - expected13 = (289.8132, 54.8663) + test_pts = np.vstack([pt1, pt2]) + reverse_test_pts = np.vstack([pt2, pt1]) + + expected_bearings = np.array([121.9895, 302.5061]) + expected_dips = np.array([-2.8654, 2.6593]) + + expected13 = np.array([[289.8132], [54.8663]]) + + act_bearings, act_dips = horizon.dip_calc(test_pts, reverse_test_pts) + assert_allclose(act_bearings, expected_bearings, rtol=1e-3) + assert_allclose(act_dips, expected_dips, rtol=1e-3) - actual12 = horizon.dip_calc(pt1, pt2) actual13 = horizon.dip_calc(pt1, pt3) - actual21 = horizon.dip_calc(pt2, pt1) - assert_allclose(expected12, actual12, rtol=1e-3) assert_allclose(expected13, actual13, rtol=1e-3) - assert_allclose(expected21, actual21, rtol=1e-3) def test_calculate_horizon_points(): - pass + test_lat_grid = np.array([[1, 2, 3], + [1, 2, 3], + [1, 2, 3]]) + + test_lon_grid = np.array([[-3, -3, -3], + [-2, -2, -2], + [-1, -1, -1]]) + + test_elev_grid = np.array([[15, 18, 43], + [212, 135, 1], + [36, 145, 5]]) + + bearings, dips = horizon.calculate_horizon_points(test_lat_grid, + test_lon_grid, + test_elev_grid, + sampling_method="grid") + + expected_bearings = np.array([0, 90, 45, 270, 135]) + rounded_bearings = np.round(bearings).astype(int) + assert(bearings.shape == dips.shape) + assert(np.all(np.in1d(expected_bearings, rounded_bearings))) + + bears, dips = horizon.calculate_horizon_points(test_lat_grid, + test_lon_grid, + test_elev_grid, + sampling_method="triangles", + sampling_param=5) + assert(bears.shape == dips.shape) + + bears, _ = horizon.calculate_horizon_points(test_lat_grid, + test_lon_grid, + test_elev_grid, + sampling_method="interpolator", + sampling_param=(10, 10)) + assert(bears.shape[0] == 100) def test_sample_using_grid(): - test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], - [[2, 1, 1], [2, 2, 2], [2, 3, 1]], - [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) - samples = horizon.sample_using_grid(test_grid) + test_lat_grid = np.array([[1, 1, 3], + [2, 1, 1], + [3, 1, 5]]) + + test_lon_grid = np.array([[1, 1, -3], + [2, -1, 1], + [3, -15, 5]]) + + test_elev_grid = np.array([[15, 18, 43], + [212, 135, 1], + [36, 145, 5]]) + + samples = horizon.sample_using_grid(test_lat_grid, + test_lon_grid, + test_elev_grid) assert(len(samples) == 8) def test_sample_using_triangles(): - test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], - [[2, 1, 1], [2, 2, 2], [2, 3, 1]], - [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) - samples = horizon.sample_using_triangles(test_grid, samples_per_triangle=2) + test_lat_grid = np.array([[1, 1, 3], + [2, 1, 1], + [3, 1, 5]]) + + test_lon_grid = np.array([[1, 1, -3], + [2, -1, 1], + [3, -15, 5]]) + + test_elev_grid = np.array([[15, 18, 43], + [212, 135, 1], + [36, 145, 5]]) + samples = horizon.sample_using_triangles(test_lat_grid, + test_lon_grid, + test_elev_grid, + samples_per_triangle=2) assert(len(samples) == 32) def test_using_interpolator(): - test_grid = np.array([[[1, 1, 3], [1, 2, 8], [1, 3, 8]], - [[2, 1, 1], [2, 2, 2], [2, 3, 1]], - [[3, 1, 5], [3, 2, 7], [3, 3, 9]]]) - samples = horizon.sample_using_interpolator(test_grid, num_samples=(5, 5)) + test_lat_grid = np.array([[1, 2, 3], + [1, 2, 3], + [1, 2, 3]]) + + test_lon_grid = np.array([[-3, -3, -3], + [-2, -2, -2], + [-1, -1, -1]]) + + test_elev_grid = np.array([[15, 18, 43], + [212, 135, 1], + [36, 145, 5]]) + samples = horizon.sample_using_interpolator(test_lat_grid, + test_lon_grid, + test_elev_grid, + num_samples=(5, 5)) assert(len(samples) == 25) @@ -74,6 +147,7 @@ def test_uniformly_sample_triangle(): area = 0.5 * np.linalg.norm(np.cross(p2-p1, p3-p1)) for point in points: + print(point) p = tools.lle_to_xyz(point) alpha = 0.5 * np.linalg.norm(np.cross(p2-p, p3-p)) / area beta = 0.5 * np.linalg.norm(np.cross(p3-p, p1-p)) / area @@ -84,14 +158,19 @@ def test_uniformly_sample_triangle(): def test_filter_points(): - bogus_horizon = [(23, 10), (23.05, 8), (22.56, 14), (55, 2)] - filtered = horizon.filter_points(bogus_horizon, bin_size=1) - assert(len(filtered) == 2) - assert(filtered[0][1] == 14) + test_azimuths = np.array([23, 23.05, 22.56, 55]) + bogus_horizon = np.array([10, 8, 14, 2]) + filtered_azimuths, filtered_angles = horizon.filter_points(test_azimuths, + bogus_horizon, + bin_size=1) + assert(filtered_azimuths.shape[0] == filtered_angles.shape[0]) + assert(filtered_angles[0] == 14) - filtered = horizon.filter_points(bogus_horizon, bin_size=.2) - assert(len(filtered) == 3) - assert(filtered[1][1] == 10) + filtered_azimuths, filtered_angles = horizon.filter_points(test_azimuths, + bogus_horizon, + bin_size=.2) + assert(filtered_azimuths.shape[0] == 3) + assert(filtered_angles[1] == 10) def test_collection_plane_dip_angle(): @@ -115,23 +194,40 @@ def test_collection_plane_dip_angle(): def test_calculate_dtf(): - zero_horizon = [] - max_horizon = [] - for i in range(-180, 181): - zero_horizon.append((i, 0.0)) - max_horizon.append((i, 90.0)) + num_points = 360 + test_azimuths = np.arange(0, num_points, dtype=np.float64) + zero_horizon = np.zeros(num_points) + max_horizon = np.full((num_points), 90.0) + uniform_horizon = np.full((num_points), 7.0) + random_horizon = np.random.random(num_points) * 7 surface_tilts = np.array([0, 5, 20, 38, 89]) surface_azimuths = np.array([0, 90, 180, 235, 355]) - adjusted = horizon.calculate_dtf(zero_horizon, + adjusted = horizon.calculate_dtf(test_azimuths, + zero_horizon, surface_tilts, surface_azimuths) expected = (1 + tools.cosd(surface_tilts)) * 0.5 assert_allclose(adjusted, expected, atol=2e-3) - adjusted = horizon.calculate_dtf(max_horizon, + adjusted = horizon.calculate_dtf(test_azimuths, + max_horizon, surface_tilts, surface_azimuths) expected = np.zeros(5) assert_allclose(adjusted, expected, atol=1e-7) + + adjusted = horizon.calculate_dtf(test_azimuths, + random_horizon, + surface_tilts, + surface_azimuths) + min_random_dtf = horizon.calculate_dtf(test_azimuths, + uniform_horizon, + surface_tilts, + surface_azimuths) + max_random_dtf = (1 + tools.cosd(surface_tilts)) * 0.5 + + mask = np.logical_and((adjusted >= min_random_dtf), + (adjusted <= max_random_dtf)) + assert(np.all(mask)) diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index b67c270ad1..5d53cccd61 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -943,11 +943,8 @@ def test_clearness_index_zenith_independent(airmass_kt): def test_adjust_direct_for_horizon(irrad_data, ephem_data): - zero_horizon = [] - max_horizon = [] - for i in range(-180, 181): - zero_horizon.append((i, 0.0)) - max_horizon.append((i, 90.0)) + zero_horizon = np.zeros(361) + max_horizon = np.full((361), 90.0) zero_adjusted = irradiance.adjust_direct_for_horizon(irrad_data["dni"], zero_horizon, @@ -966,11 +963,12 @@ def test_adjust_direct_for_horizon(irrad_data, ephem_data): def test_horizon_adjusted_isotropic(): - zero_horizon = [] - max_horizon = [] - for i in range(-180, 181): - zero_horizon.append((i, 0.0)) - max_horizon.append((i, 90.0)) + num_points = 360 + test_azimuths = np.arange(0, num_points, dtype=np.float64) + zero_horizon = np.zeros(num_points) + max_horizon = np.full((num_points), 90.0) + uniform_horizon = np.full((num_points), 7.0) + random_horizon = np.random.random(num_points) * 7 surface_tilts = np.array([0, 5, 20, 38, 89]) surface_azimuths = np.array([0, 90, 180, 235, 355]) @@ -978,15 +976,34 @@ def test_horizon_adjusted_isotropic(): adjusted = irradiance.isotropic(surface_tilts, surface_azimuths, 1.0, - horizon_profile=zero_horizon) + horizon_azimuths=test_azimuths, + horizon_angles=zero_horizon) expected = irradiance.isotropic(surface_tilts, surface_azimuths, - np.ones(5)) + 1.0) assert_allclose(adjusted, expected, atol=2e-3) adjusted = irradiance.isotropic(surface_tilts, surface_azimuths, 1.0, - horizon_profile=max_horizon) + horizon_azimuths=test_azimuths, + horizon_angles=max_horizon) expected = np.zeros(5) assert_allclose(adjusted, expected, atol=1e-7) + + adjusted = irradiance.isotropic(surface_tilts, + surface_azimuths, + 1.0, + horizon_azimuths=test_azimuths, + horizon_angles=random_horizon) + min_random_dtf = irradiance.isotropic(surface_tilts, + surface_azimuths, + 1.0, + horizon_azimuths=test_azimuths, + horizon_angles=uniform_horizon) + max_random_dtf = irradiance.isotropic(surface_tilts, + surface_azimuths, + 1.0) + mask = np.logical_and((adjusted >= min_random_dtf), + (adjusted <= max_random_dtf)) + assert(np.all(mask)) diff --git a/pvlib/tools.py b/pvlib/tools.py index be8d52ff21..fcd13c3238 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -453,16 +453,16 @@ def lle_to_xyz(point): The center of the earth is the origin in the xyz-space. The input latitude is assumed to be a common latitude (geodetic). """ - lat = point[0] - lon = point[1] - elev = point[2] + lat = np.atleast_1d(point.T[0]) + lon = np.atleast_1d(point.T[1]) + elev = np.atleast_1d(point.T[2]) a = 6378137.0 b = 6356752.0 # convert to radians - phi = lat*np.pi/180.0 - theta = lon*np.pi/180.0 + phi = np.radians(lat) + theta = np.radians(lon) # compute radius of earth at each point r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 @@ -475,7 +475,7 @@ def lle_to_xyz(point): x = h * np.cos(alpha) * np.cos(beta) y = h * np.cos(alpha) * np.sin(beta) z = h * np.sin(alpha) - v = np.array((x, y, z)) + v = np.stack([x, y, z], axis=1) return v @@ -488,24 +488,26 @@ def xyz_to_lle(point): a = 6378137.0 b = 6356752.0 - x = point[0] - y = point[1] - z = point[2] + x = np.atleast_1d(point.T[0]) + y = np.atleast_1d(point.T[1]) + z = np.atleast_1d(point.T[2]) # get corresponding point on earth's surface + t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) + t = t.reshape(point.shape[0], 1) point_s = t * point - z_s = point_s[2] + z_s = point_s.T[2] + + elev = np.linalg.norm(point-point_s, axis=1) - elev = np.linalg.norm(point-point_s) - r = np.linalg.norm(point_s) + r = np.linalg.norm(point_s, axis=1) alpha = np.arcsin(z_s / r) phi = latitude_to_geodetic(alpha) - lat = phi*180.0/np.pi - - lon = np.arctan2(y, x)*180/np.pi - return (lat, lon, elev) + lat = np.degrees(phi) + lon = np.degrees(np.arctan2(y, x)) + return np.stack([lat, lon, elev], axis=1) def polar_to_cart(rho, phi): @@ -514,11 +516,11 @@ def polar_to_cart(rho, phi): """ x = rho * np.cos(phi) y = rho * np.sin(phi) - return (x, y) + return(x, y) def round_to_nearest(x, base): """ Helper function to round x to nearest base. """ - return base * round(float(x) / base) + return base * np.round(x / base) From 345eaa062a31de054d9aa1d47ee1e62436d5c9af Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 6 Aug 2019 14:44:50 -0700 Subject: [PATCH 14/29] updated docstrings --- pvlib/horizon.py | 141 +++++++++++++++++++++++++++++++++++--------- pvlib/irradiance.py | 15 +++-- 2 files changed, 125 insertions(+), 31 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index c808087dc1..846ee2f44f 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -64,12 +64,12 @@ def dip_calc(pt1, pt2): pt1 : ndarray Nx3 array that contains lat, lon, and elev values that correspond to the origin points from which the dip angles are to be calculated. - The observer points. + The observer points. (will also take a 1darray with 3 elements) pt2 : ndarray Nx3 array that contains lat, lon, and elev values that correspond to the target points to which the dip angles are to be calculated. - The observee points. + The observee points. (will also take a 1darray with 3 elements) Returns ------- @@ -157,9 +157,13 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, Returns ------- - horizon_points: list - List of (azimuth, dip_angle) tuples that define the horizon at the - point at the center of the grid. + bearing_deg: Nx1 ndarray + The bearings from the "site" to sampled points in degrees + East of North. + + dip_angle_dip: Nx1 ndarray + The dip angles that the sampled points make with the horizontal + as observed from the the "site". """ assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) @@ -189,9 +193,27 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, def sample_using_grid(lat_grid, lon_grid, elev_grid): """ Calculates the dip angle from the site (center of the grid) - to every point on the grid and uses the results as the - horizon profile. + to every point on the grid. + + Parameters + ---------- + lat_grid : ndarray + A 2d array containing latitude values that correspond to the other + two input grids. + + lon_grid : ndarray + A 2d array containing longitude values that correspond to the other + two input grids. + + elev_grid : ndarray + A 2d array containing elevation values that correspond to the other + two input grids. + + Returns + ------- + all_samples: Nx3 ndarray + Array of lat, lon, elev points that are grid points. """ assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) @@ -199,10 +221,9 @@ def sample_using_grid(lat_grid, lon_grid, elev_grid): lons = lon_grid.flatten() elevs = elev_grid.flatten() samples = np.stack([lats, lons, elevs], axis=1) - assert(samples.shape[1] == 3) # remove site from samples + all_samples = np.delete(samples, samples.shape[0]//2, axis=0) - assert(all_samples.shape[1] == 3) return all_samples @@ -214,8 +235,17 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, Parameters ---------- - grid : ndarray - Grid that contains lat lon and elev information. + lat_grid : ndarray + A 2d array containing latitude values that correspond to the other + two input grids. + + lon_grid : ndarray + A 2d array containing longitude values that correspond to the other + two input grids. + + elev_grid : ndarray + A 2d array containing elevation values that correspond to the other + two input grids. samples_per_triangle : numeric The number of random samples to be uniformly taken from the surface @@ -223,8 +253,8 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, Returns ------- - all_samples: list - List of (azimuth, dip_angle) tuples that were sampled from the grid + all_samples: Nx3 ndarray + Array of [lat, lon, elev] points that were sampled from the grid. [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf """ @@ -302,8 +332,17 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): Parameters ---------- - grid : ndarray - Grid that contains lat lon and elev information. + lat_grid : ndarray + A 2d array containing latitude values that correspond to the other + two input grids. + + lon_grid : ndarray + A 2d array containing longitude values that correspond to the other + two input grids. + + elev_grid : ndarray + A 2d array containing elevation values that correspond to the other + two input grids. num_samples : tuple A tuple containing two integers. The first is the desired number of @@ -314,7 +353,7 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): Returns ------- all_samples: list - List of (azimuth, dip_angle) tuples taken from the polar grid + Array of [lat, lon, elev] points that were sampled using the polar grid. """ assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) @@ -352,12 +391,15 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): Parameters ---------- - pt1 : tuple - A (lat, lon, elev) tuple that defines one vertex of the triangle. - pt2 : tuple - A (lat, lon, elev) tuple that defines another vertex of the triangle. - pt3 : tuple - A (lat, lon, elev) tuple that defines the last vertex of the triangle. + pt1 : ndarray + An array conaining (lat, lon, elev) values that define one vertex + of the triangle. + pt2 : ndarray + An array conaining (lat, lon, elev) values that define another vertex + of the triangle. + pt3 : ndarray + An array conaining (lat, lon, elev) values that define the last vertex + of the triangle. num_samples : tuple The number of random samples to be uniformly taken from the surface @@ -365,8 +407,8 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): Returns ------- - points: list - List of (lat, lon, elev) tuples that lie on the surface of the + points: Nx3 ndarray + Array with N (lat, lon, elev) points that lie on the surface of the triangle. [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf @@ -393,14 +435,29 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): Parameters ---------- - horizon_points : list - List of (azimuth, dip_angle) tuples that define the horizon. + horizon_azimuths: numeric + Azimuth values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in horizon_angles. + + horizon_angle: numeric + Dip angle values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in + horizon_azimuths. bin_size : int - The bin size of azimuth values. + The width of the bins for the azimuth values. Returns ------- + filtered_azimuths: numeric + Azimuth values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in + filtered_angles. + + filtered_angles: numeric + Dip angle values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in + filtered_azimuths. """ assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) @@ -472,6 +529,36 @@ def calculate_dtf(horizon_azimuths, horizon_angles, """ Calculate the diffuse tilt factor that is adjusted with the horizon profile. + + Parameters + ---------- + horizon_azimuths: Nx1 numeric + Azimuth values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in horizon_angles. + + horizon_angles: Nx1 numeric + Dip angle values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in + horizon_azimuths. + + surface_tilt : numeric + Surface tilt angles in decimal degrees. surface_tilt must be >=0 + and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. surface_azimuth must + be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + Returns + ------- + dtf: numeric + The diffuse tilt factor that can be multiplied with the diffuse + horizontal irradiance (DHI) to get the incident irradiance from + the sky that is adjusted for the horizon profile and the tilt of + the plane. + """ assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) tilt_rad = np.radians(surface_tilt) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index c58fe04b1e..1a4d897437 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -678,9 +678,14 @@ def isotropic(surface_tilt, surface_azimuth, dhi, dhi : numeric Diffuse horizontal irradiance in W/m^2. DHI must be >=0. - horizon_profile : list (optional) - A list of (azimuth, dip_angle) tuples that defines the horizon - profile + horizon_azimuths: numeric (optional) + Azimuth values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in horizon_angles. + + horizon_angle: numeric (optional) + Dip angle values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in + horizon_azimuths. Returns ------- @@ -2939,7 +2944,9 @@ def adjust_direct_for_horizon(poa_direct, horizon_angles, The modeled direct normal irradiance in the plane of array. horizon_angles : numeric A list or vector of 361 values where the ith element corresponds - to the dip angle of the horizon at i degrees of azimuth. + to the dip angle of the horizon at i degrees of azimuth. Note that + 0 and 360 degrees of azimuth correspond to the same direction + but both are required. solar_zenith : numeric Solar zenith angle. solar_azimuth : numeric From 5405d7544821b861693443fad221048becf52a3b Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 10:30:00 -0700 Subject: [PATCH 15/29] docstring changes --- pvlib/horizon.py | 18 ++++++++++++------ pvlib/irradiance.py | 11 +++++++---- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 846ee2f44f..82a9606c68 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -440,8 +440,11 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): element in this array corresponds to the ith element in horizon_angles. horizon_angle: numeric - Dip angle values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in + Dip angle values for points that define the horizon profile. The dip + angle of the horizon is the angle that the horizon makes with the + horizontal. It is given in degrees. If the horizon appears above + the horizontal, then the dip angle is positive. The ith element in + this array corresponds to the ith element in horizon_azimuths. bin_size : int @@ -532,13 +535,16 @@ def calculate_dtf(horizon_azimuths, horizon_angles, Parameters ---------- - horizon_azimuths: Nx1 numeric + horizon_azimuths: numeric Azimuth values for points that define the horizon profile. The ith element in this array corresponds to the ith element in horizon_angles. - horizon_angles: Nx1 numeric - Dip angle values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in + horizon_angles: numeric + Dip angle values for points that define the horizon profile. The dip + angle of the horizon is the angle that the horizon makes with the + horizontal. It is given in degrees. If the horizon appears above + the horizontal, then the dip angle is positive. The ith element in + this array corresponds to the ith element in horizon_azimuths. surface_tilt : numeric diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 1a4d897437..40867a6e43 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -13,7 +13,7 @@ import numpy as np import pandas as pd -from pvlib import atmosphere, solarposition, tools, horizon +from pvlib import atmosphere, horizon, solarposition, tools from pvlib._deprecation import deprecated # see References section of grounddiffuse function @@ -682,9 +682,12 @@ def isotropic(surface_tilt, surface_azimuth, dhi, Azimuth values for points that define the horizon profile. The ith element in this array corresponds to the ith element in horizon_angles. - horizon_angle: numeric (optional) - Dip angle values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in + horizon_angles: numeric (optional) + Dip angle values for points that define the horizon profile. The dip + angle of the horizon is the angle that the horizon makes with the + horizontal. It is given in degrees. If the horizon appears above + the horizontal, then the dip angle is positive. The ith element in + this array corresponds to the ith element in horizon_azimuths. Returns From cc481c715e89a1f132dbaaf826b987a44d88030c Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 11:22:20 -0700 Subject: [PATCH 16/29] made some changes to modelchain and location due to restructuring of horizon profile --- pvlib/irradiance.py | 4 +++- pvlib/location.py | 18 ++++++++++++++++-- pvlib/modelchain.py | 27 ++++++++++++++++++--------- pvlib/test/test_modelchain.py | 19 +++++++++++++++++++ 4 files changed, 56 insertions(+), 12 deletions(-) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 40867a6e43..111e46d063 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -661,7 +661,9 @@ def isotropic(surface_tilt, surface_azimuth, dhi, surface can be found from the diffuse horizontal irradiance and the tilt angle of the surface. If a horizon profile is given as an argument to the function then a numerical integration over the visible part of - the sky dome is used as the conversion factor. + the sky dome is used as the conversion factor. Note that the horizon + correction is designed for satellite irradiance data, not irradiance + data from ground stations. Parameters ---------- diff --git a/pvlib/location.py b/pvlib/location.py index 57cd77d727..7ea3cac1eb 100644 --- a/pvlib/location.py +++ b/pvlib/location.py @@ -48,6 +48,18 @@ class Location(object): name : None or string, default None. Sets the name attribute of the Location object. + horizon_azimuths: numeric + Azimuth values for points that define the horizon profile. The ith + element in this array corresponds to the ith element in horizon_angles. + + horizon_angles: numeric + Dip angle values for points that define the horizon profile. The dip + angle of the horizon is the angle that the horizon makes with the + horizontal. It is given in degrees. If the horizon appears above + the horizontal, then the dip angle is positive. The ith element in + this array corresponds to the ith element in + horizon_azimuths. + **kwargs Arbitrary keyword arguments. Included for compatibility, but not used. @@ -58,7 +70,8 @@ class Location(object): """ def __init__(self, latitude, longitude, tz='UTC', altitude=0, - name=None, horizon_profile=None, **kwargs): + name=None, horizon_azimuths=None, horizon_angles=None, + **kwargs): self.latitude = latitude self.longitude = longitude @@ -76,7 +89,8 @@ def __init__(self, latitude, longitude, tz='UTC', altitude=0, raise TypeError('Invalid tz specification') self.altitude = altitude - self.horizon_profile = horizon_profile + self.horizon_angles = horizon_angles + self.horizon_azimuths = horizon_azimuths self.name = name def __repr__(self): diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index cc4d8c79cc..63ca5c3d1b 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -719,14 +719,22 @@ def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) direct = self.total_irrad['poa_direct'] - if self.location.horizon_profile is not None: - horizon = self.location.horizon_profile - solar_zenith = self.solar_position['apparent_zenith'] - solar_az = self.solar_position['azimuth'] - direct = pvlib.irradiance.adjust_direct_for_horizon(direct, - horizon, - solar_az, - solar_zenith) + if self.location.horizon_angles is not None: + horizon = self.location.horizon_angles + if len(horizon) != 361: + warnings.warn('The location used in modelchain has specified' + 'horizon_angles but does not contain 361 values.' + 'The horizon is not being used to correct DNI.' + 'To correct DNI with the horizon, horizon_angles' + 'must contain 361 values.', + pvlibDeprecationWarning) + else: + zenith = self.solar_position['apparent_zenith'] + solar_az = self.solar_position['azimuth'] + direct = pvlib.irradiance.adjust_direct_for_horizon(direct, + horizon, + solar_az, + zenith) self.effective_irradiance = self.spectral_modifier * ( direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) @@ -900,7 +908,8 @@ def prepare_inputs(self, times=None, weather=None): self.weather['dhi'], airmass=self.airmass['airmass_relative'], model=self.transposition_model, - horizon_profile=self.location.horizon_profile) + horizon_azimuths=self.location.horizon_azimuths, + horizon_angles=self.location.horizon_angles) if self.weather.get('wind_speed') is None: self.weather['wind_speed'] = 0 diff --git a/pvlib/test/test_modelchain.py b/pvlib/test/test_modelchain.py index 7fac62542c..1a70a2fecc 100644 --- a/pvlib/test/test_modelchain.py +++ b/pvlib/test/test_modelchain.py @@ -117,6 +117,13 @@ def pvwatts_dc_pvwatts_ac_system(sam_data): def location(): return Location(32.2, -111, altitude=700) +@pytest.fixture +def location_with_horizon(): + horizon_azimuths = np.arange(0, 361, dtype=np.float64) + horizon_angles = np.full((361), 10.0) + return Location(32.2, -111, altitude=700, + horizon_azimuths=horizon_azimuths, + horizon_angles=horizon_angles) @pytest.fixture def weather(): @@ -169,6 +176,18 @@ def test_run_model_with_irradiance(system, location): assert_series_equal(ac, expected) +def test_run_model_with_horizon(system, location_with_horizon): + mc = ModelChain(system, location_with_horizon) + times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') + irradiance = pd.DataFrame({'dni': 900, 'ghi': 600, 'dhi': 150}, + index=times) + ac = mc.run_model(times, weather=irradiance).ac + + expected = pd.Series(np.array([ 1.90054749e+02, -2.00000000e-02]), + index=times) + assert_series_equal(ac, expected) + + def test_run_model_perez(system, location): mc = ModelChain(system, location, transposition_model='perez') times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') From 15e59ebeec1640679570426f143a7ab4a13b6f41 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 11:37:18 -0700 Subject: [PATCH 17/29] added one more test to modelchain --- pvlib/test/test_modelchain.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/pvlib/test/test_modelchain.py b/pvlib/test/test_modelchain.py index 1a70a2fecc..ad68d8d6aa 100644 --- a/pvlib/test/test_modelchain.py +++ b/pvlib/test/test_modelchain.py @@ -176,14 +176,27 @@ def test_run_model_with_irradiance(system, location): assert_series_equal(ac, expected) -def test_run_model_with_horizon(system, location_with_horizon): +def test_run_model_with_full_horizon(system, location_with_horizon): mc = ModelChain(system, location_with_horizon) times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') irradiance = pd.DataFrame({'dni': 900, 'ghi': 600, 'dhi': 150}, index=times) ac = mc.run_model(times, weather=irradiance).ac - expected = pd.Series(np.array([ 1.90054749e+02, -2.00000000e-02]), + expected = pd.Series(np.array([1.90054749e+02, -2.00000000e-02]), + index=times) + assert_series_equal(ac, expected) + + +def test_run_model_with_partial_horizon(system, location): + location.horizon_angles = np.array([0, 0, 0]) + mc = ModelChain(system, location) + times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') + irradiance = pd.DataFrame({'dni': 900, 'ghi': 600, 'dhi': 150}, + index=times) + ac = mc.run_model(times, weather=irradiance).ac + + expected = pd.Series(np.array([1.90054749e+02, -2.00000000e-02]), index=times) assert_series_equal(ac, expected) From 33c0bb8b9012e7aff3c12e241ac7943692084d3a Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 12:44:54 -0700 Subject: [PATCH 18/29] threw code and some docs into horizon.rst --- docs/sphinx/source/horizon.rst | 154 +++++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 docs/sphinx/source/horizon.rst diff --git a/docs/sphinx/source/horizon.rst b/docs/sphinx/source/horizon.rst new file mode 100644 index 0000000000..c042434fbd --- /dev/null +++ b/docs/sphinx/source/horizon.rst @@ -0,0 +1,154 @@ +.. _horizon: + +Horizon +======== + +.. ipython:: python + :suppress: + + import pandas as pd + from pvlib import pvsystem + + +The :py:mod:`~pvlib.horizon` module contains many functions for horizon +profile modeling. + +The horizon profile at a location is a mapping from azimuth to elevation angle +(also called dip angle). + +A source of elevation data is needed for many of the +A common source of elevation data is NASA's SRTM [1]. An easily accessible source +of elevation data (albeit not free) is the googlemaps elevation API. There +are a few examples of how to query the googlemaps elevation API further below. + + +def fake_horizon_profile(max_dip): + """ + Creates a bogus horizon profile by randomly generating dip_angles at + integral azimuth values. Used for testing purposes. + """ + fake_profile = [] + for i in range(-180, 181): + fake_profile.append((i, random.random() * max_dip)) + return fake_profile + + + +def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): + """ + Uses the functions defined in this modules to generate a complete horizon + profile for a location (specified by lat/lon). An API key for the + googlemaps elevation API is needeed. + """ + + grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) + elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) + horizon_points = calculate_horizon_points(elev_grid, + sampling_method="interpolator", + sampling_param=(1000, 1000)) + filtered_points = filter_points(horizon_points, bin_size=1) + return filtered_points + + + +def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): + """ + Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). + Queries the googlemaps elevation API to get elevation data at each lat/lon + point. Outputs the original grid with the elevation data appended along + the third axis so the shape is grid_size+1 x grid_size+1 x 3. + """ + + import googlemaps + + in_shape = in_grid.shape + lats = in_grid.T[0].flatten() + longs = in_grid.T[1].flatten() + locations = zip(lats, longs) + gmaps = googlemaps.Client(key=GMAPS_API_KEY) + + out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) + + # Get elevation data from gmaps + elevations = [] + responses = [] + + while len(locations) > 512: + locations_to_request = locations[:512] + locations = locations[512:] + responses += gmaps.elevation(locations=locations_to_request) + responses += gmaps.elevation(locations=locations) + for entry in responses: + elevations.append(entry["elevation"]) + + for i in range(in_shape[0]): + for j in range(in_shape[1]): + lat = in_grid[i, j, 0] + lon = in_grid[i, j, 1] + elevation = elevations[i + j * in_shape[1]] + + out_grid[i, j, 0] = lat + out_grid[i, j, 1] = lon + out_grid[i, j, 2] = elevation + return out_grid + + + +def visualize(horizon_profile, pvsyst_scale=False): + """ + Plots a horizon profile with azimuth on the x-axis and dip angle on the y. + """ + azimuths = [] + dips = [] + for pair in horizon_profile: + azimuth = pair[0] + azimuths.append(azimuth) + dips.append(pair[1]) + plt.figure(figsize=(10, 6)) + if pvsyst_scale: + plt.ylim(0, 90) + plt.plot(azimuths, dips, "-") + plt.show + + +def polar_plot(horizon_profile): + """ + Plots a horizon profile on a polar plot with dip angle as the raidus and + azimuth as the theta value. An offset of 5 is added to the dip_angle to + make the plot more readable with low dip angles. + """ + azimuths = [] + dips = [] + for pair in horizon_profile: + azimuth = pair[0] + azimuths.append(np.radians(azimuth)) + dips.append(pair[1] + 5) + plt.figure(figsize=(10, 6)) + sp = plt.subplot(1, 1, 1, projection='polar') + sp.set_theta_zero_location('N') + sp.set_theta_direction(-1) + plt.plot(azimuths, dips, "o") + plt.show + +def invert_for_pvsyst(horizon_points, hemisphere="north"): + """ + Modify the azimuth values in horizon_points to match PVSyst's azimuth + convention (which is dependent on hemisphere) + """ + + # look at that northern hemisphere bias right there + # not even sorry. + assert hemisphere == "north" or hemisphere == "south" + + inverted_points = [] + for pair in horizon_points: + azimuth = pair[0] + if hemisphere == "north": + azimuth -= 180 + if azimuth < -180: + azimuth += 360 + elif hemisphere == "south": + azimuth = -azimuth + inverted_points.append((azimuth, pair[1])) + sorted_points = sorted(inverted_points, key=lambda tup: tup[0]) + return sorted_points \ No newline at end of file From 18ccd1a34e16f69bbd1fd0beacdee7be39d5c77c Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 14:07:21 -0700 Subject: [PATCH 19/29] reverted irradiance, location and modelchain (and tests) to master --- pvlib/irradiance.py | 98 +++-------------------------------- pvlib/location.py | 18 +------ pvlib/modelchain.py | 26 ++-------- pvlib/test/test_irradiance.py | 71 +------------------------ pvlib/test/test_modelchain.py | 32 ------------ 5 files changed, 14 insertions(+), 231 deletions(-) diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 111e46d063..000cf3216c 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -13,7 +13,7 @@ import numpy as np import pandas as pd -from pvlib import atmosphere, horizon, solarposition, tools +from pvlib import atmosphere, solarposition, tools from pvlib._deprecation import deprecated # see References section of grounddiffuse function @@ -365,7 +365,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth, poa_sky_diffuse = get_sky_diffuse( surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=dni_extra, airmass=airmass, model=model, - model_perez=model_perez, **kwargs) + model_perez=model_perez) poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo, surface_type) @@ -383,8 +383,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=None, airmass=None, model='isotropic', - model_perez='allsitescomposite1990', - **kwargs): + model_perez='allsitescomposite1990'): r""" Determine in-plane sky diffuse irradiance component using the specified sky diffuse irradiance model. @@ -429,12 +428,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model = model.lower() if model == 'isotropic': - if "horizon_azimuths" in kwargs and "horizon_angles" in kwargs: - sky = isotropic(surface_tilt, surface_azimuth, dhi, - horizon_azimuths=kwargs["horizon_azimuths"], - horizon_angles=kwargs["horizon_angles"]) - else: - sky = isotropic(surface_tilt, surface_azimuth, dhi) + sky = isotropic(surface_tilt, dhi) elif model == 'klucher': sky = klucher(surface_tilt, surface_azimuth, dhi, ghi, solar_zenith, solar_azimuth) @@ -645,8 +639,7 @@ def get_ground_diffuse(surface_tilt, ghi, albedo=.25, surface_type=None): get_ground_diffuse) -def isotropic(surface_tilt, surface_azimuth, dhi, - horizon_azimuths=None, horizon_angles=None): +def isotropic(surface_tilt, dhi): r''' Determine diffuse irradiance from the sky on a tilted surface using the isotropic sky model. @@ -659,11 +652,7 @@ def isotropic(surface_tilt, surface_azimuth, dhi, diffuse irradiance. Thus the diffuse irradiance from the sky (ground reflected irradiance is not included in this algorithm) on a tilted surface can be found from the diffuse horizontal irradiance and the - tilt angle of the surface. If a horizon profile is given as an argument - to the function then a numerical integration over the visible part of - the sky dome is used as the conversion factor. Note that the horizon - correction is designed for satellite irradiance data, not irradiance - data from ground stations. + tilt angle of the surface. Parameters ---------- @@ -672,26 +661,9 @@ def isotropic(surface_tilt, surface_azimuth, dhi, <=180. The tilt angle is defined as degrees from horizontal (e.g. surface facing up = 0, surface facing horizon = 90) - surface_azimuth : numeric - Surface azimuth angles in decimal degrees. surface_azimuth must - be >=0 and <=360. The azimuth convention is defined as degrees - east of north (e.g. North = 0, South=180 East = 90, West = 270). - dhi : numeric Diffuse horizontal irradiance in W/m^2. DHI must be >=0. - horizon_azimuths: numeric (optional) - Azimuth values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in horizon_angles. - - horizon_angles: numeric (optional) - Dip angle values for points that define the horizon profile. The dip - angle of the horizon is the angle that the horizon makes with the - horizontal. It is given in degrees. If the horizon appears above - the horizontal, then the dip angle is positive. The ith element in - this array corresponds to the ith element in - horizon_azimuths. - Returns ------- diffuse : numeric @@ -707,29 +679,8 @@ def isotropic(surface_tilt, surface_azimuth, dhi, heat collector. Trans. ASME 64, 91. ''' - ''' - Determine diffuse irradiance from the sky on a tilted surface using - an isotropic model that is adjusted with a horizon profile. + sky_diffuse = dhi * (1 + tools.cosd(surface_tilt)) * 0.5 - - - Returns - -------- - - sky_diffuse : numeric - The sky diffuse component of the solar radiation on a tilted - surface. - ''' - - if horizon_azimuths is not None and horizon_angles is not None: - dtf = horizon.calculate_dtf(horizon_azimuths, - horizon_angles, - surface_tilt, - surface_azimuth) - else: - dtf = (1 + tools.cosd(surface_tilt)) * 0.5 - - sky_diffuse = dhi * dtf return sky_diffuse @@ -2934,38 +2885,3 @@ def dni(ghi, dhi, zenith, clearsky_dni=None, clearsky_tolerance=1.1, (zenith < zenith_threshold_for_zero_dni) & (dni > max_dni)] = max_dni return dni - - -def adjust_direct_for_horizon(poa_direct, horizon_angles, - solar_azimuth, solar_zenith): - ''' - Adjusts modeled DNI to account for the presence of a horizon. At - each time step check to see if the sun is behind the horizon. If so, - set the DNI to zero. - - Parameters - ---------- - poa_direct : numeric - The modeled direct normal irradiance in the plane of array. - horizon_angles : numeric - A list or vector of 361 values where the ith element corresponds - to the dip angle of the horizon at i degrees of azimuth. Note that - 0 and 360 degrees of azimuth correspond to the same direction - but both are required. - solar_zenith : numeric - Solar zenith angle. - solar_azimuth : numeric - Solar azimuth angle. - - Returns - ------- - adjusted_irradiance : Series - POA direct normal irradiance after adjusting for the horizon. - ''' - - adjusted_irradiance = poa_direct - rounded_solar_azimuth = np.round(solar_azimuth).astype(int) - horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] - mask = solar_zenith > horizon_zenith - adjusted_irradiance[mask] = 0 - return adjusted_irradiance diff --git a/pvlib/location.py b/pvlib/location.py index 7ea3cac1eb..68b2062907 100644 --- a/pvlib/location.py +++ b/pvlib/location.py @@ -48,18 +48,6 @@ class Location(object): name : None or string, default None. Sets the name attribute of the Location object. - horizon_azimuths: numeric - Azimuth values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in horizon_angles. - - horizon_angles: numeric - Dip angle values for points that define the horizon profile. The dip - angle of the horizon is the angle that the horizon makes with the - horizontal. It is given in degrees. If the horizon appears above - the horizontal, then the dip angle is positive. The ith element in - this array corresponds to the ith element in - horizon_azimuths. - **kwargs Arbitrary keyword arguments. Included for compatibility, but not used. @@ -70,8 +58,7 @@ class Location(object): """ def __init__(self, latitude, longitude, tz='UTC', altitude=0, - name=None, horizon_azimuths=None, horizon_angles=None, - **kwargs): + name=None, **kwargs): self.latitude = latitude self.longitude = longitude @@ -89,8 +76,7 @@ def __init__(self, latitude, longitude, tz='UTC', altitude=0, raise TypeError('Invalid tz specification') self.altitude = altitude - self.horizon_angles = horizon_angles - self.horizon_azimuths = horizon_azimuths + self.name = name def __repr__(self): diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py index 63ca5c3d1b..1455a57ad3 100644 --- a/pvlib/modelchain.py +++ b/pvlib/modelchain.py @@ -717,27 +717,9 @@ def no_extra_losses(self): def effective_irradiance_model(self): fd = self.system.module_parameters.get('FD', 1.) - direct = self.total_irrad['poa_direct'] - - if self.location.horizon_angles is not None: - horizon = self.location.horizon_angles - if len(horizon) != 361: - warnings.warn('The location used in modelchain has specified' - 'horizon_angles but does not contain 361 values.' - 'The horizon is not being used to correct DNI.' - 'To correct DNI with the horizon, horizon_angles' - 'must contain 361 values.', - pvlibDeprecationWarning) - else: - zenith = self.solar_position['apparent_zenith'] - solar_az = self.solar_position['azimuth'] - direct = pvlib.irradiance.adjust_direct_for_horizon(direct, - horizon, - solar_az, - zenith) - self.effective_irradiance = self.spectral_modifier * ( - direct*self.aoi_modifier + fd*self.total_irrad['poa_diffuse']) + self.total_irrad['poa_direct']*self.aoi_modifier + + fd*self.total_irrad['poa_diffuse']) return self def complete_irradiance(self, times=None, weather=None): @@ -907,9 +889,7 @@ def prepare_inputs(self, times=None, weather=None): self.weather['ghi'], self.weather['dhi'], airmass=self.airmass['airmass_relative'], - model=self.transposition_model, - horizon_azimuths=self.location.horizon_azimuths, - horizon_angles=self.location.horizon_angles) + model=self.transposition_model) if self.weather.get('wind_speed') is None: self.weather['wind_speed'] = 0 diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 5d53cccd61..27598e785d 100644 --- a/pvlib/test/test_irradiance.py +++ b/pvlib/test/test_irradiance.py @@ -160,12 +160,12 @@ def test_grounddiffuse_albedo_surface(irrad_data): def test_isotropic_float(): - result = irradiance.isotropic(40, 45, 100) + result = irradiance.isotropic(40, 100) assert_allclose(result, 88.30222215594891) def test_isotropic_series(irrad_data): - result = irradiance.isotropic(40, 45, irrad_data['dhi']) + result = irradiance.isotropic(40, irrad_data['dhi']) assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4) @@ -940,70 +940,3 @@ def test_clearness_index_zenith_independent(airmass_kt): airmass) expected = pd.Series([np.nan, 0.553744437562], index=times) assert_series_equal(out, expected) - - -def test_adjust_direct_for_horizon(irrad_data, ephem_data): - zero_horizon = np.zeros(361) - max_horizon = np.full((361), 90.0) - - zero_adjusted = irradiance.adjust_direct_for_horizon(irrad_data["dni"], - zero_horizon, - ephem_data["azimuth"], - ephem_data["zenith"]) - zero_expected = irrad_data["dni"] - assert_allclose(zero_adjusted.values, zero_expected.values) - - max_adjusted = irradiance.adjust_direct_for_horizon(irrad_data["dni"], - max_horizon, - ephem_data["azimuth"], - ephem_data["zenith"]) - - max_expected = pd.Series([0, 0, 0, 0]) - assert_allclose(max_adjusted.values, max_expected.values, atol=1e-7) - - -def test_horizon_adjusted_isotropic(): - num_points = 360 - test_azimuths = np.arange(0, num_points, dtype=np.float64) - zero_horizon = np.zeros(num_points) - max_horizon = np.full((num_points), 90.0) - uniform_horizon = np.full((num_points), 7.0) - random_horizon = np.random.random(num_points) * 7 - - surface_tilts = np.array([0, 5, 20, 38, 89]) - surface_azimuths = np.array([0, 90, 180, 235, 355]) - - adjusted = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0, - horizon_azimuths=test_azimuths, - horizon_angles=zero_horizon) - expected = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0) - assert_allclose(adjusted, expected, atol=2e-3) - - adjusted = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0, - horizon_azimuths=test_azimuths, - horizon_angles=max_horizon) - expected = np.zeros(5) - assert_allclose(adjusted, expected, atol=1e-7) - - adjusted = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0, - horizon_azimuths=test_azimuths, - horizon_angles=random_horizon) - min_random_dtf = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0, - horizon_azimuths=test_azimuths, - horizon_angles=uniform_horizon) - max_random_dtf = irradiance.isotropic(surface_tilts, - surface_azimuths, - 1.0) - mask = np.logical_and((adjusted >= min_random_dtf), - (adjusted <= max_random_dtf)) - assert(np.all(mask)) diff --git a/pvlib/test/test_modelchain.py b/pvlib/test/test_modelchain.py index ad68d8d6aa..7fac62542c 100644 --- a/pvlib/test/test_modelchain.py +++ b/pvlib/test/test_modelchain.py @@ -117,13 +117,6 @@ def pvwatts_dc_pvwatts_ac_system(sam_data): def location(): return Location(32.2, -111, altitude=700) -@pytest.fixture -def location_with_horizon(): - horizon_azimuths = np.arange(0, 361, dtype=np.float64) - horizon_angles = np.full((361), 10.0) - return Location(32.2, -111, altitude=700, - horizon_azimuths=horizon_azimuths, - horizon_angles=horizon_angles) @pytest.fixture def weather(): @@ -176,31 +169,6 @@ def test_run_model_with_irradiance(system, location): assert_series_equal(ac, expected) -def test_run_model_with_full_horizon(system, location_with_horizon): - mc = ModelChain(system, location_with_horizon) - times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') - irradiance = pd.DataFrame({'dni': 900, 'ghi': 600, 'dhi': 150}, - index=times) - ac = mc.run_model(times, weather=irradiance).ac - - expected = pd.Series(np.array([1.90054749e+02, -2.00000000e-02]), - index=times) - assert_series_equal(ac, expected) - - -def test_run_model_with_partial_horizon(system, location): - location.horizon_angles = np.array([0, 0, 0]) - mc = ModelChain(system, location) - times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') - irradiance = pd.DataFrame({'dni': 900, 'ghi': 600, 'dhi': 150}, - index=times) - ac = mc.run_model(times, weather=irradiance).ac - - expected = pd.Series(np.array([1.90054749e+02, -2.00000000e-02]), - index=times) - assert_series_equal(ac, expected) - - def test_run_model_perez(system, location): mc = ModelChain(system, location, transposition_model='perez') times = pd.date_range('20160101 1200-0700', periods=2, freq='6H') From d1119abbcfd073a981abd9179269af9614ebb08f Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 14:31:41 -0700 Subject: [PATCH 20/29] changed dip angles to elevation angles --- pvlib/horizon.py | 118 +++++++++++++++++++------------------ pvlib/test/test_horizon.py | 63 ++++++++++---------- 2 files changed, 93 insertions(+), 88 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 82a9606c68..c50ccd9dc2 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -52,33 +52,36 @@ def grid_lat_lon(lat, lon, grid_radius=200, grid_step=.001): return lat_grid, lon_grid -def dip_calc(pt1, pt2): +def elevation_angle_calc(pt1, pt2): ''' - Calculates the dip angle from pt1 to pt2 where dip angle is defined as - the angle that the line connecting pt1 to pt2 makes with the plane normal - to the Earth's surface at pt2. Also computes the azimuth defined as degrees - East of North the bearing of pt2 from pt1. This uses the Haversine formula. + Calculates the elevation angle from pt1 to pt2 where elevation angle is + defined as the angle between the line connecting pt1 to pt2 and the plane + normal to the Earth's surface at pt1. A point that appears above the + horizontal has a positive elevation angle. Also computes the azimuth + defined as degrees East of North the bearing of pt2 from pt1. + This uses the Haversine formula. Parameters ---------- pt1 : ndarray Nx3 array that contains lat, lon, and elev values that correspond - to the origin points from which the dip angles are to be calculated. - The observer points. (will also take a 1darray with 3 elements) + to the origin points from which the elevation angles are to be + calculated. The observer points. pt2 : ndarray Nx3 array that contains lat, lon, and elev values that correspond - to the target points to which the dip angles are to be calculated. - The observee points. (will also take a 1darray with 3 elements) + to the target points to which the elevation angles are to be + calculated. The observee points. Returns ------- - bearing_deg: Nx1 ndarray + bearing_deg: numeric The bearings from pt1 to pt2 in degrees East of North. - dip_angle: Nx1 ndarray - The dip angles that the points in pt2 make with the horizontal - as observed from the points in pt1. + elevation_angle_deg: numeric + The elevation angles that the points in pt2 make with the horizontal + as observed from the points in pt1. Given in degrees above the + horizontal. ''' a = 6378137.0 @@ -111,9 +114,9 @@ def dip_calc(pt1, pt2): dot = np.sum(np.multiply(delta, normal), axis=1) beta = np.arccos(dot / np.linalg.norm(delta, axis=1) / np.linalg.norm(normal, axis=1)) - dip_angle = beta - np.pi/2 + elevation_angle = beta - np.pi/2 - dip_angle_deg = np.degrees(dip_angle) + elevation_angle_deg = np.degrees(elevation_angle) bearing = np.arctan2(np.sin(theta2-theta1)*np.cos(phi2), (np.cos(phi1) * np.sin(phi2) @@ -123,7 +126,7 @@ def dip_calc(pt1, pt2): mask = (bearing_deg < 0) bearing_deg[mask] += 360 - return bearing_deg, dip_angle_deg + return bearing_deg, elevation_angle_deg def calculate_horizon_points(lat_grid, lon_grid, elev_grid, @@ -161,9 +164,10 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, The bearings from the "site" to sampled points in degrees East of North. - dip_angle_dip: Nx1 ndarray - The dip angles that the sampled points make with the horizontal - as observed from the the "site". + elevation_angle_deg: numeric + The angles that the sampled points make with the horizontal + as observed from the "site". Given in degrees above the + horizontal. """ assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) @@ -185,14 +189,14 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, samples = sample_using_interpolator(lat_grid, lon_grid, elev_grid, sampling_param) - bearing_deg, dip_angle_deg = dip_calc(site, samples) + bearing_deg, elevation_angle_deg = elevation_angle_calc(site, samples) - return bearing_deg, dip_angle_deg + return bearing_deg, elevation_angle_deg def sample_using_grid(lat_grid, lon_grid, elev_grid): """ - Calculates the dip angle from the site (center of the grid) + Calculates the Elevation angle from the site (center of the grid) to every point on the grid. Parameters @@ -231,7 +235,7 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, samples_per_triangle=10): """ Creates triangles using nearest neighbors for every grid point and randomly - samples each of these triangles to find dip angles for the horizon profile. + samples each of these triangles to find Elevation angles for the horizon profile. Parameters ---------- @@ -326,7 +330,7 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): """ Creates a "grid" using polar coordinates and uses the scipy's grid interpolator to estimate elevation values at each point on the polar grid - from the input (rectangular) grid that has true elevation values. Dip + from the input (rectangular) grid that has true elevation values. Elevation calculations are done at each point on the polar grid and the results are returned. @@ -439,12 +443,11 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): Azimuth values for points that define the horizon profile. The ith element in this array corresponds to the ith element in horizon_angles. - horizon_angle: numeric - Dip angle values for points that define the horizon profile. The dip - angle of the horizon is the angle that the horizon makes with the - horizontal. It is given in degrees. If the horizon appears above - the horizontal, then the dip angle is positive. The ith element in - this array corresponds to the ith element in + horizon_angles: numeric + Elevation angle values for points that define the horizon profile. The + elevation angle of the horizon is the angle that the horizon makes with + the horizontal. It is given in degrees above the horizontal. The ith + element in this array corresponds to the ith element in horizon_azimuths. bin_size : int @@ -458,9 +461,9 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): filtered_angles. filtered_angles: numeric - Dip angle values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in - filtered_azimuths. + elevation angle values for points that define the horizon profile given + in degrees above the horizontal. The ith element in this array + corresponds to the ith element in filtered_azimuths. """ assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) @@ -468,13 +471,13 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): wedges = {} for i in range(horizon_angles.shape[0]): azimuth = horizon_azimuths[i] - dip = horizon_angles[i] + elevation = horizon_angles[i] azimuth_wedge = tools.round_to_nearest(azimuth, bin_size) if azimuth_wedge in wedges: - wedges[azimuth_wedge] = max(dip, wedges[azimuth_wedge]) + wedges[azimuth_wedge] = max(elevation, wedges[azimuth_wedge]) else: - wedges[azimuth_wedge] = dip + wedges[azimuth_wedge] = elevation filtered_azimuths = [] filtered_angles = [] @@ -487,10 +490,10 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): return filtered_azimuths, filtered_angles -def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): +def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ - Determine the dip angle created by the surface of a tilted plane - in a given direction. The dip angle is limited to be non-negative. + Determine the elevation angle created by the surface of a tilted plane + in a given direction. The angle is limited to be non-negative. Parameters ---------- @@ -505,15 +508,17 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): east of north (e.g. North = 0, South=180 East = 90, West = 270). direction : numeric - The direction along which the dip angle is to be calculated in + The direction along which the elevation angle is to be calculated in decimal degrees. The convention is defined as degrees east of north (e.g. North = 0, South=180 East = 90, West = 270). Returns -------- - dip_angle : numeric - The dip angle in direction created by the tilted plane. + elevation_angle : numeric + The angle between the surface of the tilted plane and the horizontal + when looking in the specified direction. Given in degrees above the + horizontal and limited to be non-negative. """ tilt = np.radians(surface_tilt) @@ -521,10 +526,10 @@ def collection_plane_dip_angle(surface_tilt, surface_azimuth, direction): declination = np.degrees(np.arctan(1.0/np.tan(tilt)/np.cos(bearing))) mask = (declination <= 0) - dip = 90.0 - declination - dip[mask] = 0.0 + elevation_angle = 90.0 - declination + elevation_angle[mask] = 0.0 - return dip + return elevation_angle def calculate_dtf(horizon_azimuths, horizon_angles, @@ -540,11 +545,10 @@ def calculate_dtf(horizon_azimuths, horizon_angles, element in this array corresponds to the ith element in horizon_angles. horizon_angles: numeric - Dip angle values for points that define the horizon profile. The dip - angle of the horizon is the angle that the horizon makes with the - horizontal. It is given in degrees. If the horizon appears above - the horizontal, then the dip angle is positive. The ith element in - this array corresponds to the ith element in + Elevation angle values for points that define the horizon profile. The + elevation angle of the horizon is the angle that the horizon makes with + the horizontal. It is given in degrees above the horizontal. The ith + element in this array corresponds to the ith element in horizon_azimuths. surface_tilt : numeric @@ -578,14 +582,14 @@ def calculate_dtf(horizon_azimuths, horizon_angles, num_points = horizon_azimuths.shape[0] for i in range(horizon_azimuths.shape[0]): az = np.radians(horizon_azimuths[i]) - horizon_dip = np.radians(horizon_angles[i]) - temp = np.radians(collection_plane_dip_angle(surface_tilt, - surface_azimuth, - horizon_azimuths[i])) - dip = np.maximum(horizon_dip, temp) + horizon_elev = np.radians(horizon_angles[i]) + temp = np.radians(collection_plane_elev_angle(surface_tilt, + surface_azimuth, + horizon_azimuths[i])) + elev = np.maximum(horizon_elev, temp) first_term = .5 * (a*np.cos(az) + b*np.sin(az)) * \ - (np.pi/2 - dip - np.sin(dip) * np.cos(dip)) - second_term = .5 * c * np.cos(dip)**2 + (np.pi/2 - elev - np.sin(elev) * np.cos(elev)) + second_term = .5 * c * np.cos(elev)**2 dtf += 2 * (first_term + second_term) / num_points return dtf diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 9049f00927..5c837d0087 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -19,7 +19,7 @@ def test_grid_lat_lon(): assert_allclose(lat_grid[30][72], lat + 22*grid_step) -def test_dip_calc(): +def test_elev_calc(): pt1 = np.array((71.23, -34.70, 1234)) pt2 = np.array((71.12, -34.16, 124)) pt3 = np.array((71.29, -35.23, 30044)) @@ -28,15 +28,16 @@ def test_dip_calc(): reverse_test_pts = np.vstack([pt2, pt1]) expected_bearings = np.array([121.9895, 302.5061]) - expected_dips = np.array([-2.8654, 2.6593]) + expected_elevs = np.array([-2.8654, 2.6593]) expected13 = np.array([[289.8132], [54.8663]]) - act_bearings, act_dips = horizon.dip_calc(test_pts, reverse_test_pts) + act_bearings, act_elevs = horizon.elevation_angle_calc(test_pts, + reverse_test_pts) assert_allclose(act_bearings, expected_bearings, rtol=1e-3) - assert_allclose(act_dips, expected_dips, rtol=1e-3) + assert_allclose(act_elevs, expected_elevs, rtol=1e-3) - actual13 = horizon.dip_calc(pt1, pt3) + actual13 = horizon.elevation_angle_calc(pt1, pt3) assert_allclose(expected13, actual13, rtol=1e-3) @@ -53,29 +54,29 @@ def test_calculate_horizon_points(): [212, 135, 1], [36, 145, 5]]) - bearings, dips = horizon.calculate_horizon_points(test_lat_grid, - test_lon_grid, - test_elev_grid, - sampling_method="grid") + dirs, elevs = horizon.calculate_horizon_points(test_lat_grid, + test_lon_grid, + test_elev_grid, + sampling_method="grid") - expected_bearings = np.array([0, 90, 45, 270, 135]) - rounded_bearings = np.round(bearings).astype(int) - assert(bearings.shape == dips.shape) - assert(np.all(np.in1d(expected_bearings, rounded_bearings))) + expected_dirs = np.array([0, 90, 45, 270, 135]) + rounded_dirs = np.round(dirs).astype(int) + assert(dirs.shape == elevs.shape) + assert(np.all(np.in1d(expected_dirs, rounded_dirs))) - bears, dips = horizon.calculate_horizon_points(test_lat_grid, + dirs, elevs = horizon.calculate_horizon_points(test_lat_grid, test_lon_grid, test_elev_grid, sampling_method="triangles", sampling_param=5) - assert(bears.shape == dips.shape) + assert(dirs.shape == elevs.shape) - bears, _ = horizon.calculate_horizon_points(test_lat_grid, - test_lon_grid, - test_elev_grid, - sampling_method="interpolator", - sampling_param=(10, 10)) - assert(bears.shape[0] == 100) + dirs, _ = horizon.calculate_horizon_points(test_lat_grid, + test_lon_grid, + test_elev_grid, + sampling_method="interpolator", + sampling_param=(10, 10)) + assert(dirs.shape[0] == 100) def test_sample_using_grid(): @@ -173,7 +174,7 @@ def test_filter_points(): assert(filtered_angles[1] == 10) -def test_collection_plane_dip_angle(): +def test_collection_plane_elev_angle(): surface_tilts = np.array([0, 5, 20, 38, 89]) surface_azimuths = np.array([0, 90, 180, 235, 355]) directions_easy = np.array([78, 270, 0, 145, 355]) @@ -182,15 +183,15 @@ def test_collection_plane_dip_angle(): expected_easy = np.array([0, 5, 20, 0, 0]) expected_hard = np.array([0, 3.21873120519, 10.3141048156, 21.3377447931, 0]) - dips_easy = horizon.collection_plane_dip_angle(surface_tilts, - surface_azimuths, - directions_easy) - assert_allclose(dips_easy, expected_easy) - - dips_hard = horizon.collection_plane_dip_angle(surface_tilts, - surface_azimuths, - directions_hard) - assert_allclose(dips_hard, expected_hard) + elevs_easy = horizon.collection_plane_elev_angle(surface_tilts, + surface_azimuths, + directions_easy) + assert_allclose(elevs_easy, expected_easy) + + elevs_hard = horizon.collection_plane_elev_angle(surface_tilts, + surface_azimuths, + directions_hard) + assert_allclose(elevs_hard, expected_hard) def test_calculate_dtf(): From 57c2035892895ba37bc776fd8650ed9bc8da16c7 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Wed, 7 Aug 2019 14:36:29 -0700 Subject: [PATCH 21/29] removed horizon.rst for now --- docs/sphinx/source/horizon.rst | 154 --------------------------------- pvlib/horizon.py | 3 +- 2 files changed, 2 insertions(+), 155 deletions(-) delete mode 100644 docs/sphinx/source/horizon.rst diff --git a/docs/sphinx/source/horizon.rst b/docs/sphinx/source/horizon.rst deleted file mode 100644 index c042434fbd..0000000000 --- a/docs/sphinx/source/horizon.rst +++ /dev/null @@ -1,154 +0,0 @@ -.. _horizon: - -Horizon -======== - -.. ipython:: python - :suppress: - - import pandas as pd - from pvlib import pvsystem - - -The :py:mod:`~pvlib.horizon` module contains many functions for horizon -profile modeling. - -The horizon profile at a location is a mapping from azimuth to elevation angle -(also called dip angle). - -A source of elevation data is needed for many of the -A common source of elevation data is NASA's SRTM [1]. An easily accessible source -of elevation data (albeit not free) is the googlemaps elevation API. There -are a few examples of how to query the googlemaps elevation API further below. - - -def fake_horizon_profile(max_dip): - """ - Creates a bogus horizon profile by randomly generating dip_angles at - integral azimuth values. Used for testing purposes. - """ - fake_profile = [] - for i in range(-180, 181): - fake_profile.append((i, random.random() * max_dip)) - return fake_profile - - - -def horizon_from_gmaps(lat, lon, GMAPS_API_KEY): - """ - Uses the functions defined in this modules to generate a complete horizon - profile for a location (specified by lat/lon). An API key for the - googlemaps elevation API is needeed. - """ - - grid = grid_lat_lon(lat, lon, grid_size=400, grid_step=.002) - elev_grid = grid_elevations_from_gmaps(grid, GMAPS_API_KEY) - horizon_points = calculate_horizon_points(elev_grid, - sampling_method="interpolator", - sampling_param=(1000, 1000)) - filtered_points = filter_points(horizon_points, bin_size=1) - return filtered_points - - - -def grid_elevations_from_gmaps(in_grid, GMAPS_API_KEY): - """ - Takes in a grid of lat lon values (shape: grid_size+1 x grid_size+1 x 2). - Queries the googlemaps elevation API to get elevation data at each lat/lon - point. Outputs the original grid with the elevation data appended along - the third axis so the shape is grid_size+1 x grid_size+1 x 3. - """ - - import googlemaps - - in_shape = in_grid.shape - lats = in_grid.T[0].flatten() - longs = in_grid.T[1].flatten() - locations = zip(lats, longs) - gmaps = googlemaps.Client(key=GMAPS_API_KEY) - - out_grid = np.ndarray((in_shape[0], in_shape[1], 3)) - - # Get elevation data from gmaps - elevations = [] - responses = [] - - while len(locations) > 512: - locations_to_request = locations[:512] - locations = locations[512:] - responses += gmaps.elevation(locations=locations_to_request) - responses += gmaps.elevation(locations=locations) - for entry in responses: - elevations.append(entry["elevation"]) - - for i in range(in_shape[0]): - for j in range(in_shape[1]): - lat = in_grid[i, j, 0] - lon = in_grid[i, j, 1] - elevation = elevations[i + j * in_shape[1]] - - out_grid[i, j, 0] = lat - out_grid[i, j, 1] = lon - out_grid[i, j, 2] = elevation - return out_grid - - - -def visualize(horizon_profile, pvsyst_scale=False): - """ - Plots a horizon profile with azimuth on the x-axis and dip angle on the y. - """ - azimuths = [] - dips = [] - for pair in horizon_profile: - azimuth = pair[0] - azimuths.append(azimuth) - dips.append(pair[1]) - plt.figure(figsize=(10, 6)) - if pvsyst_scale: - plt.ylim(0, 90) - plt.plot(azimuths, dips, "-") - plt.show - - -def polar_plot(horizon_profile): - """ - Plots a horizon profile on a polar plot with dip angle as the raidus and - azimuth as the theta value. An offset of 5 is added to the dip_angle to - make the plot more readable with low dip angles. - """ - azimuths = [] - dips = [] - for pair in horizon_profile: - azimuth = pair[0] - azimuths.append(np.radians(azimuth)) - dips.append(pair[1] + 5) - plt.figure(figsize=(10, 6)) - sp = plt.subplot(1, 1, 1, projection='polar') - sp.set_theta_zero_location('N') - sp.set_theta_direction(-1) - plt.plot(azimuths, dips, "o") - plt.show - -def invert_for_pvsyst(horizon_points, hemisphere="north"): - """ - Modify the azimuth values in horizon_points to match PVSyst's azimuth - convention (which is dependent on hemisphere) - """ - - # look at that northern hemisphere bias right there - # not even sorry. - assert hemisphere == "north" or hemisphere == "south" - - inverted_points = [] - for pair in horizon_points: - azimuth = pair[0] - if hemisphere == "north": - azimuth -= 180 - if azimuth < -180: - azimuth += 360 - elif hemisphere == "south": - azimuth = -azimuth - inverted_points.append((azimuth, pair[1])) - sorted_points = sorted(inverted_points, key=lambda tup: tup[0]) - return sorted_points \ No newline at end of file diff --git a/pvlib/horizon.py b/pvlib/horizon.py index c50ccd9dc2..0400c6503e 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -235,7 +235,8 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, samples_per_triangle=10): """ Creates triangles using nearest neighbors for every grid point and randomly - samples each of these triangles to find Elevation angles for the horizon profile. + samples each of these triangles to find elevation angles for the horizon + profile. Parameters ---------- From 00017e289ca254023948312ca50a66216d764df0 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 8 Aug 2019 11:23:38 -0700 Subject: [PATCH 22/29] added DNI correction to horizon.py --- pvlib/horizon.py | 50 ++++++++++++++++++++++++++++++++++++++ pvlib/test/test_horizon.py | 33 +++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 0400c6503e..8d3078968f 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -9,6 +9,8 @@ import itertools import numpy as np +import warnings + from scipy.interpolate import RegularGridInterpolator from pvlib import tools @@ -594,3 +596,51 @@ def calculate_dtf(horizon_azimuths, horizon_angles, second_term = .5 * c * np.cos(elev)**2 dtf += 2 * (first_term + second_term) / num_points return dtf + + +def DNI_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): + ''' + Calculates an adjustment to DNI based on a horizon profile. The adjustment + is a vector of binary values with the same length as the provided + solar position values. Where the sun is below the horizon, the adjustment + vector is 0 and it is 1 elsewhere. The horizon profile must be given as a + vector with 361 values where the ith value corresponds to the ith degree + of azimuth (0-360). + + + Parameters + ---------- + horizon_angles: numeric + Elevation angle values for points that define the horizon profile. The + elevation angle of the horizon is the angle that the horizon makes with + the horizontal. It is given in degrees above the horizontal. The ith + element in this array corresponds to the ith degree of azimuth. + + solar_zenith : numeric + Solar zenith angle. + + solar_azimuth : numeric + Solar azimuth angle. + + Returns + ------- + adjustment : numeric + A vector of binary values with the same shape as the inputted solar + position values. 0 when the sun is below the horizon and 1 elsewhere. + ''' + adjustment = np.ones(solar_zenith.shape) + + if (horizon_angles.shape[0] != 361): + warnings.warn('For DNI adjustment, horizon_angles needs to contain' + 'exactly 361 values (for each degree of azimuth 0-360).' + 'Since the provided horizon_angles contains {} values,' + 'no adjustment is calculated. A vector of ones is' + 'returned.'.format(horizon_angles.shape[0]), + UserWarning) + return adjustment + + rounded_solar_azimuth = np.round(solar_azimuth).astype(int) + horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] + mask = solar_zenith > horizon_zenith + adjustment[mask] = 0 + return adjustment diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 5c837d0087..27a079d582 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -1,7 +1,10 @@ import numpy as np +import pandas as pd from numpy.testing import assert_allclose from pvlib import horizon, tools +from test_irradiance import ephem_data, times + def test_grid_lat_lon(): grid_radius = 50 @@ -232,3 +235,33 @@ def test_calculate_dtf(): mask = np.logical_and((adjusted >= min_random_dtf), (adjusted <= max_random_dtf)) assert(np.all(mask)) + + +def test_DNI_horizon_adjustment(ephem_data): + zero_horizon = np.zeros(361) + max_horizon = np.full((361), 90.0) + + zero_adjusted = horizon.DNI_horizon_adjustment(zero_horizon, + ephem_data["zenith"], + ephem_data["azimuth"]) + zero_expected = np.array([0, 1, 1, 1]) + assert_allclose(zero_adjusted, zero_expected, atol=1e-7) + + max_adjusted = horizon.DNI_horizon_adjustment(max_horizon, + ephem_data["zenith"], + ephem_data["azimuth"]) + + max_expected = np.array([0, 0, 0, 0]) + assert_allclose(max_adjusted, max_expected, atol=1e-7) + + test_horizon = np.zeros(361) + test_horizon[145] = 75 + test_horizon[287] = 20 + test_horizon[67] = 10 + + adjusted = horizon.DNI_horizon_adjustment(test_horizon, + ephem_data["zenith"], + ephem_data["azimuth"]) + + expected = np.array([0, 0, 1, 0]) + assert_allclose(adjusted, expected, atol=1e-7) \ No newline at end of file From a87b79786f0308561533e9921290e6af8001fb33 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 8 Aug 2019 11:58:08 -0700 Subject: [PATCH 23/29] added a test case to get 100% of diff hit --- pvlib/test/test_horizon.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 27a079d582..68702b5e7d 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -239,29 +239,32 @@ def test_calculate_dtf(): def test_DNI_horizon_adjustment(ephem_data): zero_horizon = np.zeros(361) - max_horizon = np.full((361), 90.0) - + zero_expected = np.array([0, 1, 1, 1]) zero_adjusted = horizon.DNI_horizon_adjustment(zero_horizon, ephem_data["zenith"], ephem_data["azimuth"]) - zero_expected = np.array([0, 1, 1, 1]) assert_allclose(zero_adjusted, zero_expected, atol=1e-7) + max_horizon = np.full((361), 90.0) + max_expected = np.array([0, 0, 0, 0]) max_adjusted = horizon.DNI_horizon_adjustment(max_horizon, ephem_data["zenith"], ephem_data["azimuth"]) - - max_expected = np.array([0, 0, 0, 0]) assert_allclose(max_adjusted, max_expected, atol=1e-7) test_horizon = np.zeros(361) test_horizon[145] = 75 test_horizon[287] = 20 test_horizon[67] = 10 - + expected = np.array([0, 0, 1, 0]) adjusted = horizon.DNI_horizon_adjustment(test_horizon, ephem_data["zenith"], ephem_data["azimuth"]) + assert_allclose(adjusted, expected, atol=1e-7) - expected = np.array([0, 0, 1, 0]) - assert_allclose(adjusted, expected, atol=1e-7) \ No newline at end of file + bad_horizon = np.ones(360) + bad_expected = np.array([1, 1, 1, 1]) + bad_adjusted = horizon.DNI_horizon_adjustment(bad_horizon, + ephem_data["zenith"], + ephem_data["azimuth"]) + assert_allclose(bad_adjusted, bad_expected, atol=1e-7) From aeb5f955fe702774d017a671a979998153246dfe Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Thu, 8 Aug 2019 16:42:31 -0700 Subject: [PATCH 24/29] minor test fix --- pvlib/test/test_horizon.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 68702b5e7d..e9e030614f 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -156,9 +156,9 @@ def test_uniformly_sample_triangle(): alpha = 0.5 * np.linalg.norm(np.cross(p2-p, p3-p)) / area beta = 0.5 * np.linalg.norm(np.cross(p3-p, p1-p)) / area gamma = 1 - alpha - beta - assert(0 <= alpha <= 1) - assert(0 <= beta <= 1) - assert(0 <= gamma <= 1) + assert(-.02 <= alpha <= 1.02) + assert(-.02 <= beta <= 1.02) + assert(-.02 <= gamma <= 1.02) def test_filter_points(): From b022f9ee52e5e8b60141a5f8e0ae932d7c28f6e9 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Fri, 9 Aug 2019 14:25:18 -0700 Subject: [PATCH 25/29] docstring changes and improvement to filter_points --- pvlib/horizon.py | 193 +++++++++++++++++++------------------ pvlib/test/test_horizon.py | 53 +++++----- 2 files changed, 124 insertions(+), 122 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 8d3078968f..63c03a926e 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -11,8 +11,6 @@ import numpy as np import warnings -from scipy.interpolate import RegularGridInterpolator - from pvlib import tools @@ -30,6 +28,14 @@ def grid_lat_lon(lat, lon, grid_radius=200, grid_step=.001): lon : numeric The longitude of the location that is to be the center of the grid. + grid_radius : numeric + The number of points to generate in each of the cardinal directions + for the grid. The resulting grid will be a square with + (2xgrid_radius)+1 points along each side. Unitless. + + grid_step : numeric + The degrees of latitude/longitude between adjacent points in the grid. + Returns ------- lat_grid: 2d-array @@ -54,26 +60,31 @@ def grid_lat_lon(lat, lon, grid_radius=200, grid_step=.001): return lat_grid, lon_grid -def elevation_angle_calc(pt1, pt2): +def elevation_and_azimuth(pt1, pt2): ''' - Calculates the elevation angle from pt1 to pt2 where elevation angle is + Calculates the elevation angle from pt1 to pt2. Elevation angle is defined as the angle between the line connecting pt1 to pt2 and the plane - normal to the Earth's surface at pt1. A point that appears above the + tangent to the Earth's surface at pt1. A point that appears above the horizontal has a positive elevation angle. Also computes the azimuth - defined as degrees East of North the bearing of pt2 from pt1. + defined as degrees East of North of the bearing from pt1 to pt2. This uses the Haversine formula. Parameters ---------- pt1 : ndarray - Nx3 array that contains lat, lon, and elev values that correspond - to the origin points from which the elevation angles are to be - calculated. The observer points. + Nx3 array that contains latitude, longitude, and elevation values + that correspond to the origin (observer) points from which the + elevation angles are to be calculated. Longitude should be given in + degrees East of the Prime Meridian and latitude in degrees North of the + Equator. Units are [deg, deg, meters] + pt2 : ndarray - Nx3 array that contains lat, lon, and elev values that correspond - to the target points to which the elevation angles are to be - calculated. The observee points. + Nx3 array that contains latitude, longitude, and elevation values + that correspond to the target (observee) points to which the elevation + angles are to be calculated. Longitude should be given in + degrees East of the Prime Meridian and latitude in degrees North of the + Equator. Units are [deg, deg, meters] Returns ------- @@ -132,74 +143,74 @@ def elevation_angle_calc(pt1, pt2): def calculate_horizon_points(lat_grid, lon_grid, elev_grid, - sampling_method="grid", sampling_param=400): + sampling_method="grid", num_samples=None): """ - Calculates a horizon profile from a three grids containing lat, lon, - and elevation values. The "site" is assumed to be at the center of the - grid. + Calculates a horizon profile viewed from the center of a latitude, + longitude, elevation grid. Parameters ---------- lat_grid : ndarray - A 2d array containing latitude values that correspond to the other - two input grids. + A 2d array of latitude values for the grid. lon_grid : ndarray - A 2d array containing longitude values that correspond to the other - two input grids. + A 2d array of longitude values for the grid. elev_grid : ndarray - A 2d array containing elevation values that correspond to the other - two input grids. + A 2d array of elevation values for the grid. - sampling_method : string + sampling_method : string, default "grid" A string that specifies the sampling method used to generate the - horizon profile. + horizon profile. Acceptable values are: "grid", "triangles", "polar". - sampling_param : variable + num_samples : variable, default None A parameter that is passed into the function specified by sampling_method. + If the sampling method is "triangles" this corresponds + to the number of samples taken from each triangle. + See _sampling_using_triangles for more info. + If the sampling method is "polar" this should be a tuple with 2 values + that define the number of points along each polar axis to sample. + See _sampling_using_interpolator for more info. Returns ------- bearing_deg: Nx1 ndarray - The bearings from the "site" to sampled points in degrees - East of North. + The bearings from the grid_center to horizon points. elevation_angle_deg: numeric The angles that the sampled points make with the horizontal - as observed from the "site". Given in degrees above the + as observed from the grid center. Given in degrees above the horizontal. """ - assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) grid_shape = lat_grid.shape grid_center_i = (grid_shape[0] - 1) // 2 grid_center_j = (grid_shape[1] - 1) // 2 - site_lat = lat_grid[grid_center_i, grid_center_j] - site_lon = lon_grid[grid_center_i, grid_center_j] - site_elev = elev_grid[grid_center_i, grid_center_j] - site = np.array([site_lat, site_lon, site_elev]) + center_lat = lat_grid[grid_center_i, grid_center_j] + center_lon = lon_grid[grid_center_i, grid_center_j] + center_elev = elev_grid[grid_center_i, grid_center_j] + center = np.array([center_lat, center_lon, center_elev]) if sampling_method == "grid": - samples = sample_using_grid(lat_grid, lon_grid, elev_grid) + samples = _sample_using_grid(lat_grid, lon_grid, elev_grid) elif sampling_method == "triangles": - samples = sample_using_triangles(lat_grid, lon_grid, elev_grid, - sampling_param) - elif sampling_method == "interpolator": - samples = sample_using_interpolator(lat_grid, lon_grid, elev_grid, - sampling_param) + samples = _sample_using_triangles(lat_grid, lon_grid, elev_grid, + num_samples) + elif sampling_method == "polar": + samples = _sample_using_interpolator(lat_grid, lon_grid, elev_grid, + num_samples) - bearing_deg, elevation_angle_deg = elevation_angle_calc(site, samples) + bearing_deg, elevation_angle_deg = elevation_and_azimuth(center, samples) return bearing_deg, elevation_angle_deg -def sample_using_grid(lat_grid, lon_grid, elev_grid): +def _sample_using_grid(lat_grid, lon_grid, elev_grid): """ - Calculates the Elevation angle from the site (center of the grid) - to every point on the grid. + Returns every point on the grid excluding the grid center as samples + for horizon calculations. Parameters ---------- @@ -221,24 +232,22 @@ def sample_using_grid(lat_grid, lon_grid, elev_grid): all_samples: Nx3 ndarray Array of lat, lon, elev points that are grid points. """ - assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) lats = lat_grid.flatten() lons = lon_grid.flatten() elevs = elev_grid.flatten() samples = np.stack([lats, lons, elevs], axis=1) - # remove site from samples + # remove grid center from samples all_samples = np.delete(samples, samples.shape[0]//2, axis=0) return all_samples -def sample_using_triangles(lat_grid, lon_grid, elev_grid, - samples_per_triangle=10): +def _sample_using_triangles(lat_grid, lon_grid, elev_grid, + samples_per_triangle=10): """ Creates triangles using nearest neighbors for every grid point and randomly - samples each of these triangles to find elevation angles for the horizon - profile. + samples each of these triangles. Parameters ---------- @@ -265,7 +274,6 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf """ - assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) # start with empty array all_samples = np.array([], dtype=np.float64).reshape(0, 3) @@ -329,13 +337,11 @@ def sample_using_triangles(lat_grid, lon_grid, elev_grid, return np.array(all_samples) -def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): +def _sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): """ Creates a "grid" using polar coordinates and uses the scipy's grid interpolator to estimate elevation values at each point on the polar grid - from the input (rectangular) grid that has true elevation values. Elevation - calculations are done at each point on the polar grid and the results - are returned. + from the input (rectangular) grid that has true elevation values. Parameters ---------- @@ -363,7 +369,6 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): Array of [lat, lon, elev] points that were sampled using the polar grid. """ - assert(lat_grid.shape == lon_grid.shape == elev_grid.shape) lats = lat_grid[0] lons = lon_grid.T[0] @@ -373,8 +378,13 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): grid_shape = lat_grid.shape grid_center_i = (grid_shape[0] - 1) // 2 grid_center_j = (grid_shape[1] - 1) // 2 - site_lat = lat_grid[grid_center_i, grid_center_j] - site_lon = lon_grid[grid_center_i, grid_center_j] + center_lat = lat_grid[grid_center_i, grid_center_j] + center_lon = lon_grid[grid_center_i, grid_center_j] + + try: + from scipy.interpolate import RegularGridInterpolator + except ImportError: + raise ImportError('The polar sampling function requires scipy') interpolator = RegularGridInterpolator((lats, lons), elev_grid.T) @@ -383,7 +393,7 @@ def sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): polar_pts = np.array(list(itertools.product(r, theta))) pts = np.array([tools.polar_to_cart(e[0], e[1]) for e in polar_pts]) - pts += np.array((site_lat, site_lon)) + pts += np.array((center_lat, center_lon)) total_num_samples = num_samples[0]*num_samples[1] interpolated_elevs = interpolator(pts).reshape(total_num_samples, 1) @@ -434,27 +444,28 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): return random_pts -def filter_points(horizon_azimuths, horizon_angles, bin_size=1): +def filter_points(azimuths, elevation_angles, bin_size=1): """ - Bins the horizon_points by azimuth values. The azimuth value of each - point in horizon_points is rounded to the nearest bin and then the - max value in each bin is returned. + Bins the horizon points by azimuth values. The azimuth value of each + point is rounded to the nearest bin and then the + max elevation angle in each bin is returned. Parameters ---------- - horizon_azimuths: numeric + azimuths: numeric Azimuth values for points that define the horizon profile. The ith - element in this array corresponds to the ith element in horizon_angles. + element in this array corresponds to the ith element in + elevation_angles. - horizon_angles: numeric + elevation_angles: numeric Elevation angle values for points that define the horizon profile. The elevation angle of the horizon is the angle that the horizon makes with the horizontal. It is given in degrees above the horizontal. The ith element in this array corresponds to the ith element in - horizon_azimuths. + azimuths. bin_size : int - The width of the bins for the azimuth values. + The width of the bins for the azimuth values. (degrees) Returns ------- @@ -469,34 +480,25 @@ def filter_points(horizon_azimuths, horizon_angles, bin_size=1): corresponds to the ith element in filtered_azimuths. """ - assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) + assert(azimuths.shape[0] == elevation_angles.shape[0]) - wedges = {} - for i in range(horizon_angles.shape[0]): - azimuth = horizon_azimuths[i] - elevation = horizon_angles[i] - azimuth_wedge = tools.round_to_nearest(azimuth, bin_size) + rounded_azimuths = tools.round_to_nearest(azimuths, bin_size) + bins = np.unique(rounded_azimuths) - if azimuth_wedge in wedges: - wedges[azimuth_wedge] = max(elevation, wedges[azimuth_wedge]) - else: - wedges[azimuth_wedge] = elevation + filtered = np.column_stack((bins, np.nan * bins)) - filtered_azimuths = [] - filtered_angles = [] - for key in sorted(wedges.keys()): - filtered_azimuths.append(key) - filtered_angles.append(wedges[key]) + for i in range(filtered.shape[0]): + idx = (rounded_azimuths == filtered[i, 0]) + filtered[i, 1] = np.max(elevation_angles[idx]) - filtered_angles = np.array(filtered_angles) - filtered_azimuths = np.array(filtered_azimuths) - return filtered_azimuths, filtered_angles + return filtered[:, 0], filtered[:, 1] def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ Determine the elevation angle created by the surface of a tilted plane - in a given direction. The angle is limited to be non-negative. + intersecting the plane tangent to the Earth's surface in a given direction. + The angle is limited to be non-negative. Parameters ---------- @@ -598,14 +600,14 @@ def calculate_dtf(horizon_azimuths, horizon_angles, return dtf -def DNI_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): +def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' - Calculates an adjustment to DNI based on a horizon profile. The adjustment - is a vector of binary values with the same length as the provided - solar position values. Where the sun is below the horizon, the adjustment - vector is 0 and it is 1 elsewhere. The horizon profile must be given as a - vector with 361 values where the ith value corresponds to the ith degree - of azimuth (0-360). + Calculates an adjustment to direct normal irradiance based on a horizon + profile. The adjustment is a vector of binary values with the same length + as the provided solar position values. Where the sun is below the horizon, + the adjustment vector is 0 and it is 1 elsewhere. The horizon profile must + be given as a vector with 360 values where the ith value corresponds to the + ith degree of azimuth (0-359). Parameters @@ -630,9 +632,9 @@ def DNI_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' adjustment = np.ones(solar_zenith.shape) - if (horizon_angles.shape[0] != 361): + if (horizon_angles.shape[0] != 360): warnings.warn('For DNI adjustment, horizon_angles needs to contain' - 'exactly 361 values (for each degree of azimuth 0-360).' + 'exactly 360 values (for each degree of azimuth 0-359).' 'Since the provided horizon_angles contains {} values,' 'no adjustment is calculated. A vector of ones is' 'returned.'.format(horizon_angles.shape[0]), @@ -640,6 +642,7 @@ def DNI_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): return adjustment rounded_solar_azimuth = np.round(solar_azimuth).astype(int) + rounded_solar_azimuth[rounded_solar_azimuth == 360] = 0 horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] mask = solar_zenith > horizon_zenith adjustment[mask] = 0 diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index e9e030614f..80a0e4efa2 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -1,5 +1,4 @@ import numpy as np -import pandas as pd from numpy.testing import assert_allclose from pvlib import horizon, tools @@ -35,12 +34,12 @@ def test_elev_calc(): expected13 = np.array([[289.8132], [54.8663]]) - act_bearings, act_elevs = horizon.elevation_angle_calc(test_pts, - reverse_test_pts) + act_bearings, act_elevs = horizon.elevation_and_azimuth(test_pts, + reverse_test_pts) assert_allclose(act_bearings, expected_bearings, rtol=1e-3) assert_allclose(act_elevs, expected_elevs, rtol=1e-3) - actual13 = horizon.elevation_angle_calc(pt1, pt3) + actual13 = horizon.elevation_and_azimuth(pt1, pt3) assert_allclose(expected13, actual13, rtol=1e-3) @@ -71,14 +70,14 @@ def test_calculate_horizon_points(): test_lon_grid, test_elev_grid, sampling_method="triangles", - sampling_param=5) + num_samples=5) assert(dirs.shape == elevs.shape) dirs, _ = horizon.calculate_horizon_points(test_lat_grid, test_lon_grid, test_elev_grid, - sampling_method="interpolator", - sampling_param=(10, 10)) + sampling_method="polar", + num_samples=(10, 10)) assert(dirs.shape[0] == 100) @@ -95,9 +94,9 @@ def test_sample_using_grid(): [212, 135, 1], [36, 145, 5]]) - samples = horizon.sample_using_grid(test_lat_grid, - test_lon_grid, - test_elev_grid) + samples = horizon._sample_using_grid(test_lat_grid, + test_lon_grid, + test_elev_grid) assert(len(samples) == 8) @@ -113,10 +112,10 @@ def test_sample_using_triangles(): test_elev_grid = np.array([[15, 18, 43], [212, 135, 1], [36, 145, 5]]) - samples = horizon.sample_using_triangles(test_lat_grid, - test_lon_grid, - test_elev_grid, - samples_per_triangle=2) + samples = horizon._sample_using_triangles(test_lat_grid, + test_lon_grid, + test_elev_grid, + samples_per_triangle=2) assert(len(samples) == 32) @@ -132,10 +131,10 @@ def test_using_interpolator(): test_elev_grid = np.array([[15, 18, 43], [212, 135, 1], [36, 145, 5]]) - samples = horizon.sample_using_interpolator(test_lat_grid, - test_lon_grid, - test_elev_grid, - num_samples=(5, 5)) + samples = horizon._sample_using_interpolator(test_lat_grid, + test_lon_grid, + test_elev_grid, + num_samples=(5, 5)) assert(len(samples) == 25) @@ -237,34 +236,34 @@ def test_calculate_dtf(): assert(np.all(mask)) -def test_DNI_horizon_adjustment(ephem_data): - zero_horizon = np.zeros(361) +def test_dni_horizon_adjustment(ephem_data): + zero_horizon = np.zeros(360) zero_expected = np.array([0, 1, 1, 1]) - zero_adjusted = horizon.DNI_horizon_adjustment(zero_horizon, + zero_adjusted = horizon.dni_horizon_adjustment(zero_horizon, ephem_data["zenith"], ephem_data["azimuth"]) assert_allclose(zero_adjusted, zero_expected, atol=1e-7) - max_horizon = np.full((361), 90.0) + max_horizon = np.full((360), 90.0) max_expected = np.array([0, 0, 0, 0]) - max_adjusted = horizon.DNI_horizon_adjustment(max_horizon, + max_adjusted = horizon.dni_horizon_adjustment(max_horizon, ephem_data["zenith"], ephem_data["azimuth"]) assert_allclose(max_adjusted, max_expected, atol=1e-7) - test_horizon = np.zeros(361) + test_horizon = np.zeros(360) test_horizon[145] = 75 test_horizon[287] = 20 test_horizon[67] = 10 expected = np.array([0, 0, 1, 0]) - adjusted = horizon.DNI_horizon_adjustment(test_horizon, + adjusted = horizon.dni_horizon_adjustment(test_horizon, ephem_data["zenith"], ephem_data["azimuth"]) assert_allclose(adjusted, expected, atol=1e-7) - bad_horizon = np.ones(360) + bad_horizon = np.ones(361) bad_expected = np.array([1, 1, 1, 1]) - bad_adjusted = horizon.DNI_horizon_adjustment(bad_horizon, + bad_adjusted = horizon.dni_horizon_adjustment(bad_horizon, ephem_data["zenith"], ephem_data["azimuth"]) assert_allclose(bad_adjusted, bad_expected, atol=1e-7) From d07e66f8b8084a7f461814b17d3aec058d1e0eb4 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Mon, 12 Aug 2019 13:33:54 -0700 Subject: [PATCH 26/29] added tests for functions added in tools.py. Some docstring changes as well. --- pvlib/horizon.py | 14 ++++++++------ pvlib/test/test_tools.py | 37 +++++++++++++++++++++++++++++++++++++ pvlib/tools.py | 18 ++++++++++++++---- 3 files changed, 59 insertions(+), 10 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 63c03a926e..20f2cd2a3b 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,8 +1,6 @@ """ The ``horizon`` module contains functions for horizon profile modeling. -There are various geometric utilities that are useful in horizon calculations -as well as a method that uses the googlemaps elevation API to create a -horizon profile. +There are various geometric utilities that are useful in horizon calculations. """ from __future__ import division @@ -68,6 +66,7 @@ def elevation_and_azimuth(pt1, pt2): horizontal has a positive elevation angle. Also computes the azimuth defined as degrees East of North of the bearing from pt1 to pt2. This uses the Haversine formula. + The method used to calculate the elevation angle is discussed in [1]. Parameters ---------- @@ -97,7 +96,9 @@ def elevation_and_azimuth(pt1, pt2): horizontal. ''' + # Equatorial Radius of the Earth (ellipsoid model) in meters a = 6378137.0 + # Polar Radius of the Earth (ellipsoid model) in meters b = 6356752.0 lat1 = np.atleast_1d(pt1.T[0]) @@ -341,7 +342,7 @@ def _sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): """ Creates a "grid" using polar coordinates and uses the scipy's grid interpolator to estimate elevation values at each point on the polar grid - from the input (rectangular) grid that has true elevation values. + from the input (rectangular) grid that has true elevation values. Parameters ---------- @@ -447,8 +448,9 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): def filter_points(azimuths, elevation_angles, bin_size=1): """ Bins the horizon points by azimuth values. The azimuth value of each - point is rounded to the nearest bin and then the - max elevation angle in each bin is returned. + point is rounded to the nearest bin and then the max elevation angle + in each bin is returned. The bins will have azimuth values of n*bin_size + where n is some integer. Parameters ---------- diff --git a/pvlib/test/test_tools.py b/pvlib/test/test_tools.py index 42eaedca58..a2d1a9622e 100644 --- a/pvlib/test/test_tools.py +++ b/pvlib/test/test_tools.py @@ -1,5 +1,7 @@ import pytest +import numpy as np +from numpy.testing import assert_allclose from pvlib import tools @@ -12,3 +14,38 @@ def test_build_kwargs(keys, input_dict, expected): kwargs = tools._build_kwargs(keys, input_dict) assert kwargs == expected + + +def test_latitude_conversions(): + geocentric_lats = np.radians(np.array([0, 30.2, 45, 90])) + geodetic_lats = np.radians(np.array([0, 30.36759, 45.19243, 90])) + temp = tools.latitude_to_geodetic(geocentric_lats) + assert_allclose(np.degrees(temp), np.degrees(geodetic_lats)) + temp2 = tools.latitude_to_geocentric(geodetic_lats) + assert_allclose(np.degrees(temp2), np.degrees(geocentric_lats)) + + +def test_basis_conversions(): + test_lle = np.array([[30, -45, 900], + [20, 85, 10], + [-50, 45, 5]]) + test_xyz = np.array([[3909619.9, -3909619.9, 3170821.2], + [522572.4, 5973029.8, 2167700], + [2904700.9, 2904700.9, -4862792.5]]) + + actual_xyz = tools.lle_to_xyz(test_lle) + assert_allclose(actual_xyz, test_xyz) + + actual_lle = tools.xyz_to_lle(test_xyz) + assert_allclose(actual_lle, test_lle, rtol=1e-2) + + +def test_polar_to_cart(): + test_rho = np.array([10, 10, 50, 20]) + test_phi = np.radians(np.array([180, -30, 45, 270])) + expected_x = np.array([-10, 5*3**.5, 50.0/2**.5, 0]) + expected_y = np.array([0, -5, 50.0/2**.5, -20]) + actual_x, actual_y = tools.polar_to_cart(test_rho, + test_phi) + assert_allclose(actual_x, expected_x, atol=1e-7) + assert_allclose(actual_y, expected_y, atol=1e-7) diff --git a/pvlib/tools.py b/pvlib/tools.py index fcd13c3238..c68dfa0845 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -430,20 +430,26 @@ def _golden_sect_DataFrame(params, VL, VH, func): def latitude_to_geocentric(phi): """ Converts a geodetic (common) latitude to a geocentric latitude. + Latitude must be given in radians. [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf """ - a = 6378.137 - b = 6356.752 + # Equatorial Radius of the Earth (ellipsoid model) in meters + a = 6378137.0 + # Polar Radius of the Earth (ellipsoid model) in meters + b = 6356752.0 return np.arctan(b**2/a**2*np.tan(phi)) def latitude_to_geodetic(phi): """ Converts a geocentric latitude to a geodetic (common) latitude. + Latitude must be given in radians. [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf """ - a = 6378.137 - b = 6356.752 + # Equatorial Radius of the Earth (ellipsoid model) in meters + a = 6378137.0 + # Polar Radius of the Earth (ellipsoid model) in meters + b = 6356752.0 return np.arctan(a**2/b**2*np.tan(phi)) @@ -457,7 +463,9 @@ def lle_to_xyz(point): lon = np.atleast_1d(point.T[1]) elev = np.atleast_1d(point.T[2]) + # Equatorial Radius of the Earth (ellipsoid model) in meters a = 6378137.0 + # Polar Radius of the Earth (ellipsoid model) in meters b = 6356752.0 # convert to radians @@ -485,7 +493,9 @@ def xyz_to_lle(point): The center of the earth is the origin in the xyz-space. The output latitude is assumed to be a common latitude (geodetic). """ + # Equatorial Radius of the Earth (ellipsoid model) in meters a = 6378137.0 + # Polar Radius of the Earth (ellipsoid model) in meters b = 6356752.0 x = np.atleast_1d(point.T[0]) From 24115201cf50c09d95886ba6ccd9c38e537ff2f6 Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 13 Aug 2019 14:11:34 -0700 Subject: [PATCH 27/29] minor improvements and docstring changes --- pvlib/horizon.py | 130 +++++++++++++++++++++++-------------- pvlib/test/test_horizon.py | 19 +++--- pvlib/tools.py | 36 +++++----- setup.py | 6 +- 4 files changed, 112 insertions(+), 79 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 20f2cd2a3b..18a4879541 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -7,7 +7,6 @@ import itertools import numpy as np -import warnings from pvlib import tools @@ -77,7 +76,6 @@ def elevation_and_azimuth(pt1, pt2): degrees East of the Prime Meridian and latitude in degrees North of the Equator. Units are [deg, deg, meters] - pt2 : ndarray Nx3 array that contains latitude, longitude, and elevation values that correspond to the target (observee) points to which the elevation @@ -95,18 +93,23 @@ def elevation_and_azimuth(pt1, pt2): as observed from the points in pt1. Given in degrees above the horizontal. + Examples + ________ + site_loc = np.array([[37, 34, 100]]) + target_locs = np.array([[38, 34, 63], + [36, 35, 231], + [36, 35, 21]]) + bearing, elev_angles = elevation_and_azimuth(site_loc, target_locs) ''' # Equatorial Radius of the Earth (ellipsoid model) in meters a = 6378137.0 # Polar Radius of the Earth (ellipsoid model) in meters b = 6356752.0 - lat1 = np.atleast_1d(pt1.T[0]) - lon1 = np.atleast_1d(pt1.T[1]) - elev1 = np.atleast_1d(pt1.T[2]) - lat2 = np.atleast_1d(pt2.T[0]) - lon2 = np.atleast_1d(pt2.T[1]) - elev2 = np.atleast_1d(pt2.T[2]) + lat1 = pt1.T[0] + lon1 = pt1.T[1] + lat2 = pt2.T[0] + lon2 = pt2.T[1] # convert to radians phi1 = np.radians(lat1) @@ -114,15 +117,17 @@ def elevation_and_azimuth(pt1, pt2): phi2 = np.radians(lat2) theta2 = np.radians(lon2) - v1 = tools.lle_to_xyz(np.stack([lat1, lon1, elev1], axis=1)) - v2 = tools.lle_to_xyz(np.stack([lat2, lon2, elev2], axis=1)) - x1 = np.atleast_1d(v1.T[0]) - y1 = np.atleast_1d(v1.T[1]) - z1 = np.atleast_1d(v1.T[2]) - - delta = np.atleast_2d(np.subtract(v1, v2)) + v1 = tools.lle_to_xyz(pt1) + v2 = tools.lle_to_xyz(pt2) + x1 = v1.T[0] + y1 = v1.T[1] + z1 = v1.T[2] - normal = np.atleast_2d(np.stack([2*x1/a**2, 2*y1/a**2, 2*z1/b**2], axis=1)) + delta = np.subtract(v1, v2) + a_sqrd = a**2 + b_sqrd = b**2 + normal = 2 * np.stack([x1/a_sqrd, y1/a_sqrd, z1/b_sqrd], + axis=1) # Take the dot product of corresponding vectors dot = np.sum(np.multiply(delta, normal), axis=1) @@ -163,6 +168,7 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, sampling_method : string, default "grid" A string that specifies the sampling method used to generate the horizon profile. Acceptable values are: "grid", "triangles", "polar". + See Notes for brief descriptions of each. num_samples : variable, default None A parameter that is passed into the function specified by @@ -172,7 +178,7 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, See _sampling_using_triangles for more info. If the sampling method is "polar" this should be a tuple with 2 values that define the number of points along each polar axis to sample. - See _sampling_using_interpolator for more info. + See Notes for more info. Returns ------- @@ -184,6 +190,20 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, as observed from the grid center. Given in degrees above the horizontal. + Notes + _____ + Sampling methods: + "grid" - Uses every point on the grid exclusing the grid + center as samples for hrizon calculations. + + "triangles" - Creates triangles using nearest neighbors for + every grid point and randomly samples the surface of each of these + triangles. num_samples sets the number of samples taken from each triangle. + + "polar" - Creates a polar "grid" and uses scipy's grid interpolator to + estimate elevation values at each point on the polar grid from the true + elevation data. num_samples sets the number of points along each polar + axis (radial and angular). """ grid_shape = lat_grid.shape @@ -202,6 +222,8 @@ def calculate_horizon_points(lat_grid, lon_grid, elev_grid, elif sampling_method == "polar": samples = _sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples) + else: + raise ValueError('Invalid sampling method: %s', sampling_method) bearing_deg, elevation_angle_deg = elevation_and_azimuth(center, samples) @@ -273,7 +295,7 @@ def _sample_using_triangles(lat_grid, lon_grid, elev_grid, all_samples: Nx3 ndarray Array of [lat, lon, elev] points that were sampled from the grid. - [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf + [1] Osada et al. (2002) ACM Transactions on Graphics. 21(4) 807-832 """ # start with empty array @@ -291,10 +313,10 @@ def _sample_using_triangles(lat_grid, lon_grid, elev_grid, top = np.array([lat_grid[i-1, j], lon_grid[i-1, j], elev_grid[i-1, j]]) - samples = uniformly_sample_triangle(center, - top, - left, - samples_per_triangle) + samples = _uniformly_sample_triangle(center, + top, + left, + samples_per_triangle) all_samples = np.vstack([all_samples, samples]) if i != 0 and j != lat_grid.shape[1] - 1: @@ -304,10 +326,10 @@ def _sample_using_triangles(lat_grid, lon_grid, elev_grid, top = np.array([lat_grid[i-1, j], lon_grid[i-1, j], elev_grid[i-1, j]]) - samples = uniformly_sample_triangle(center, - top, - right, - samples_per_triangle) + samples = _uniformly_sample_triangle(center, + top, + right, + samples_per_triangle) all_samples = np.vstack([all_samples, samples]) if i != lat_grid.shape[0] - 1 and j != 0: @@ -317,10 +339,10 @@ def _sample_using_triangles(lat_grid, lon_grid, elev_grid, bottom = np.array([lat_grid[i+1, j], lon_grid[i+1, j], elev_grid[i+1, j]]) - samples = uniformly_sample_triangle(center, - bottom, - left, - samples_per_triangle) + samples = _uniformly_sample_triangle(center, + bottom, + left, + samples_per_triangle) all_samples = np.vstack([all_samples, samples]) if i != lat_grid.shape[0] - 1 and j != lat_grid.shape[1] - 1: @@ -330,17 +352,17 @@ def _sample_using_triangles(lat_grid, lon_grid, elev_grid, bottom = np.array([lat_grid[i+1, j], lon_grid[i+1, j], elev_grid[i+1, j]]) - samples = uniformly_sample_triangle(center, - bottom, - right, - samples_per_triangle) + samples = _uniformly_sample_triangle(center, + bottom, + right, + samples_per_triangle) all_samples = np.vstack([all_samples, samples]) return np.array(all_samples) def _sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): """ - Creates a "grid" using polar coordinates and uses the scipy's grid + Creates a "grid" using polar coordinates and uses scipy's grid interpolator to estimate elevation values at each point on the polar grid from the input (rectangular) grid that has true elevation values. @@ -402,7 +424,7 @@ def _sample_using_interpolator(lat_grid, lon_grid, elev_grid, num_samples): return samples -def uniformly_sample_triangle(p1, p2, p3, num_samples): +def _uniformly_sample_triangle(p1, p2, p3, num_samples): """ Randomly sample the surface of a triangle defined by three (lat, lon, elev) points uniformly [1]. @@ -429,7 +451,7 @@ def uniformly_sample_triangle(p1, p2, p3, num_samples): Array with N (lat, lon, elev) points that lie on the surface of the triangle. - [1] http://graphics.stanford.edu/courses/cs468-08-fall/pdf/osada.pdf + [1] Osada et al. (2002) ACM Transactions on Graphics. 21(4) 807-832 """ c1 = tools.lle_to_xyz(p1) c2 = tools.lle_to_xyz(p2) @@ -482,7 +504,9 @@ def filter_points(azimuths, elevation_angles, bin_size=1): corresponds to the ith element in filtered_azimuths. """ - assert(azimuths.shape[0] == elevation_angles.shape[0]) + if azimuths.shape[0] != elevation_angles.shape[0]: + raise ValueError('azimuths and elevation_angles must be of the same' + 'length.') rounded_azimuths = tools.round_to_nearest(azimuths, bin_size) bins = np.unique(rounded_azimuths) @@ -500,7 +524,7 @@ def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ Determine the elevation angle created by the surface of a tilted plane intersecting the plane tangent to the Earth's surface in a given direction. - The angle is limited to be non-negative. + The angle is limited to be non-negative. This comes from Equation 10 in [1] Parameters ---------- @@ -527,6 +551,9 @@ def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): when looking in the specified direction. Given in degrees above the horizontal and limited to be non-negative. + + [1] doi.org/10.1016/j.solener.2014.09.037 + """ tilt = np.radians(surface_tilt) bearing = np.radians(direction - surface_azimuth - 180.0) @@ -542,8 +569,9 @@ def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): def calculate_dtf(horizon_azimuths, horizon_angles, surface_tilt, surface_azimuth): """ - Calculate the diffuse tilt factor that is adjusted with the horizon - profile. + Calculate the diffuse tilt factor for a tilted plane that is adjusted + with for horizon profile. The idea for a diffuse tilt factor is explained + in [1]. Parameters ---------- @@ -576,6 +604,17 @@ def calculate_dtf(horizon_azimuths, horizon_angles, the sky that is adjusted for the horizon profile and the tilt of the plane. + Notes + _____ + + The dtf in this method is calculated by approximating the surface integral + over the visible section of the sky dome. The integrand of the surface + integral is the cosine of the angle between the incoming radiation and the + vector normal to the surface. The method calculates a sum of integrations + from the "peak" of the sky dome down to the elevation angle of the horizon. + + [1] Goss et al. (2014) Solar Energy 110, 410-419 + """ assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) tilt_rad = np.radians(surface_tilt) @@ -635,13 +674,8 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): adjustment = np.ones(solar_zenith.shape) if (horizon_angles.shape[0] != 360): - warnings.warn('For DNI adjustment, horizon_angles needs to contain' - 'exactly 360 values (for each degree of azimuth 0-359).' - 'Since the provided horizon_angles contains {} values,' - 'no adjustment is calculated. A vector of ones is' - 'returned.'.format(horizon_angles.shape[0]), - UserWarning) - return adjustment + raise ValueError('horizon_angles must contain exactly 360 values' + '(for each degree of azimuth 0-359).') rounded_solar_azimuth = np.round(solar_azimuth).astype(int) rounded_solar_azimuth[rounded_solar_azimuth == 360] = 0 diff --git a/pvlib/test/test_horizon.py b/pvlib/test/test_horizon.py index 80a0e4efa2..c6ff856c65 100644 --- a/pvlib/test/test_horizon.py +++ b/pvlib/test/test_horizon.py @@ -22,9 +22,9 @@ def test_grid_lat_lon(): def test_elev_calc(): - pt1 = np.array((71.23, -34.70, 1234)) - pt2 = np.array((71.12, -34.16, 124)) - pt3 = np.array((71.29, -35.23, 30044)) + pt1 = np.array([[71.23, -34.70, 1234]]) + pt2 = np.array([[71.12, -34.16, 124]]) + pt3 = np.array([[71.29, -35.23, 30044]]) test_pts = np.vstack([pt1, pt2]) reverse_test_pts = np.vstack([pt2, pt1]) @@ -142,7 +142,7 @@ def test_uniformly_sample_triangle(): pt1 = np.array((71.23, -34.70, 1234)) pt2 = np.array((69.12, -38.16, 124)) pt3 = np.array((78.23, -36.23, 344)) - points = horizon.uniformly_sample_triangle(pt1, pt2, pt3, 5) + points = horizon._uniformly_sample_triangle(pt1, pt2, pt3, 5) p1 = tools.lle_to_xyz(pt1) p2 = tools.lle_to_xyz(pt2) @@ -263,7 +263,10 @@ def test_dni_horizon_adjustment(ephem_data): bad_horizon = np.ones(361) bad_expected = np.array([1, 1, 1, 1]) - bad_adjusted = horizon.dni_horizon_adjustment(bad_horizon, - ephem_data["zenith"], - ephem_data["azimuth"]) - assert_allclose(bad_adjusted, bad_expected, atol=1e-7) + try: + bad_adjusted = horizon.dni_horizon_adjustment(bad_horizon, + ephem_data["zenith"], + ephem_data["azimuth"]) + assert(False) + except ValueError: + pass diff --git a/pvlib/tools.py b/pvlib/tools.py index c68dfa0845..4ba26f724b 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -9,6 +9,10 @@ import pandas as pd import pytz +# Ellipsoid model of the Earth +earth_radius_equatorial = 6378137.0 +earth_radius_polar = 6356752.0 + def cosd(angle): """ @@ -433,10 +437,9 @@ def latitude_to_geocentric(phi): Latitude must be given in radians. [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf """ - # Equatorial Radius of the Earth (ellipsoid model) in meters - a = 6378137.0 - # Polar Radius of the Earth (ellipsoid model) in meters - b = 6356752.0 + + a = earth_radius_equatorial + b = earth_radius_polar return np.arctan(b**2/a**2*np.tan(phi)) @@ -446,10 +449,9 @@ def latitude_to_geodetic(phi): Latitude must be given in radians. [1] https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf """ - # Equatorial Radius of the Earth (ellipsoid model) in meters - a = 6378137.0 - # Polar Radius of the Earth (ellipsoid model) in meters - b = 6356752.0 + + a = earth_radius_equatorial + b = earth_radius_polar return np.arctan(a**2/b**2*np.tan(phi)) @@ -463,18 +465,16 @@ def lle_to_xyz(point): lon = np.atleast_1d(point.T[1]) elev = np.atleast_1d(point.T[2]) - # Equatorial Radius of the Earth (ellipsoid model) in meters - a = 6378137.0 - # Polar Radius of the Earth (ellipsoid model) in meters - b = 6356752.0 + a_sqrd = earth_radius_equatorial**2 + b_sqrd = earth_radius_polar**2 # convert to radians phi = np.radians(lat) theta = np.radians(lon) # compute radius of earth at each point - r = (a**2 * np.cos(phi))**2 + (b**2 * np.sin(phi))**2 - r = r / (a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2) + r = (a_sqrd * np.cos(phi))**2 + (b_sqrd * np.sin(phi))**2 + r = r / (a_sqrd * np.cos(phi)**2 + b_sqrd * np.sin(phi)**2) r = np.sqrt(r) h = r + elev @@ -493,10 +493,8 @@ def xyz_to_lle(point): The center of the earth is the origin in the xyz-space. The output latitude is assumed to be a common latitude (geodetic). """ - # Equatorial Radius of the Earth (ellipsoid model) in meters - a = 6378137.0 - # Polar Radius of the Earth (ellipsoid model) in meters - b = 6356752.0 + a_sqrd = earth_radius_equatorial**2 + b_sqrd = earth_radius_polar**2 x = np.atleast_1d(point.T[0]) y = np.atleast_1d(point.T[1]) @@ -504,7 +502,7 @@ def xyz_to_lle(point): # get corresponding point on earth's surface - t = np.sqrt((a*b)**2/(b**2*(x**2+y**2)+a**2*z**2)) + t = np.sqrt(a_sqrd*b_sqrd/(b_sqrd*(x**2+y**2)+a_sqrd*z**2)) t = t.reshape(point.shape[0], 1) point_s = t * point z_s = point_s.T[2] diff --git a/setup.py b/setup.py index de4125b588..f2e320ccbb 100755 --- a/setup.py +++ b/setup.py @@ -37,12 +37,10 @@ MAINTAINER_EMAIL = 'holmgren@email.arizona.edu' URL = 'https://github.com/pvlib/pvlib-python' -INSTALL_REQUIRES = ['numpy >= 1.13.3', +INSTALL_REQUIRES = ['numpy >= 1.10.4', 'pandas >= 0.18.1', 'pytz', - 'requests', - 'scipy', - 'matplotlib'] + 'requests'] TESTS_REQUIRE = ['nose', 'pytest', 'pytest-cov', 'pytest-mock', 'pytest-timeout'] EXTRAS_REQUIRE = { From 210fd1cf035ab29593c749d5d16ac145713959af Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 13 Aug 2019 15:04:59 -0700 Subject: [PATCH 28/29] docstring changes --- pvlib/horizon.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 18a4879541..05e3d54337 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -2,7 +2,6 @@ The ``horizon`` module contains functions for horizon profile modeling. There are various geometric utilities that are useful in horizon calculations. """ -from __future__ import division import itertools @@ -65,7 +64,7 @@ def elevation_and_azimuth(pt1, pt2): horizontal has a positive elevation angle. Also computes the azimuth defined as degrees East of North of the bearing from pt1 to pt2. This uses the Haversine formula. - The method used to calculate the elevation angle is discussed in [1]. + The trigonometry used to calculate the elevation angle is described in [1]. Parameters ---------- @@ -100,6 +99,9 @@ def elevation_and_azimuth(pt1, pt2): [36, 35, 231], [36, 35, 21]]) bearing, elev_angles = elevation_and_azimuth(site_loc, target_locs) + + + [1] https://aty.sdsu.edu/explain/atmos_refr/dip.html ''' # Equatorial Radius of the Earth (ellipsoid model) in meters a = 6378137.0 @@ -612,11 +614,15 @@ def calculate_dtf(horizon_azimuths, horizon_angles, integral is the cosine of the angle between the incoming radiation and the vector normal to the surface. The method calculates a sum of integrations from the "peak" of the sky dome down to the elevation angle of the horizon. + A similar method is used in [2] although it accounts for albedo and doesn't + account for the horizon. [1] Goss et al. (2014) Solar Energy 110, 410-419 - + [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 """ - assert(horizon_azimuths.shape[0] == horizon_angles.shape[0]) + if horizon_azimuths.shape[0] != horizon_angles.shape[0]: + raise ValueError('azimuths and elevation_angles must be of the same' + 'length.') tilt_rad = np.radians(surface_tilt) plane_az_rad = np.radians(surface_azimuth) a = np.sin(tilt_rad) * np.cos(plane_az_rad) From a630acced71786827146226ee606c3602a671f1c Mon Sep 17 00:00:00 2001 From: Joseph Palakapilly Date: Tue, 13 Aug 2019 15:08:37 -0700 Subject: [PATCH 29/29] reference update --- pvlib/horizon.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 05e3d54337..a5d1122cbf 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -614,10 +614,9 @@ def calculate_dtf(horizon_azimuths, horizon_angles, integral is the cosine of the angle between the incoming radiation and the vector normal to the surface. The method calculates a sum of integrations from the "peak" of the sky dome down to the elevation angle of the horizon. - A similar method is used in [2] although it accounts for albedo and doesn't - account for the horizon. + A similar method is used in section II of [1] although it is looking at + both ground and sky diffuse irradiation. - [1] Goss et al. (2014) Solar Energy 110, 410-419 [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 """ if horizon_azimuths.shape[0] != horizon_angles.shape[0]: