From 539bc5c051b8f57aba28c90c25e7ecdc27db0f25 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 26 Jan 2022 12:07:21 -0700 Subject: [PATCH 01/94] added horizon model first draft --- pvlib/horizon.py | 460 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 460 insertions(+) create mode 100644 pvlib/horizon.py diff --git a/pvlib/horizon.py b/pvlib/horizon.py new file mode 100644 index 0000000000..75157029ea --- /dev/null +++ b/pvlib/horizon.py @@ -0,0 +1,460 @@ +''' +The 'horizon' module contains function definitions that +retrive & calculate the surrounding horizon using DEM +elevation data +''' +import numpy as np +import gdal +from skimage.draw import line, line_aa +import pandas as pd +from scipy.linalg import polar +from scipy.signal import resample +from osgeo import gdal_array +import geoio + +def latlong(ds): + r''' + From a gdal dataset, retrive the geotransform + and return an latitude and longitude coordinates + for a DEM GeoTiff. + + Parameters + ---------- + ds : tuple + Geotransform parameters given by a osgeo.gdal.Dataset + + Returns + ------- + lat,long : tuple + Tuple of np.array of latitude,longitude corresponding to the pixels in the GeoTiff + + ''' + + width = ds.RasterXSize + height = ds.RasterYSize + + gt = ds.GetGeoTransform() + minx = gt[0] + miny = gt[3] + width*gt[4] + height*gt[5] + maxx = gt[0] + width*gt[1] + height*gt[2] + maxy = gt[3] + lat = np.linspace(miny,maxy,height) + long = np.linspace(minx,maxx,width) + return (lat,long) + +def load_DEM(filepath): + r''' + Loads a DEM from a .hgt file, segements into the GeoTiff + (elevation) and corresponding latitude and longitude + + Parameters + ---------- + filepath : string + Absolute or relative path to the .hgt file + + Returns + ------- + elevation : np.array + 2D numpy array, pixel value is elevation in meters + lat : np.array + latitudes corresponding to the pixels along the 0th dim of elevation + long : np.array + longitudes corresponding to the pixels along the 1th dim of elevation + + ''' + dataset = gdal.Open(filepath) + rasterArray = gdal_array.LoadFile(filepath) + elevation = rasterArray + lat, long = latlong(dataset) + return elevation, lat, long + +def get_pixel_coords(lat, lon, DEM_path): + r''' + + Gets pixel coordinates from the raster, given latitude and longitude + + Parameters + ---------- + lat : float + Latitude of the point + long: float + Longitude of the point + DEM_path: string + Path of the DEM .hgt file + + Returns + ------- + coords : tuple + Tuple with two elements, containing the x and y coordinate of the raster corresponding to latitiude and longitude + ''' + img = geoio.GeoImage(DEM_path) + return map(int, img.proj_to_raster(lon,lat)) + +def point_symmetry(x, y): + r""" + Reflect a point 8 ways to form a circle. + """ + return [( x, y), + ( y, x), + (-x, y), + (-y, x), + ( x, -y), + ( y, -x), + (-x, -y), + (-y, -x)] + +def bresenham_circ(r, center = (0,0)): + r""" + Use midpoint algorithum to build a rasterized circle. + Modified from https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm + to add circles not centered at the origin + + Parameters + ---------- + r : numeric + Radius of the circle + center: tuple + Center of the point (Cartesian) + + Returns + ------- + points : np.array + Array of shape (n,2) of points that form a rasterized circle. n depends on the radius of the circle. + """ + points = [] + x = 0 + y = -r + F_M = 1 - r + # Initial value for (0,-r) for 2x + 3 = 0x + 3 = 3. + d_e = 3 + # Initial value for (0,-r) for 2(x + y) + 5 = 0 - 2y + 5 = -2y + 5. + d_ne = -(r << 1) + 5 + points.extend(point_symmetry(x, y)) + while x < -y: + if F_M < 0: + F_M += d_e + else: + F_M += d_ne + d_ne += 2 + y += 1 + d_e += 2 + d_ne += 2 + x += 1 + points.extend(point_symmetry(x, y)) + points = np.array(points) + newx = points[:,0] + center[0] + newy = points[:,1] + center[1] + return np.vstack((newx,newy)) + +def pol2cart(r,theta): + r''' + Converts polar to cartesian + Parameters + ---------- + r : np.array + r values for the points + theta: np.array + theta values for the points + + Returns + ------- + x : np.array + x values for the points + y : np.array + y values for the points + + ''' + + z = r * np.exp(1j * theta) + x, y = z.real, z.imag + + return x, y + +def cart2pol(x, y): + r''' + Converts cartesian to polar + Parameters + ---------- + x : np.array + x values for the points + y : np.array + y values for the points + + Returns + ------- + r : np.array + r values for the points + theta: np.array + theta values for the points + + ''' + + z = x + y * 1j + r,theta = np.abs(z), np.angle(z) + + return r,theta + + +def sort_circ(pts, center = (0,0), az_len = 360): + r''' + Sort and resample points on a circle such that the zeroth element is + due east and the points move around the circle counter clockwise. + While in polar domain, resample points using FFT to obtain desired number of bins, typically 360, for degrees around the circle. + + Parameters + ---------- + pts : np.array + Array of shape (n, 2) of points in Cartesian system + center : tuple + Center (x,y) of the circle rasterized by pts + az_len : numeric + Desired number of points in the rasterized circle. + + Returns + ------- + points : np.array + Sorted and resampled points that comprise the circle + ''' + pts[0] -= center[0] + pts[1] -= center[1] + r, theta = cart2pol(pts[0], pts[1]) + stacked = np.vstack((theta,r)) + sort = np.sort(stacked, axis = 1) + sort = resample(sort, az_len, axis = 1) + theta = sort[0] + r = sort[1] + x,y = pol2cart(r,theta) + x += center[0] + y += center[1] + return np.vstack((x,y)).astype(int) + +def horizon_map(dem_pixel, dem_path, dem_res = 30.0, + view_distance=500, az_len = 36): + """ + Finds the horizon at point on a dem in pixel coordinates dem_pixel + + Parameters + ---------- + dem_pixel : tuple (int,int) + Point on the DEM expressed in pixel coordinates + dem_path : string + Path to a downloaded DEM .hgt + dem_res : float + Resolution of the DEM. The default is SRTM 30m + view_distance : int + Radius of the area of consideration. + az_len : int + Number of samples on the perimeter of the circle. + + Returns + ------- + azimuth : np.array + Numpy array of shape (az_len). Linearly spaced from [0,360] + elevation_angles: np.array + The calculated elevation values at each point of azimuth + profile: + The largest elevation in meters on the line between the observer and the highest point on the horizon within view_distance + + + + """ + # retrive data from the DEM + elevation, lat, long = load_DEM(dem_path) + + azimuth = np.linspace(0,360,az_len) + + # Use Midpoint Circle Algo/ Bresenham's circle + pts = bresenham_circ(view_distance,dem_pixel) + + # sort circle points and resample to desired number of points (az_len) + pts = sort_circ(pts, center = dem_pixel, az_len = az_len) + x0,y0 = dem_pixel + profile = np.zeros(azimuth.shape) + elevation_angles = np.zeros(azimuth.shape) + + for az in range(az_len): + + source_lon = pts[1][az] + source_lat = pts[0][az] + # check for outmapping + # for now, if out of the DEM, assign to the border + if source_lon >= y_bound: + source_lon = y_bound-1 + elif source_lon <= 0: + source_lon = 0 + if source_lat >= y_bound: + source_lat = y_bound-1 + elif source_lat <= 0: + source_lat = 0 + + target_lon = int(y0) + target_lat = int(x0) + + # draw a line from observer to horizon and record points + rr,cc = line(source_lon, source_lat, target_lon, target_lat) + + # get the highest elevation on the line + elvs = elevation[rr,cc] + elvs_on_line = elevation[rr,cc] + idx = np.stack((rr,cc), axis = -1) + highest_elv = np.max(elvs_on_line) + highest_point = idx[np.argmax(elvs_on_line)] + high_x, high_y = tuple(highest_point) + + # convert from altitude in m to elevation degrees. + xdist = np.abs(highest_point[0] - x0) + ydist = highest_elv - elevation[y0][x0] + elv_ang = np.arctan(ydist/xdist) + elevation_angles[az] = np.rad2deg(elv_ang) + profile[az] = highest_elv + + return azimuth, elevation_angles, profile + + +def calculate_dtf(horizon_azimuths, horizon_angles, + surface_tilt, surface_azimuth): + """ + Author: JPalakapillyKWH + 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 + ---------- + 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 + 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 + 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. + 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. + A similar method is used in section II of [1] although it is looking at + both ground and sky diffuse irradiation. + [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 + """ + 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) + b = np.sin(tilt_rad) * np.sin(plane_az_rad) + c = np.cos(tilt_rad) + + # this gets either a float or an array of zeros + dtf = np.multiply(0.0, surface_tilt) + num_points = horizon_azimuths.shape[0] + for i in range(horizon_azimuths.shape[0]): + az = np.radians(horizon_azimuths[i]) + 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 - 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 + +def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): + """ + Author: JPalakapillyKWH + 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. This comes from Equation 10 in [1] + 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 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 + -------- + 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. + [1] doi.org/10.1016/j.solener.2014.09.037 + """ + 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) + elevation_angle = 90.0 - declination + elevation_angle[mask] = 0.0 + + return elevation_angle + + +def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): + ''' + Author: JPalakapillyKWH + 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 + ---------- + 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] != 360): + 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 + horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] + mask = solar_zenith > horizon_zenith + adjustment[mask] = 0 + return adjustment \ No newline at end of file From 85754b5566303a5400b9a7e123eb4d7ccd042447 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 26 Jan 2022 13:05:10 -0700 Subject: [PATCH 02/94] added get horizon function --- pvlib/iotools/pvgis.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 64b23139eb..ba957f9e0f 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -659,3 +659,34 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta + + +def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): + r''' + Get horizon data from PVGIS + Parameters + ---------- + latitude : float + Latitude in degrees north + longitude : float + Longitude in degrees east + proxies: dict + Dictionary of proxies to access through a corporate network + url: string + Base URL for PVGIS + + Returns + ------- + df : pd.DataFrame + Pandas dataframe of the retrived horizon + ''' + res = requests.get( url + 'printhorizon?lat=45&lon=8', proxies = proxies, verify = False) + res.raise_for_status() + string = str(io.BytesIO(res.content).read().decode('UTF-8')) + # the horizon data is given in a different format then the others + # numpy has an easier time decoding it + array = np.genfromtxt(io.StringIO(string), skip_header = 4, skip_footer = 7) + df = pd.DataFrame(array) + # Get the column names + df.columns = [s for s in string.splitlines()[3].split('\t') if s] + return df \ No newline at end of file From 4985e9cd827ab000e3e9148cfb62f54877483500 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 31 Jan 2022 10:01:23 -0700 Subject: [PATCH 03/94] fixing stickler items & updates --- pvlib/horizon.py | 111 ++++++++++++++++++++--------------------- pvlib/iotools/pvgis.py | 3 +- 2 files changed, 54 insertions(+), 60 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 75157029ea..1a94c8106b 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,20 +1,17 @@ ''' -The 'horizon' module contains function definitions that -retrive & calculate the surrounding horizon using DEM +The 'horizon' module contains function definitions that +retrive & calculate the surrounding horizon using DEM elevation data ''' import numpy as np import gdal -from skimage.draw import line, line_aa -import pandas as pd -from scipy.linalg import polar +from skimage.draw import line from scipy.signal import resample from osgeo import gdal_array -import geoio - +import geoio +import numpy as np def latlong(ds): - r''' - From a gdal dataset, retrive the geotransform + r'''From a gdal dataset, retrive the geotransform and return an latitude and longitude coordinates for a DEM GeoTiff. @@ -42,9 +39,9 @@ def latlong(ds): long = np.linspace(minx,maxx,width) return (lat,long) + def load_DEM(filepath): - r''' - Loads a DEM from a .hgt file, segements into the GeoTiff + r'''Loads a DEM from a .hgt file, segements into the GeoTiff (elevation) and corresponding latitude and longitude Parameters @@ -68,10 +65,9 @@ def load_DEM(filepath): lat, long = latlong(dataset) return elevation, lat, long + def get_pixel_coords(lat, lon, DEM_path): - r''' - - Gets pixel coordinates from the raster, given latitude and longitude + r'''Gets pixel coordinates from the raster, given latitude and longitude Parameters ---------- @@ -85,12 +81,14 @@ def get_pixel_coords(lat, lon, DEM_path): Returns ------- coords : tuple - Tuple with two elements, containing the x and y coordinate of the raster corresponding to latitiude and longitude + Tuple with two elements, containing the x and y coordinate + of the raster corresponding to latitiude and longitude ''' img = geoio.GeoImage(DEM_path) return map(int, img.proj_to_raster(lon,lat)) -def point_symmetry(x, y): + +def _point_symmetry(x, y): r""" Reflect a point 8 ways to form a circle. """ @@ -103,10 +101,11 @@ def point_symmetry(x, y): (-x, -y), (-y, -x)] -def bresenham_circ(r, center = (0,0)): - r""" - Use midpoint algorithum to build a rasterized circle. - Modified from https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm + +def _bresenham_circ(r, center = (0,0)): + r""" Use midpoint algorithum to build a rasterized circle. + Modified from + https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm to add circles not centered at the origin Parameters @@ -119,7 +118,8 @@ def bresenham_circ(r, center = (0,0)): Returns ------- points : np.array - Array of shape (n,2) of points that form a rasterized circle. n depends on the radius of the circle. + Array of shape (n,2) of points that form a rasterized circle. + n depends on the radius of the circle. """ points = [] x = 0 @@ -129,7 +129,7 @@ def bresenham_circ(r, center = (0,0)): d_e = 3 # Initial value for (0,-r) for 2(x + y) + 5 = 0 - 2y + 5 = -2y + 5. d_ne = -(r << 1) + 5 - points.extend(point_symmetry(x, y)) + points.extend(_point_symmetry(x, y)) while x < -y: if F_M < 0: F_M += d_e @@ -140,15 +140,15 @@ def bresenham_circ(r, center = (0,0)): d_e += 2 d_ne += 2 x += 1 - points.extend(point_symmetry(x, y)) + points.extend(_point_symmetry(x, y)) points = np.array(points) newx = points[:,0] + center[0] newy = points[:,1] + center[1] return np.vstack((newx,newy)) -def pol2cart(r,theta): - r''' - Converts polar to cartesian + +def _pol2cart(r,theta): + r'''Converts polar to cartesian Parameters ---------- r : np.array @@ -167,12 +167,11 @@ def pol2cart(r,theta): z = r * np.exp(1j * theta) x, y = z.real, z.imag - return x, y -def cart2pol(x, y): - r''' - Converts cartesian to polar + +def _cart2pol(x, y): + r'''Converts cartesian to polar Parameters ---------- x : np.array @@ -190,16 +189,15 @@ def cart2pol(x, y): ''' z = x + y * 1j - r,theta = np.abs(z), np.angle(z) + r, theta = np.abs(z), np.angle(z) + return r, theta - return r,theta - -def sort_circ(pts, center = (0,0), az_len = 360): - r''' - Sort and resample points on a circle such that the zeroth element is +def _sort_circ(pts, center = (0,0), az_len = 360): + r'''Sort and resample points on a circle such that the zeroth element is due east and the points move around the circle counter clockwise. - While in polar domain, resample points using FFT to obtain desired number of bins, typically 360, for degrees around the circle. + While in polar domain, resample points using FFT to obtain desired number of bins, + typically 360, for degrees around the circle. Parameters ---------- @@ -217,21 +215,21 @@ def sort_circ(pts, center = (0,0), az_len = 360): ''' pts[0] -= center[0] pts[1] -= center[1] - r, theta = cart2pol(pts[0], pts[1]) - stacked = np.vstack((theta,r)) + r, theta = _cart2pol(pts[0], pts[1]) + stacked = np.vstack((theta, r)) sort = np.sort(stacked, axis = 1) sort = resample(sort, az_len, axis = 1) theta = sort[0] r = sort[1] - x,y = pol2cart(r,theta) + x, y = _pol2cart(r,theta) x += center[0] y += center[1] - return np.vstack((x,y)).astype(int) + return np.vstack((x, y)).astype(int) + def horizon_map(dem_pixel, dem_path, dem_res = 30.0, view_distance=500, az_len = 36): - """ - Finds the horizon at point on a dem in pixel coordinates dem_pixel + r"""Finds the horizon at point on a dem in pixel coordinates dem_pixel Parameters ---------- @@ -253,9 +251,8 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, elevation_angles: np.array The calculated elevation values at each point of azimuth profile: - The largest elevation in meters on the line between the observer and the highest point on the horizon within view_distance - - + The largest elevation in meters on the line between the observer + and the highest point on the horizon within view_distance """ # retrive data from the DEM @@ -264,14 +261,14 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, azimuth = np.linspace(0,360,az_len) # Use Midpoint Circle Algo/ Bresenham's circle - pts = bresenham_circ(view_distance,dem_pixel) + pts = _bresenham_circ(view_distance,dem_pixel) # sort circle points and resample to desired number of points (az_len) - pts = sort_circ(pts, center = dem_pixel, az_len = az_len) + pts = _sort_circ(pts, center = dem_pixel, az_len = az_len) x0,y0 = dem_pixel profile = np.zeros(azimuth.shape) elevation_angles = np.zeros(azimuth.shape) - + x_bound, y_bound = elevation.shape for az in range(az_len): source_lon = pts[1][az] @@ -291,12 +288,12 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, target_lat = int(x0) # draw a line from observer to horizon and record points - rr,cc = line(source_lon, source_lat, target_lon, target_lat) + rr, cc = line(source_lon, source_lat, target_lon, target_lat) # get the highest elevation on the line - elvs = elevation[rr,cc] + elvs = elevation[rr, cc] elvs_on_line = elevation[rr,cc] - idx = np.stack((rr,cc), axis = -1) + idx = np.stack((rr, cc), axis = -1) highest_elv = np.max(elvs_on_line) highest_point = idx[np.argmax(elvs_on_line)] high_x, high_y = tuple(highest_point) @@ -313,8 +310,7 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, def calculate_dtf(horizon_azimuths, horizon_angles, surface_tilt, surface_azimuth): - """ - Author: JPalakapillyKWH + r"""Author: JPalakapillyKWH 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]. @@ -381,9 +377,9 @@ def calculate_dtf(horizon_azimuths, horizon_angles, dtf += 2 * (first_term + second_term) / num_points return dtf + def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): - """ - Author: JPalakapillyKWH + r"""Author: JPalakapillyKWH 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. This comes from Equation 10 in [1] @@ -421,8 +417,7 @@ def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): - ''' - Author: JPalakapillyKWH + r'''Author: JPalakapillyKWH 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, diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index ba957f9e0f..841726c74e 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -19,6 +19,7 @@ from pathlib import Path import requests import pandas as pd +import numpy as np from pvlib.iotools import read_epw, parse_epw import warnings from pvlib._deprecation import pvlibDeprecationWarning @@ -659,8 +660,6 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta - - def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): r''' Get horizon data from PVGIS From 132ab4e2b200f9dcbc51152b439c27e6f09f3cbf Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 31 Jan 2022 10:07:26 -0700 Subject: [PATCH 04/94] stickler whitespace --- pvlib/horizon.py | 68 ++++++++++++++++++++---------------------- pvlib/iotools/pvgis.py | 3 +- 2 files changed, 34 insertions(+), 37 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 1a94c8106b..7fe7091ebc 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -9,17 +9,18 @@ from scipy.signal import resample from osgeo import gdal_array import geoio -import numpy as np + + def latlong(ds): r'''From a gdal dataset, retrive the geotransform and return an latitude and longitude coordinates for a DEM GeoTiff. - + Parameters ---------- ds : tuple Geotransform parameters given by a osgeo.gdal.Dataset - + Returns ------- lat,long : tuple @@ -29,10 +30,9 @@ def latlong(ds): width = ds.RasterXSize height = ds.RasterYSize - gt = ds.GetGeoTransform() minx = gt[0] - miny = gt[3] + width*gt[4] + height*gt[5] + miny = gt[3] + width*gt[4] + height*gt[5] maxx = gt[0] + width*gt[1] + height*gt[2] maxy = gt[3] lat = np.linspace(miny,maxy,height) @@ -43,12 +43,11 @@ def latlong(ds): def load_DEM(filepath): r'''Loads a DEM from a .hgt file, segements into the GeoTiff (elevation) and corresponding latitude and longitude - + Parameters ---------- filepath : string Absolute or relative path to the .hgt file - Returns ------- elevation : np.array @@ -57,9 +56,8 @@ def load_DEM(filepath): latitudes corresponding to the pixels along the 0th dim of elevation long : np.array longitudes corresponding to the pixels along the 1th dim of elevation - ''' - dataset = gdal.Open(filepath) + dataset = gdal.Open(filepath) rasterArray = gdal_array.LoadFile(filepath) elevation = rasterArray lat, long = latlong(dataset) @@ -68,7 +66,7 @@ def load_DEM(filepath): def get_pixel_coords(lat, lon, DEM_path): r'''Gets pixel coordinates from the raster, given latitude and longitude - + Parameters ---------- lat : float @@ -77,7 +75,7 @@ def get_pixel_coords(lat, lon, DEM_path): Longitude of the point DEM_path: string Path of the DEM .hgt file - + Returns ------- coords : tuple @@ -89,8 +87,8 @@ def get_pixel_coords(lat, lon, DEM_path): def _point_symmetry(x, y): - r""" - Reflect a point 8 ways to form a circle. + r""" + Reflect a point 8 ways to form a circle. """ return [( x, y), ( y, x), @@ -104,22 +102,22 @@ def _point_symmetry(x, y): def _bresenham_circ(r, center = (0,0)): r""" Use midpoint algorithum to build a rasterized circle. - Modified from + Modified from https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm to add circles not centered at the origin - + Parameters ---------- r : numeric Radius of the circle center: tuple Center of the point (Cartesian) - + Returns ------- points : np.array Array of shape (n,2) of points that form a rasterized circle. - n depends on the radius of the circle. + n depends on the radius of the circle. """ points = [] x = 0 @@ -155,7 +153,7 @@ def _pol2cart(r,theta): r values for the points theta: np.array theta values for the points - + Returns ------- x : np.array @@ -198,7 +196,7 @@ def _sort_circ(pts, center = (0,0), az_len = 360): due east and the points move around the circle counter clockwise. While in polar domain, resample points using FFT to obtain desired number of bins, typically 360, for degrees around the circle. - + Parameters ---------- pts : np.array @@ -206,7 +204,7 @@ def _sort_circ(pts, center = (0,0), az_len = 360): center : tuple Center (x,y) of the circle rasterized by pts az_len : numeric - Desired number of points in the rasterized circle. + Desired number of points in the rasterized circle. Returns ------- @@ -230,7 +228,7 @@ def _sort_circ(pts, center = (0,0), az_len = 360): def horizon_map(dem_pixel, dem_path, dem_res = 30.0, view_distance=500, az_len = 36): r"""Finds the horizon at point on a dem in pixel coordinates dem_pixel - + Parameters ---------- dem_pixel : tuple (int,int) @@ -240,10 +238,10 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, dem_res : float Resolution of the DEM. The default is SRTM 30m view_distance : int - Radius of the area of consideration. + Radius of the area of consideration. az_len : int Number of samples on the perimeter of the circle. - + Returns ------- azimuth : np.array @@ -251,18 +249,18 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, elevation_angles: np.array The calculated elevation values at each point of azimuth profile: - The largest elevation in meters on the line between the observer + The largest elevation in meters on the line between the observer and the highest point on the horizon within view_distance - + """ # retrive data from the DEM elevation, lat, long = load_DEM(dem_path) - + azimuth = np.linspace(0,360,az_len) - + # Use Midpoint Circle Algo/ Bresenham's circle pts = _bresenham_circ(view_distance,dem_pixel) - + # sort circle points and resample to desired number of points (az_len) pts = _sort_circ(pts, center = dem_pixel, az_len = az_len) x0,y0 = dem_pixel @@ -270,7 +268,6 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, elevation_angles = np.zeros(azimuth.shape) x_bound, y_bound = elevation.shape for az in range(az_len): - source_lon = pts[1][az] source_lat = pts[0][az] # check for outmapping @@ -283,13 +280,13 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, source_lat = y_bound-1 elif source_lat <= 0: source_lat = 0 - + target_lon = int(y0) target_lat = int(x0) - + # draw a line from observer to horizon and record points rr, cc = line(source_lon, source_lat, target_lon, target_lat) - + # get the highest elevation on the line elvs = elevation[rr, cc] elvs_on_line = elevation[rr,cc] @@ -297,14 +294,13 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, highest_elv = np.max(elvs_on_line) highest_point = idx[np.argmax(elvs_on_line)] high_x, high_y = tuple(highest_point) - - # convert from altitude in m to elevation degrees. + + # convert from altitude in m to elevation degrees. xdist = np.abs(highest_point[0] - x0) - ydist = highest_elv - elevation[y0][x0] + ydist = highest_elv - elevation[y0][x0] elv_ang = np.arctan(ydist/xdist) elevation_angles[az] = np.rad2deg(elv_ang) profile[az] = highest_elv - return azimuth, elevation_angles, profile diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 841726c74e..32de347c85 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -660,6 +660,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta + def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): r''' Get horizon data from PVGIS @@ -673,7 +674,7 @@ def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): Dictionary of proxies to access through a corporate network url: string Base URL for PVGIS - + Returns ------- df : pd.DataFrame From 64bbc975d658a2cc8739f6d4cfab1c31cecc0fa8 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 31 Jan 2022 10:18:05 -0700 Subject: [PATCH 05/94] finishing up stickler --- pvlib/horizon.py | 81 ++++++++++++++++++++++++------------------ pvlib/iotools/pvgis.py | 10 +++--- 2 files changed, 52 insertions(+), 39 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 7fe7091ebc..1e7a6e38b1 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -24,7 +24,8 @@ def latlong(ds): Returns ------- lat,long : tuple - Tuple of np.array of latitude,longitude corresponding to the pixels in the GeoTiff + Tuple of np.array of latitude,longitude + corresponding to the pixels in the GeoTiff ''' @@ -35,9 +36,9 @@ def latlong(ds): miny = gt[3] + width*gt[4] + height*gt[5] maxx = gt[0] + width*gt[1] + height*gt[2] maxy = gt[3] - lat = np.linspace(miny,maxy,height) - long = np.linspace(minx,maxx,width) - return (lat,long) + lat = np.linspace(miny, maxy, height) + long = np.linspace(minx, maxx, width) + return (lat, long) def load_DEM(filepath): @@ -83,20 +84,30 @@ def get_pixel_coords(lat, lon, DEM_path): of the raster corresponding to latitiude and longitude ''' img = geoio.GeoImage(DEM_path) - return map(int, img.proj_to_raster(lon,lat)) + return map(int, img.proj_to_raster(lon, lat)) def _point_symmetry(x, y): - r""" - Reflect a point 8 ways to form a circle. + r"""Reflect a point 8 ways to form a circle. + Parameters + ---------- + x : numeric + x point of the circle + y: numeric + y point of the circle + + Returns + ------- + points : list + List of reflected points """ - return [( x, y), - ( y, x), - (-x, y), - (-y, x), - ( x, -y), - ( y, -x), - (-x, -y), + return [(x, y), + (y, x), + (-x, y), + (-y, x), + (x, -y), + (y, -x), + (-x, -y), (-y, -x)] @@ -140,12 +151,12 @@ def _bresenham_circ(r, center = (0,0)): x += 1 points.extend(_point_symmetry(x, y)) points = np.array(points) - newx = points[:,0] + center[0] - newy = points[:,1] + center[1] - return np.vstack((newx,newy)) + newx = points[:, 0] + center[0] + newy = points[:, 1] + center[1] + return np.vstack((newx, newy)) -def _pol2cart(r,theta): +def _pol2cart(r, theta): r'''Converts polar to cartesian Parameters ---------- @@ -191,11 +202,12 @@ def _cart2pol(x, y): return r, theta -def _sort_circ(pts, center = (0,0), az_len = 360): +def _sort_circ(pts, center=(0,0), az_len=360): r'''Sort and resample points on a circle such that the zeroth element is due east and the points move around the circle counter clockwise. - While in polar domain, resample points using FFT to obtain desired number of bins, - typically 360, for degrees around the circle. + While in polar domain, resample points using FFT to + obtain desired number of bins,typically 360, for + degrees around the circle. Parameters ---------- @@ -215,18 +227,18 @@ def _sort_circ(pts, center = (0,0), az_len = 360): pts[1] -= center[1] r, theta = _cart2pol(pts[0], pts[1]) stacked = np.vstack((theta, r)) - sort = np.sort(stacked, axis = 1) - sort = resample(sort, az_len, axis = 1) + sort = np.sort(stacked, axis=1) + sort = resample(sort, az_len, axis=1) theta = sort[0] r = sort[1] - x, y = _pol2cart(r,theta) + x, y = _pol2cart(r, theta) x += center[0] y += center[1] return np.vstack((x, y)).astype(int) -def horizon_map(dem_pixel, dem_path, dem_res = 30.0, - view_distance=500, az_len = 36): +def horizon_map(dem_pixel, dem_path, dem_res=30.0, + view_distance=500, az_len=36): r"""Finds the horizon at point on a dem in pixel coordinates dem_pixel Parameters @@ -256,14 +268,14 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, # retrive data from the DEM elevation, lat, long = load_DEM(dem_path) - azimuth = np.linspace(0,360,az_len) + azimuth = np.linspace(0, 360, az_len) # Use Midpoint Circle Algo/ Bresenham's circle - pts = _bresenham_circ(view_distance,dem_pixel) + pts = _bresenham_circ(view_distance, dem_pixel) # sort circle points and resample to desired number of points (az_len) - pts = _sort_circ(pts, center = dem_pixel, az_len = az_len) - x0,y0 = dem_pixel + pts = _sort_circ(pts, center=dem_pixel, az_len=az_len) + x0, y0 = dem_pixel profile = np.zeros(azimuth.shape) elevation_angles = np.zeros(azimuth.shape) x_bound, y_bound = elevation.shape @@ -288,15 +300,14 @@ def horizon_map(dem_pixel, dem_path, dem_res = 30.0, rr, cc = line(source_lon, source_lat, target_lon, target_lat) # get the highest elevation on the line - elvs = elevation[rr, cc] - elvs_on_line = elevation[rr,cc] - idx = np.stack((rr, cc), axis = -1) + elvs_on_line = elevation[rr, cc] + idx = np.stack((rr, cc), axis=-1) highest_elv = np.max(elvs_on_line) highest_point = idx[np.argmax(elvs_on_line)] high_x, high_y = tuple(highest_point) # convert from altitude in m to elevation degrees. - xdist = np.abs(highest_point[0] - x0) + xdist = np.abs(highest_point[0]-x0) ydist = highest_elv - elevation[y0][x0] elv_ang = np.arctan(ydist/xdist) elevation_angles[az] = np.rad2deg(elv_ang) @@ -448,4 +459,4 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] mask = solar_zenith > horizon_zenith adjustment[mask] = 0 - return adjustment \ No newline at end of file + return adjustment diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 32de347c85..3f922e29c6 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -661,7 +661,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta -def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): +def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): r''' Get horizon data from PVGIS Parameters @@ -680,13 +680,15 @@ def get_pvgis_horizon(latitude, longitude, proxies = None, url = URL): df : pd.DataFrame Pandas dataframe of the retrived horizon ''' - res = requests.get( url + 'printhorizon?lat=45&lon=8', proxies = proxies, verify = False) + res = requests.get(url + 'printhorizon?lat=45&lon=8', + proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) # the horizon data is given in a different format then the others # numpy has an easier time decoding it - array = np.genfromtxt(io.StringIO(string), skip_header = 4, skip_footer = 7) + array = np.genfromtxt(io.StringIO(string), + skip_header=4, skip_footer=7) df = pd.DataFrame(array) # Get the column names df.columns = [s for s in string.splitlines()[3].split('\t') if s] - return df \ No newline at end of file + return df From 19e1073b51dde09c4f2dad2fc8ac98763ab77d37 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 31 Jan 2022 10:20:49 -0700 Subject: [PATCH 06/94] finishing up stickler, 2 --- pvlib/horizon.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 1e7a6e38b1..c010fd4503 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -101,17 +101,17 @@ def _point_symmetry(x, y): points : list List of reflected points """ - return [(x, y), - (y, x), - (-x, y), - (-y, x), - (x, -y), - (y, -x), - (-x, -y), + return [(x, y), + (y, x), + (-x, y), + (-y, x), + (x, -y), + (y, -x), + (-x, -y), (-y, -x)] -def _bresenham_circ(r, center = (0,0)): +def _bresenham_circ(r, center=(0, 0)): r""" Use midpoint algorithum to build a rasterized circle. Modified from https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm @@ -202,11 +202,11 @@ def _cart2pol(x, y): return r, theta -def _sort_circ(pts, center=(0,0), az_len=360): +def _sort_circ(pts, center=(0, 0), az_len=360): r'''Sort and resample points on a circle such that the zeroth element is due east and the points move around the circle counter clockwise. While in polar domain, resample points using FFT to - obtain desired number of bins,typically 360, for + obtain desired number of bins,typically 360, for degrees around the circle. Parameters @@ -275,7 +275,7 @@ def horizon_map(dem_pixel, dem_path, dem_res=30.0, # sort circle points and resample to desired number of points (az_len) pts = _sort_circ(pts, center=dem_pixel, az_len=az_len) - x0, y0 = dem_pixel + x0, y0 = dem_pixel profile = np.zeros(azimuth.shape) elevation_angles = np.zeros(azimuth.shape) x_bound, y_bound = elevation.shape From d563a81a85a1e4debb14d4b1e88e4ab81027ef36 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 1 Feb 2022 10:56:09 -0700 Subject: [PATCH 07/94] added optional deps, changed argument of horizon_map from path to array --- pvlib/horizon.py | 20 ++++++++------------ setup.py | 2 +- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index c010fd4503..238b36d063 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -4,12 +4,7 @@ elevation data ''' import numpy as np -import gdal -from skimage.draw import line from scipy.signal import resample -from osgeo import gdal_array -import geoio - def latlong(ds): r'''From a gdal dataset, retrive the geotransform @@ -58,6 +53,8 @@ def load_DEM(filepath): long : np.array longitudes corresponding to the pixels along the 1th dim of elevation ''' + import gdal + from osgeo import gdal_array dataset = gdal.Open(filepath) rasterArray = gdal_array.LoadFile(filepath) elevation = rasterArray @@ -83,7 +80,8 @@ def get_pixel_coords(lat, lon, DEM_path): Tuple with two elements, containing the x and y coordinate of the raster corresponding to latitiude and longitude ''' - img = geoio.GeoImage(DEM_path) + from geoio import GeoImage + img = GeoImage(DEM_path) return map(int, img.proj_to_raster(lon, lat)) @@ -237,7 +235,7 @@ def _sort_circ(pts, center=(0, 0), az_len=360): return np.vstack((x, y)).astype(int) -def horizon_map(dem_pixel, dem_path, dem_res=30.0, +def horizon_map(dem_pixel, elevation, dem_res=30.0, view_distance=500, az_len=36): r"""Finds the horizon at point on a dem in pixel coordinates dem_pixel @@ -245,8 +243,8 @@ def horizon_map(dem_pixel, dem_path, dem_res=30.0, ---------- dem_pixel : tuple (int,int) Point on the DEM expressed in pixel coordinates - dem_path : string - Path to a downloaded DEM .hgt + elevation : np.array + nxn DEM of elevation values dem_res : float Resolution of the DEM. The default is SRTM 30m view_distance : int @@ -265,9 +263,7 @@ def horizon_map(dem_pixel, dem_path, dem_res=30.0, and the highest point on the horizon within view_distance """ - # retrive data from the DEM - elevation, lat, long = load_DEM(dem_path) - + from skimage.draw import line azimuth = np.linspace(0, 360, az_len) # Use Midpoint Circle Algo/ Bresenham's circle diff --git a/setup.py b/setup.py index b2cbf486d5..d654aa4705 100755 --- a/setup.py +++ b/setup.py @@ -55,7 +55,7 @@ EXTRAS_REQUIRE = { 'optional': ['cython', 'ephem', 'netcdf4', 'nrel-pysam', 'numba', 'pvfactors', 'siphon', 'statsmodels', - 'cftime >= 1.1.1'], + 'cftime >= 1.1.1', 'geoio', 'gdal', 'osgeo', 'scikit-image'], 'doc': ['ipython', 'matplotlib', 'sphinx == 3.1.2', 'pydata-sphinx-theme', 'sphinx-gallery', 'docutils == 0.15.2', 'pillow', 'netcdf4', 'siphon', From e9856da4ae0956780b1beff9dabfc0ff9500c64d Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 1 Feb 2022 10:57:26 -0700 Subject: [PATCH 08/94] stickler --- pvlib/horizon.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 238b36d063..2c5c248bec 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -6,6 +6,7 @@ import numpy as np from scipy.signal import resample + def latlong(ds): r'''From a gdal dataset, retrive the geotransform and return an latitude and longitude coordinates From c392e29b1c6b247b3bdb0105e894433b9ee950c7 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:27:53 -0700 Subject: [PATCH 09/94] Update pvlib/iotools/pvgis.py Co-authored-by: Karel De Brabandere --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 3f922e29c6..d31cb87ccf 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -680,7 +680,7 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): df : pd.DataFrame Pandas dataframe of the retrived horizon ''' - res = requests.get(url + 'printhorizon?lat=45&lon=8', + res = requests.get(url +f 'printhorizon?lat={latitude}&lon={longitude}', proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) From e8d6100e8cd71a49aa28105a9b4d6690920656d0 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:29:19 -0700 Subject: [PATCH 10/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 2c5c248bec..25ae59c9f3 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -38,7 +38,7 @@ def latlong(ds): def load_DEM(filepath): - r'''Loads a DEM from a .hgt file, segements into the GeoTiff + r'''Loads a DEM from a .hgt file, segments into the GeoTiff (elevation) and corresponding latitude and longitude Parameters From 3dca8bb5d9831ac50ad0a66c9c6f8a2a17ea63b1 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:29:29 -0700 Subject: [PATCH 11/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 25ae59c9f3..4e01fae076 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -45,6 +45,7 @@ def load_DEM(filepath): ---------- filepath : string Absolute or relative path to the .hgt file + Returns ------- elevation : np.array From 93d914d5f52a6fd3bc11338bff705c5de2b2450d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:30:17 -0700 Subject: [PATCH 12/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 4e01fae076..d03071c76c 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -89,6 +89,7 @@ def get_pixel_coords(lat, lon, DEM_path): def _point_symmetry(x, y): r"""Reflect a point 8 ways to form a circle. + Parameters ---------- x : numeric From 5c1bc2dae835da5ebb1d7a6a81b2990aedf660d8 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:30:45 -0700 Subject: [PATCH 13/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index d03071c76c..f43337df3a 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -164,7 +164,7 @@ def _pol2cart(r, theta): r : np.array r values for the points theta: np.array - theta values for the points + theta values for the points. [radian] Returns ------- From eb8f8c54309b4759159440fc28155aeaa033d1ff Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:30:55 -0700 Subject: [PATCH 14/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index f43337df3a..05021edd7e 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -194,7 +194,7 @@ def _cart2pol(x, y): r : np.array r values for the points theta: np.array - theta values for the points + theta values for the points. [radian] ''' From 5a3c9cc8d16b755b14030febff9f8870c33ebd7b Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:31:24 -0700 Subject: [PATCH 15/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 05021edd7e..a78236b095 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -8,8 +8,8 @@ def latlong(ds): - r'''From a gdal dataset, retrive the geotransform - and return an latitude and longitude coordinates + r'''From a gdal dataset, retrieve the geotransform + and return latitude and longitude coordinates for a DEM GeoTiff. Parameters From 515ed85fb1cb58cf9bf723cb309b6aad8f167308 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:33:10 -0700 Subject: [PATCH 16/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index a78236b095..a73f295699 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -252,7 +252,7 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, Resolution of the DEM. The default is SRTM 30m view_distance : int Radius of the area of consideration. - az_len : int + az_len : int, default 360 Number of samples on the perimeter of the circle. Returns From 1c3a3988ed343525ff675483213e9d2656912bbf Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:33:24 -0700 Subject: [PATCH 17/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index a73f295699..e70bc2d4a5 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -207,7 +207,7 @@ def _sort_circ(pts, center=(0, 0), az_len=360): r'''Sort and resample points on a circle such that the zeroth element is due east and the points move around the circle counter clockwise. While in polar domain, resample points using FFT to - obtain desired number of bins,typically 360, for + obtain desired number of bins, typically 360, for degrees around the circle. Parameters From a191d2df6213a530428b9c80d4a93c64c04fb7a6 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:33:35 -0700 Subject: [PATCH 18/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index e70bc2d4a5..5522842afa 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -240,7 +240,7 @@ def _sort_circ(pts, center=(0, 0), az_len=360): def horizon_map(dem_pixel, elevation, dem_res=30.0, view_distance=500, az_len=36): - r"""Finds the horizon at point on a dem in pixel coordinates dem_pixel + r"""Finds the horizon for the point ``dem_pixel``. Parameters ---------- From 82ba2a529fa80c7279615fa705ce1a94d1236ab0 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:39:07 -0700 Subject: [PATCH 19/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 5522842afa..89d7830000 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -216,7 +216,7 @@ def _sort_circ(pts, center=(0, 0), az_len=360): Array of shape (n, 2) of points in Cartesian system center : tuple Center (x,y) of the circle rasterized by pts - az_len : numeric + az_len : int, default 360 Desired number of points in the rasterized circle. Returns From ca2eec3ebc649466e4934f2b022f88f81ca9f6a2 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:39:18 -0700 Subject: [PATCH 20/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 89d7830000..bacf889138 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -308,7 +308,7 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, # convert from altitude in m to elevation degrees. xdist = np.abs(highest_point[0]-x0) ydist = highest_elv - elevation[y0][x0] - elv_ang = np.arctan(ydist/xdist) + elv_ang = np.arctan2(ydist, xdist) elevation_angles[az] = np.rad2deg(elv_ang) profile[az] = highest_elv return azimuth, elevation_angles, profile From 996b3f90b08ae70446aee31a702b4c14c1867183 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:39:44 -0700 Subject: [PATCH 21/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index bacf889138..d24ad14add 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -248,8 +248,8 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, Point on the DEM expressed in pixel coordinates elevation : np.array nxn DEM of elevation values - dem_res : float - Resolution of the DEM. The default is SRTM 30m + dem_res : float, default 30 + Resolution of the DEM. The default is SRTM's 30m. [m] view_distance : int Radius of the area of consideration. az_len : int, default 360 From 371718614c37e8ab49eece897a61815d084e9173 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:40:19 -0700 Subject: [PATCH 22/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index d24ad14add..1800090e09 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -120,7 +120,7 @@ def _bresenham_circ(r, center=(0, 0)): Parameters ---------- - r : numeric + r : float Radius of the circle center: tuple Center of the point (Cartesian) From f9a2193b9ac928ce387aae408df7d2f0479da5e1 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 7 Feb 2022 09:40:35 -0700 Subject: [PATCH 23/94] Update pvlib/horizon.py Co-authored-by: Cliff Hansen --- pvlib/horizon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 1800090e09..901b77ac2f 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -250,8 +250,8 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, nxn DEM of elevation values dem_res : float, default 30 Resolution of the DEM. The default is SRTM's 30m. [m] - view_distance : int - Radius of the area of consideration. + view_distance : float, default 500 + Radius of the area of consideration. [m] az_len : int, default 360 Number of samples on the perimeter of the circle. From 989b700b2753dd7f9b3cb31f9524771d1a950ee2 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 09:56:22 -0700 Subject: [PATCH 24/94] updates --- pvlib/horizon.py | 167 +++-------------------------------------- pvlib/iotools/pvgis.py | 2 +- pvlib/shading.py | 158 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 171 insertions(+), 156 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 901b77ac2f..4a0e312cd8 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,7 +1,8 @@ ''' The 'horizon' module contains function definitions that retrive & calculate the surrounding horizon using DEM -elevation data +elevation data. +Optional dependencies for this module include gdal, osgeo, geoio, and scikit-image. ''' import numpy as np from scipy.signal import resample @@ -23,8 +24,11 @@ def latlong(ds): Tuple of np.array of latitude,longitude corresponding to the pixels in the GeoTiff + Notes + ------ + Latitude and longitude are in decimal degrees and negative longitude is west of the prime meridian. ''' - + import gdal width = ds.RasterXSize height = ds.RasterYSize gt = ds.GetGeoTransform() @@ -64,14 +68,14 @@ def load_DEM(filepath): return elevation, lat, long -def get_pixel_coords(lat, lon, DEM_path): +def get_pixel_coords(latitude, longitude, DEM_path): r'''Gets pixel coordinates from the raster, given latitude and longitude Parameters ---------- - lat : float + latitude : float Latitude of the point - long: float + longitude: float Longitude of the point DEM_path: string Path of the DEM .hgt file @@ -251,14 +255,14 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, dem_res : float, default 30 Resolution of the DEM. The default is SRTM's 30m. [m] view_distance : float, default 500 - Radius of the area of consideration. [m] + Radius of the area of consideration. [pixels] az_len : int, default 360 Number of samples on the perimeter of the circle. Returns ------- azimuth : np.array - Numpy array of shape (az_len). Linearly spaced from [0,360] + Numpy array of shape (az_len). Linearly spaced from [0,360) elevation_angles: np.array The calculated elevation values at each point of azimuth profile: @@ -267,7 +271,7 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, """ from skimage.draw import line - azimuth = np.linspace(0, 360, az_len) + azimuth = np.linspace(0, 359, num=az_len) # Use Midpoint Circle Algo/ Bresenham's circle pts = _bresenham_circ(view_distance, dem_pixel) @@ -312,150 +316,3 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, elevation_angles[az] = np.rad2deg(elv_ang) profile[az] = highest_elv return azimuth, elevation_angles, profile - - -def calculate_dtf(horizon_azimuths, horizon_angles, - surface_tilt, surface_azimuth): - r"""Author: JPalakapillyKWH - 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 - ---------- - 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 - 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 - 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. - 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. - A similar method is used in section II of [1] although it is looking at - both ground and sky diffuse irradiation. - [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 - """ - 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) - b = np.sin(tilt_rad) * np.sin(plane_az_rad) - c = np.cos(tilt_rad) - - # this gets either a float or an array of zeros - dtf = np.multiply(0.0, surface_tilt) - num_points = horizon_azimuths.shape[0] - for i in range(horizon_azimuths.shape[0]): - az = np.radians(horizon_azimuths[i]) - 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 - 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 - - -def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): - r"""Author: JPalakapillyKWH - 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. This comes from Equation 10 in [1] - 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 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 - -------- - 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. - [1] doi.org/10.1016/j.solener.2014.09.037 - """ - 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) - elevation_angle = 90.0 - declination - elevation_angle[mask] = 0.0 - - return elevation_angle - - -def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): - r'''Author: JPalakapillyKWH - 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 - ---------- - 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] != 360): - 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 - horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] - mask = solar_zenith > horizon_zenith - adjustment[mask] = 0 - return adjustment diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index d31cb87ccf..93448ffb3f 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -680,7 +680,7 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): df : pd.DataFrame Pandas dataframe of the retrived horizon ''' - res = requests.get(url +f 'printhorizon?lat={latitude}&lon={longitude}', + res = requests.get(url +f'printhorizon?lat={latitude}&lon={longitude}', proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) diff --git a/pvlib/shading.py b/pvlib/shading.py index 9479eb1739..89b5ae1cc7 100644 --- a/pvlib/shading.py +++ b/pvlib/shading.py @@ -191,3 +191,161 @@ def sky_diffuse_passias(masking_angle): Available at https://www.nrel.gov/docs/fy18osti/67399.pdf """ return 1 - cosd(masking_angle/2)**2 + + +def calculate_dtf(horizon_azimuths, horizon_angles, + surface_tilt, surface_azimuth): + r"""Author: JPalakapillyKWH + 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 + ---------- + 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 + 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 + 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. + + 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. + A similar method is used in section II of [1] although it is looking at + both ground and sky diffuse irradiation. + [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 + + This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + """ + 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) + b = np.sin(tilt_rad) * np.sin(plane_az_rad) + c = np.cos(tilt_rad) + + # this gets either a float or an array of zeros + dtf = np.multiply(0.0, surface_tilt) + num_points = horizon_azimuths.shape[0] + for i in range(horizon_azimuths.shape[0]): + az = np.radians(horizon_azimuths[i]) + 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 - 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 + + +def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): + r""" + 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. This comes from Equation 10 in [1] + 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 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 + -------- + 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. + [1] doi.org/10.1016/j.solener.2014.09.037 + + Notes + ------ + This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + """ + 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) + elevation_angle = 90.0 - declination + elevation_angle[mask] = 0.0 + + return elevation_angle + + +def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): + r''' + 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 + ---------- + 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. + + Notes + ------ + This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + ''' + adjustment = np.ones(solar_zenith.shape) + + if (horizon_angles.shape[0] != 360): + 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 + horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] + mask = solar_zenith > horizon_zenith + adjustment[mask] = 0 + return adjustment From 02082ee041f645d5176380155418c50e5a3f4c35 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 10:01:27 -0700 Subject: [PATCH 25/94] cleaning up stickler --- pvlib/horizon.py | 9 +++++---- pvlib/shading.py | 15 +++++++++------ 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 4a0e312cd8..b252fa274c 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -2,7 +2,8 @@ The 'horizon' module contains function definitions that retrive & calculate the surrounding horizon using DEM elevation data. -Optional dependencies for this module include gdal, osgeo, geoio, and scikit-image. +Optional dependencies for this module include +gdal, osgeo, geoio, and scikit-image. ''' import numpy as np from scipy.signal import resample @@ -26,9 +27,9 @@ def latlong(ds): Notes ------ - Latitude and longitude are in decimal degrees and negative longitude is west of the prime meridian. + Latitude and longitude are in decimal degrees + and negative longitude is west of the prime meridian. ''' - import gdal width = ds.RasterXSize height = ds.RasterYSize gt = ds.GetGeoTransform() @@ -88,7 +89,7 @@ def get_pixel_coords(latitude, longitude, DEM_path): ''' from geoio import GeoImage img = GeoImage(DEM_path) - return map(int, img.proj_to_raster(lon, lat)) + return map(int, img.proj_to_raster(longitude, latitude)) def _point_symmetry(x, y): diff --git a/pvlib/shading.py b/pvlib/shading.py index 89b5ae1cc7..557770a67f 100644 --- a/pvlib/shading.py +++ b/pvlib/shading.py @@ -225,7 +225,7 @@ def calculate_dtf(horizon_azimuths, horizon_angles, horizontal irradiance (DHI) to get the incident irradiance from 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 @@ -237,7 +237,8 @@ def calculate_dtf(horizon_azimuths, horizon_angles, both ground and sky diffuse irradiation. [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 - This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + This function was written by @JPalakapillyKWH + in an uncompleted pvlib-python pull request #758. """ if horizon_azimuths.shape[0] != horizon_angles.shape[0]: raise ValueError('azimuths and elevation_angles must be of the same' @@ -292,10 +293,11 @@ 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 - + Notes ------ - This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + This function was written by @JPalakapillyKWH + in an uncompleted pvlib-python pull request #758. """ tilt = np.radians(surface_tilt) bearing = np.radians(direction - surface_azimuth - 180.0) @@ -332,10 +334,11 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): 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. - + Notes ------ - This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. + This function was written by @JPalakapillyKWH + in an uncompleted pvlib-python pull request #758. ''' adjustment = np.ones(solar_zenith.shape) From 206f979ed5cdaa3694d2de939509888903854b6b Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 10:04:37 -0700 Subject: [PATCH 26/94] cleaning up stickler 2 --- pvlib/horizon.py | 4 ++-- pvlib/iotools/pvgis.py | 2 +- pvlib/shading.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index b252fa274c..3dc2002a7c 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -2,8 +2,8 @@ The 'horizon' module contains function definitions that retrive & calculate the surrounding horizon using DEM elevation data. -Optional dependencies for this module include -gdal, osgeo, geoio, and scikit-image. +Optional dependencies for this module include +gdal, osgeo, geoio, and scikit-image. ''' import numpy as np from scipy.signal import resample diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 93448ffb3f..f60e6500d0 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -680,7 +680,7 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): df : pd.DataFrame Pandas dataframe of the retrived horizon ''' - res = requests.get(url +f'printhorizon?lat={latitude}&lon={longitude}', + res = requests.get(url + f'printhorizon?lat={ latitude }&lon={ longitude }', proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) diff --git a/pvlib/shading.py b/pvlib/shading.py index 557770a67f..f9cf62287e 100644 --- a/pvlib/shading.py +++ b/pvlib/shading.py @@ -236,8 +236,8 @@ def calculate_dtf(horizon_azimuths, horizon_angles, A similar method is used in section II of [1] although it is looking at both ground and sky diffuse irradiation. [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 - - This function was written by @JPalakapillyKWH + + This function was written by @JPalakapillyKWH in an uncompleted pvlib-python pull request #758. """ if horizon_azimuths.shape[0] != horizon_angles.shape[0]: From 156e9645044d1342ab47ac42c32e5d2ea6868a13 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 10:08:26 -0700 Subject: [PATCH 27/94] cleaning up stickler 3 --- pvlib/iotools/pvgis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index f60e6500d0..dcc7880831 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -680,7 +680,8 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): df : pd.DataFrame Pandas dataframe of the retrived horizon ''' - res = requests.get(url + f'printhorizon?lat={ latitude }&lon={ longitude }', + res = requests.get(url + + f'printhorizon?lat={ latitude }&lon={ longitude }', proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) From c674bb5f14c1dc85d993331662bcd1732287bfc2 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 15:23:37 -0700 Subject: [PATCH 28/94] updated math for angle calculation --- pvlib/horizon.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 3dc2002a7c..2945c011d5 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -311,7 +311,10 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, high_x, high_y = tuple(highest_point) # convert from altitude in m to elevation degrees. - xdist = np.abs(highest_point[0]-x0) + #xdist = np.abs(highest_point[0]-x0) + x1, y1 = highest_point[0], highest_point[1] + abs_dist = np.array([x1, y1]) - np.array([x0, y0]) + xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) ydist = highest_elv - elevation[y0][x0] elv_ang = np.arctan2(ydist, xdist) elevation_angles[az] = np.rad2deg(elv_ang) From 3a2cab9a2427eb806de2a57b53eae920d0a738a3 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 7 Feb 2022 15:24:33 -0700 Subject: [PATCH 29/94] removed old comment --- pvlib/horizon.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 2945c011d5..faa9e612bf 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -311,7 +311,6 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, high_x, high_y = tuple(highest_point) # convert from altitude in m to elevation degrees. - #xdist = np.abs(highest_point[0]-x0) x1, y1 = highest_point[0], highest_point[1] abs_dist = np.array([x1, y1]) - np.array([x0, y0]) xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) From 78f6c1df9dd6724fc636ba3b9225521106a6f3e0 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 8 Feb 2022 10:27:32 -0700 Subject: [PATCH 30/94] updated angle maximization method --- pvlib/horizon.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index faa9e612bf..f4f486f77e 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -307,15 +307,16 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, elvs_on_line = elevation[rr, cc] idx = np.stack((rr, cc), axis=-1) highest_elv = np.max(elvs_on_line) - highest_point = idx[np.argmax(elvs_on_line)] - high_x, high_y = tuple(highest_point) - - # convert from altitude in m to elevation degrees. - x1, y1 = highest_point[0], highest_point[1] - abs_dist = np.array([x1, y1]) - np.array([x0, y0]) - xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) - ydist = highest_elv - elevation[y0][x0] - elv_ang = np.arctan2(ydist, xdist) - elevation_angles[az] = np.rad2deg(elv_ang) + largest_ang = 0 + for point in idx: + y1,x1 = point + # convert from altitude in m to elevation degrees. + abs_dist = np.array([x1, y1]) - np.array([x0, y0]) + xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) + ydist = elevation[y1][x1] - elevation[y0][x0] + elv_ang = np.arctan2(ydist, xdist) + if elv_ang > largest_ang: + largest_ang = elv_ang + elevation_angles[az] = np.rad2deg(largest_ang) profile[az] = highest_elv return azimuth, elevation_angles, profile From 641bf537867e6ad72ed975574f6b1258f1e96a3c Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 8 Feb 2022 10:28:30 -0700 Subject: [PATCH 31/94] one of these days I'll remember stickler before pushing --- pvlib/horizon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index f4f486f77e..2d990e4697 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -309,7 +309,7 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, highest_elv = np.max(elvs_on_line) largest_ang = 0 for point in idx: - y1,x1 = point + y1, x1 = point # convert from altitude in m to elevation degrees. abs_dist = np.array([x1, y1]) - np.array([x0, y0]) xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) From 3d310c62b533081e9b19764d5d5f5ee27e3d0d04 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 8 Feb 2022 10:41:54 -0700 Subject: [PATCH 32/94] changed variable names to make clear that we're working in pixels --- pvlib/horizon.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 2d990e4697..ec53435ac0 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -284,35 +284,35 @@ def horizon_map(dem_pixel, elevation, dem_res=30.0, elevation_angles = np.zeros(azimuth.shape) x_bound, y_bound = elevation.shape for az in range(az_len): - source_lon = pts[1][az] - source_lat = pts[0][az] + target_y = pts[1][az] + target_x = pts[0][az] # check for outmapping # for now, if out of the DEM, assign to the border - if source_lon >= y_bound: - source_lon = y_bound-1 - elif source_lon <= 0: - source_lon = 0 - if source_lat >= y_bound: - source_lat = y_bound-1 - elif source_lat <= 0: - source_lat = 0 - - target_lon = int(y0) - target_lat = int(x0) + if target_y >= y_bound: + target_y = y_bound-1 + elif target_y <= 0: + target_y = 0 + if target_x >= y_bound: + target_x = y_bound-1 + elif target_x <= 0: + target_x = 0 + + observer_pixel_y = int(y0) + observer_pixel_x = int(x0) # draw a line from observer to horizon and record points - rr, cc = line(source_lon, source_lat, target_lon, target_lat) + rr, cc = line(target_y, target_x, observer_pixel_y, observer_pixel_x) # get the highest elevation on the line elvs_on_line = elevation[rr, cc] - idx = np.stack((rr, cc), axis=-1) + points_on_line = np.stack((rr, cc), axis=-1) highest_elv = np.max(elvs_on_line) largest_ang = 0 - for point in idx: + for point in points_on_line: y1, x1 = point # convert from altitude in m to elevation degrees. abs_dist = np.array([x1, y1]) - np.array([x0, y0]) - xdist = np.linalg.norm(abs_dist) * np.sqrt(dem_res) + xdist = np.linalg.norm(abs_dist) * dem_res ydist = elevation[y1][x1] - elevation[y0][x0] elv_ang = np.arctan2(ydist, xdist) if elv_ang > largest_ang: From c7b65b3885fe7e398cf8e17c9ba6aab7a2e76ccb Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 9 Feb 2022 12:35:45 -0700 Subject: [PATCH 33/94] added module to download SRTM DEM data --- pvlib/iotools/cgiar.py | 78 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 pvlib/iotools/cgiar.py diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py new file mode 100644 index 0000000000..1cc0bb8b91 --- /dev/null +++ b/pvlib/iotools/cgiar.py @@ -0,0 +1,78 @@ +""" +Get, read, and parse data from `CGIAR GeoPortal `_. + +For more information, see the following links: +* `FAQ `_ +* `Interpolation Method `_ + +""" +import math +from zipfile import ZipFile +import requests +import io + +def _lat_lon_to_query(longitude, latitude, srtm_arc_sec = 3): + r'''Converts latitude, longitude from the format used as + input to the other functions to format used by CGIAR + ---------- + longitude : numeric + longitude value, negative W of prime meridian + latitude: numeric + latitude value + srtm_arc_sec: numeric {1,3} + The resolution (arc-seconds) of the desired DEM. Either 1 or 3 + + Returns + ------- + rounded_lon : numeric + Rounded/adjusted longitude value + rounded_lat : numeric + Rounded/adjusted latitude value + + ''' + rounded = (int(math.floor(longitude)), int(math.floor(latitude))) + if srtm_arc_sec == 1: + return rounded + elif srtm_arc_sec == 3: + rounded_lon, rounded_lat = rounded + return (rounded_lon + 180) // 5 + 1, (64 - rounded_lat) // 5 + else: + raise Execption("Please use SRTM 1 Arc Second or SRTM 3 Arc Second") + + +def download_SRTM(latitude, longitude, srtm_arc_sec=3, path_to_save="./", proxies=None): + r'''Downloads a SRTM DEM tile from CGIAR, + saves it to a path, and loads it as an array + ---------- + latitude: numeric + latitude value to be included in the DEM + longitude : numeric + longitude value to be included in the DEM, negative W of prime meridian + srtm_arc_sec: numeric, {1,3} + The resolution (arc-seconds) of the desired DEM. Either 1 or 3 + path_to_save: string + The base path to save the DEM as a .tif file + proxies: dict + Proxy table for a corporate network + + Returns + ------- + img : np.array + Numpy array filled with elevation values in [m] + path : string + Path and filename where the DEM .tif filed was saved + + ''' + import skimage + long, lat = _lat_lon_to_query(longitude, latitude, srtm_arc_sec) + base_url = 'https://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/' + query_URL = base_url + f'srtm_{ long:02d }_{ lat:02d }.zip' + res = requests.get(query_URL, proxies=proxies, verify=False) + res.raise_for_status() + zipfile = ZipFile(io.BytesIO(res.content)) + ext = '.tif' + files = zipfile.namelist() + file = [f for f in files if ext in f][0] + path= zipfile.extract(file, path=path_to_save) + img = skimage.io.imread(path) + return img, path From 288390ba9429f43ff940cd013e53ca5fcb4657ca Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 9 Feb 2022 12:39:58 -0700 Subject: [PATCH 34/94] spelling/spacing --- pvlib/iotools/cgiar.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 1cc0bb8b91..3315b5b29b 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -1,9 +1,10 @@ """ -Get, read, and parse data from `CGIAR GeoPortal `_. +Get, read, and parse data from +`CGIAR GeoPortal `_. For more information, see the following links: * `FAQ `_ -* `Interpolation Method `_ + """ import math @@ -11,7 +12,8 @@ import requests import io -def _lat_lon_to_query(longitude, latitude, srtm_arc_sec = 3): + +def _lat_lon_to_query(longitude, latitude, srtm_arc_sec=3): r'''Converts latitude, longitude from the format used as input to the other functions to format used by CGIAR ---------- @@ -20,7 +22,8 @@ def _lat_lon_to_query(longitude, latitude, srtm_arc_sec = 3): latitude: numeric latitude value srtm_arc_sec: numeric {1,3} - The resolution (arc-seconds) of the desired DEM. Either 1 or 3 + The resolution (arc-seconds) of the desired DEM. + Either 1 or 3 Returns ------- @@ -37,19 +40,22 @@ def _lat_lon_to_query(longitude, latitude, srtm_arc_sec = 3): rounded_lon, rounded_lat = rounded return (rounded_lon + 180) // 5 + 1, (64 - rounded_lat) // 5 else: - raise Execption("Please use SRTM 1 Arc Second or SRTM 3 Arc Second") + raise Exception("Please use SRTM 1 Arc Second or SRTM 3 Arc Second") + - -def download_SRTM(latitude, longitude, srtm_arc_sec=3, path_to_save="./", proxies=None): +def download_SRTM(latitude, longitude, srtm_arc_sec=3, + path_to_save="./", proxies=None): r'''Downloads a SRTM DEM tile from CGIAR, saves it to a path, and loads it as an array ---------- latitude: numeric latitude value to be included in the DEM longitude : numeric - longitude value to be included in the DEM, negative W of prime meridian + longitude value to be included in the DEM, + negative W of prime meridian srtm_arc_sec: numeric, {1,3} - The resolution (arc-seconds) of the desired DEM. Either 1 or 3 + The resolution (arc-seconds) of the desired DEM. + Either 1 or 3 path_to_save: string The base path to save the DEM as a .tif file proxies: dict From e8d8795b713e0bfd2fba414150edceb0aff90559 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 9 Feb 2022 12:42:01 -0700 Subject: [PATCH 35/94] stickler --- pvlib/iotools/cgiar.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 3315b5b29b..fbc2ae8359 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -1,5 +1,5 @@ """ -Get, read, and parse data from +Get, read, and parse data from `CGIAR GeoPortal `_. For more information, see the following links: @@ -43,7 +43,7 @@ def _lat_lon_to_query(longitude, latitude, srtm_arc_sec=3): raise Exception("Please use SRTM 1 Arc Second or SRTM 3 Arc Second") -def download_SRTM(latitude, longitude, srtm_arc_sec=3, +def download_SRTM(latitude, longitude, srtm_arc_sec=3, path_to_save="./", proxies=None): r'''Downloads a SRTM DEM tile from CGIAR, saves it to a path, and loads it as an array @@ -71,14 +71,15 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, ''' import skimage long, lat = _lat_lon_to_query(longitude, latitude, srtm_arc_sec) - base_url = 'https://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/' + org_url = 'https://srtm.csi.cgiar.org/' + base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' query_URL = base_url + f'srtm_{ long:02d }_{ lat:02d }.zip' res = requests.get(query_URL, proxies=proxies, verify=False) res.raise_for_status() zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [f for f in files if ext in f][0] + file = [ f for f in files if ext in f ][0] path= zipfile.extract(file, path=path_to_save) img = skimage.io.imread(path) return img, path From a2496e5cbd46c8b47cb231b4df9701a8b0ac0563 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 9 Feb 2022 12:43:00 -0700 Subject: [PATCH 36/94] stickler3 --- pvlib/iotools/cgiar.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index fbc2ae8359..c77dc51bce 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -79,7 +79,7 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [ f for f in files if ext in f ][0] - path= zipfile.extract(file, path=path_to_save) + file = [ f for f in files if ext in f][0] + path = zipfile.extract(file, path=path_to_save) img = skimage.io.imread(path) return img, path From d36c50800c873d0a63dea5d107c9a6a3c7569c27 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Fri, 4 Mar 2022 14:19:42 -0700 Subject: [PATCH 37/94] fixed f string formatting stickler doesn't like it but it doesn't work with spaces --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index dcc7880831..edb98b46a4 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -681,7 +681,7 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): Pandas dataframe of the retrived horizon ''' res = requests.get(url + - f'printhorizon?lat={ latitude }&lon={ longitude }', + f'printhorizon?lat={latitude}&lon={longitude}', proxies=proxies, verify=False) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) From 60ba9cb557ce4bc3b9600a8962eb4de5fb20e4f2 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 13 Oct 2022 08:44:18 -0700 Subject: [PATCH 38/94] added horizon retrival, functions, dem retrival --- pvlib/horizon.py | 417 +++++++++--------------------- pvlib/iotools/cgiar.py | 11 +- pvlib/iotools/pvgis.py | 9 +- pvlib/tests/iotools/test_pvgis.py | 17 +- 4 files changed, 153 insertions(+), 301 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index ec53435ac0..089b00bc86 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -1,322 +1,153 @@ ''' The 'horizon' module contains function definitions that retrive & calculate the surrounding horizon using DEM -elevation data. +elevation data and its impact on performance. Optional dependencies for this module include gdal, osgeo, geoio, and scikit-image. ''' import numpy as np -from scipy.signal import resample +from pathlib import Path - -def latlong(ds): - r'''From a gdal dataset, retrieve the geotransform - and return latitude and longitude coordinates - for a DEM GeoTiff. - - Parameters - ---------- - ds : tuple - Geotransform parameters given by a osgeo.gdal.Dataset - - Returns - ------- - lat,long : tuple - Tuple of np.array of latitude,longitude - corresponding to the pixels in the GeoTiff - - Notes - ------ - Latitude and longitude are in decimal degrees - and negative longitude is west of the prime meridian. - ''' - width = ds.RasterXSize - height = ds.RasterYSize - gt = ds.GetGeoTransform() - minx = gt[0] - miny = gt[3] + width*gt[4] + height*gt[5] - maxx = gt[0] + width*gt[1] + height*gt[2] - maxy = gt[3] - lat = np.linspace(miny, maxy, height) - long = np.linspace(minx, maxx, width) - return (lat, long) - - -def load_DEM(filepath): - r'''Loads a DEM from a .hgt file, segments into the GeoTiff - (elevation) and corresponding latitude and longitude - - Parameters - ---------- - filepath : string - Absolute or relative path to the .hgt file - - Returns - ------- - elevation : np.array - 2D numpy array, pixel value is elevation in meters - lat : np.array - latitudes corresponding to the pixels along the 0th dim of elevation - long : np.array - longitudes corresponding to the pixels along the 1th dim of elevation +def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' - import gdal - from osgeo import gdal_array - dataset = gdal.Open(filepath) - rasterArray = gdal_array.LoadFile(filepath) - elevation = rasterArray - lat, long = latlong(dataset) - return elevation, lat, long - - -def get_pixel_coords(latitude, longitude, DEM_path): - r'''Gets pixel coordinates from the raster, given latitude and longitude - + 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 ---------- - latitude : float - Latitude of the point - longitude: float - Longitude of the point - DEM_path: string - Path of the DEM .hgt file - + 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 ------- - coords : tuple - Tuple with two elements, containing the x and y coordinate - of the raster corresponding to latitiude and longitude + 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. ''' - from geoio import GeoImage - img = GeoImage(DEM_path) - return map(int, img.proj_to_raster(longitude, latitude)) + adjustment = np.ones(solar_zenith.shape) + if (horizon_angles.shape[0] != 360): + raise ValueError('horizon_angles must contain exactly 360 values' + '(for each degree of azimuth 0-359).') -def _point_symmetry(x, y): - r"""Reflect a point 8 ways to form a circle. + 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 + return adjustment - Parameters - ---------- - x : numeric - x point of the circle - y: numeric - y point of the circle - - Returns - ------- - points : list - List of reflected points +def calculate_dtf(horizon_azimuths, horizon_angles, + surface_tilt, surface_azimuth): """ - return [(x, y), - (y, x), - (-x, y), - (-y, x), - (x, -y), - (y, -x), - (-x, -y), - (-y, -x)] - - -def _bresenham_circ(r, center=(0, 0)): - r""" Use midpoint algorithum to build a rasterized circle. - Modified from - https://funloop.org/post/2021-03-15-bresenham-circle-drawing-algorithm.html#bresenhams-algorithm - to add circles not centered at the origin - + 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 ---------- - r : float - Radius of the circle - center: tuple - Center of the point (Cartesian) - + 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 + 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 + 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 ------- - points : np.array - Array of shape (n,2) of points that form a rasterized circle. - n depends on the radius of the circle. + 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. + 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. + A similar method is used in section II of [1] although it is looking at + both ground and sky diffuse irradiation. + [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 """ - points = [] - x = 0 - y = -r - F_M = 1 - r - # Initial value for (0,-r) for 2x + 3 = 0x + 3 = 3. - d_e = 3 - # Initial value for (0,-r) for 2(x + y) + 5 = 0 - 2y + 5 = -2y + 5. - d_ne = -(r << 1) + 5 - points.extend(_point_symmetry(x, y)) - while x < -y: - if F_M < 0: - F_M += d_e - else: - F_M += d_ne - d_ne += 2 - y += 1 - d_e += 2 - d_ne += 2 - x += 1 - points.extend(_point_symmetry(x, y)) - points = np.array(points) - newx = points[:, 0] + center[0] - newy = points[:, 1] + center[1] - return np.vstack((newx, newy)) - - -def _pol2cart(r, theta): - r'''Converts polar to cartesian - Parameters - ---------- - r : np.array - r values for the points - theta: np.array - theta values for the points. [radian] - - Returns - ------- - x : np.array - x values for the points - y : np.array - y values for the points - - ''' - - z = r * np.exp(1j * theta) - x, y = z.real, z.imag - return x, y - - -def _cart2pol(x, y): - r'''Converts cartesian to polar - Parameters - ---------- - x : np.array - x values for the points - y : np.array - y values for the points - - Returns - ------- - r : np.array - r values for the points - theta: np.array - theta values for the points. [radian] - - ''' - - z = x + y * 1j - r, theta = np.abs(z), np.angle(z) - return r, theta - - -def _sort_circ(pts, center=(0, 0), az_len=360): - r'''Sort and resample points on a circle such that the zeroth element is - due east and the points move around the circle counter clockwise. - While in polar domain, resample points using FFT to - obtain desired number of bins, typically 360, for - degrees around the circle. - - Parameters - ---------- - pts : np.array - Array of shape (n, 2) of points in Cartesian system - center : tuple - Center (x,y) of the circle rasterized by pts - az_len : int, default 360 - Desired number of points in the rasterized circle. - - Returns - ------- - points : np.array - Sorted and resampled points that comprise the circle - ''' - pts[0] -= center[0] - pts[1] -= center[1] - r, theta = _cart2pol(pts[0], pts[1]) - stacked = np.vstack((theta, r)) - sort = np.sort(stacked, axis=1) - sort = resample(sort, az_len, axis=1) - theta = sort[0] - r = sort[1] - x, y = _pol2cart(r, theta) - x += center[0] - y += center[1] - return np.vstack((x, y)).astype(int) - - -def horizon_map(dem_pixel, elevation, dem_res=30.0, - view_distance=500, az_len=36): - r"""Finds the horizon for the point ``dem_pixel``. - + 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) + b = np.sin(tilt_rad) * np.sin(plane_az_rad) + c = np.cos(tilt_rad) + + # this gets either a float or an array of zeros + dtf = np.multiply(0.0, surface_tilt) + num_points = horizon_azimuths.shape[0] + for i in range(horizon_azimuths.shape[0]): + az = np.radians(horizon_azimuths[i]) + 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 - 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 + +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. This comes from Equation 10 in [1] Parameters ---------- - dem_pixel : tuple (int,int) - Point on the DEM expressed in pixel coordinates - elevation : np.array - nxn DEM of elevation values - dem_res : float, default 30 - Resolution of the DEM. The default is SRTM's 30m. [m] - view_distance : float, default 500 - Radius of the area of consideration. [pixels] - az_len : int, default 360 - Number of samples on the perimeter of the circle. - + 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 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 - ------- - azimuth : np.array - Numpy array of shape (az_len). Linearly spaced from [0,360) - elevation_angles: np.array - The calculated elevation values at each point of azimuth - profile: - The largest elevation in meters on the line between the observer - and the highest point on the horizon within view_distance - + -------- + 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. + [1] doi.org/10.1016/j.solener.2014.09.037 """ - from skimage.draw import line - azimuth = np.linspace(0, 359, num=az_len) - - # Use Midpoint Circle Algo/ Bresenham's circle - pts = _bresenham_circ(view_distance, dem_pixel) - - # sort circle points and resample to desired number of points (az_len) - pts = _sort_circ(pts, center=dem_pixel, az_len=az_len) - x0, y0 = dem_pixel - profile = np.zeros(azimuth.shape) - elevation_angles = np.zeros(azimuth.shape) - x_bound, y_bound = elevation.shape - for az in range(az_len): - target_y = pts[1][az] - target_x = pts[0][az] - # check for outmapping - # for now, if out of the DEM, assign to the border - if target_y >= y_bound: - target_y = y_bound-1 - elif target_y <= 0: - target_y = 0 - if target_x >= y_bound: - target_x = y_bound-1 - elif target_x <= 0: - target_x = 0 - - observer_pixel_y = int(y0) - observer_pixel_x = int(x0) + tilt = np.radians(surface_tilt) + bearing = np.radians(direction - surface_azimuth - 180.0) - # draw a line from observer to horizon and record points - rr, cc = line(target_y, target_x, observer_pixel_y, observer_pixel_x) + declination = np.degrees(np.arctan(1.0/np.tan(tilt)/np.cos(bearing))) + mask = (declination <= 0) + elevation_angle = 90.0 - declination + elevation_angle[mask] = 0.0 - # get the highest elevation on the line - elvs_on_line = elevation[rr, cc] - points_on_line = np.stack((rr, cc), axis=-1) - highest_elv = np.max(elvs_on_line) - largest_ang = 0 - for point in points_on_line: - y1, x1 = point - # convert from altitude in m to elevation degrees. - abs_dist = np.array([x1, y1]) - np.array([x0, y0]) - xdist = np.linalg.norm(abs_dist) * dem_res - ydist = elevation[y1][x1] - elevation[y0][x0] - elv_ang = np.arctan2(ydist, xdist) - if elv_ang > largest_ang: - largest_ang = elv_ang - elevation_angles[az] = np.rad2deg(largest_ang) - profile[az] = highest_elv - return azimuth, elevation_angles, profile + return elevation_angle diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index c77dc51bce..5de4b36e21 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -70,16 +70,19 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, ''' import skimage + # convert lat/lon to format cgiar uses long, lat = _lat_lon_to_query(longitude, latitude, srtm_arc_sec) org_url = 'https://srtm.csi.cgiar.org/' base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' - query_URL = base_url + f'srtm_{ long:02d }_{ lat:02d }.zip' - res = requests.get(query_URL, proxies=proxies, verify=False) + query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' + res = requests.get(query_URL, proxies = proxies, verify = False) res.raise_for_status() + + # unzip download zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [ f for f in files if ext in f][0] - path = zipfile.extract(file, path=path_to_save) + file = [ f for f in files if ext in f ][0] + path = zipfile.extract(file, path = path_to_save) img = skimage.io.imread(path) return img, path diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index edb98b46a4..e6108828b5 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -688,8 +688,11 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): # the horizon data is given in a different format then the others # numpy has an easier time decoding it array = np.genfromtxt(io.StringIO(string), - skip_header=4, skip_footer=7) + skip_header = 4, skip_footer = 7) df = pd.DataFrame(array) - # Get the column names - df.columns = [s for s in string.splitlines()[3].split('\t') if s] + + # Set the column names + df.columns = ['horizon_azimuth', 'horizon_angles', + 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', + 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] return df diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index a7b0cbbd0a..52e9ef93e2 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -9,6 +9,7 @@ import requests from pvlib.iotools import get_pvgis_tmy, read_pvgis_tmy from pvlib.iotools import get_pvgis_hourly, read_pvgis_hourly +from pvlib.iotools import get_pvgis_horizon from ..conftest import (DATA_DIR, RERUNS, RERUNS_DELAY, assert_frame_equal, fail_on_pvlib_version) from pvlib._deprecation import pvlibDeprecationWarning @@ -65,6 +66,15 @@ [1586.9, 80.76, 83.95, 9.04, 13.28, 2.79, 1.36, 0.0], [713.3, 5.18, 70.57, 7.31, 18.56, 3.66, 1.27, 0.0]] +data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, + 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, + 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, + 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, + 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.1, 1.9, 3.8, 5.0, + 6.5, 9.2, 9.9] + inputs_radiation_csv = {'latitude': 45.0, 'longitude': 8.0, 'elevation': 250.0, 'radiation_database': 'PVGIS-SARAH', 'Slope': '30 deg.', 'Azimuth': '0 deg.'} @@ -508,7 +518,11 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): actual, _, _, _ = get_pvgis_tmy(45, 8, map_variables=True) assert all([c in pvgis_tmy_mapped_columns for c in actual.columns]) - +@pytest.mark.remote_data +def test_read_pvgis_horizon(): + azimuth, horizon = get_pvgis_horizon(35.171051, -106.465158) + assert all(np.isclose(horizon, data_horizon_abq)) + def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' actual, _, _, _ = read_pvgis_tmy(fn, map_variables=True) @@ -588,3 +602,4 @@ def test_read_pvgis_tmy_exception(): with pytest.raises(ValueError, match=err_msg): read_pvgis_tmy('filename', pvgis_format=bad_outputformat, map_variables=False) + \ No newline at end of file From 709493e96600326fb218c6882174c70cfd75af95 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 13 Oct 2022 08:59:47 -0700 Subject: [PATCH 39/94] stickler --- pvlib/horizon.py | 4 +- pvlib/iotools/cgiar.py | 8 +- pvlib/iotools/pvgis.py | 58 ++++++-------- pvlib/tests/iotools/test_pvgis.py | 127 +++++++++++++++--------------- 4 files changed, 97 insertions(+), 100 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 089b00bc86..8558a6d7a8 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -6,7 +6,7 @@ gdal, osgeo, geoio, and scikit-image. ''' import numpy as np -from pathlib import Path + def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' @@ -46,6 +46,7 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): adjustment[mask] = 0 return adjustment + def calculate_dtf(horizon_azimuths, horizon_angles, surface_tilt, surface_azimuth): """ @@ -115,6 +116,7 @@ def calculate_dtf(horizon_azimuths, horizon_angles, dtf += 2 * (first_term + second_term) / num_points return dtf + def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ Determine the elevation angle created by the surface of a tilted plane diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 5de4b36e21..3701f85682 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -75,14 +75,14 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, org_url = 'https://srtm.csi.cgiar.org/' base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' - res = requests.get(query_URL, proxies = proxies, verify = False) + res = requests.get(query_URL, proxies=proxies, verify=False) res.raise_for_status() - + # unzip download zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [ f for f in files if ext in f ][0] - path = zipfile.extract(file, path = path_to_save) + file = [f for f in files if ext in f][0] + path = zipfile.extract(file, path=path_to_save) img = skimage.io.imread(path) return img, path diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 9c904c852b..005410f6c4 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -27,7 +27,7 @@ URL = 'https://re.jrc.ec.europa.eu/api/' # Dictionary mapping PVGIS names to pvlib names -VARIABLE_MAP = { +PVGIS_VARIABLE_MAP = { 'G(h)': 'ghi', 'Gb(n)': 'dni', 'Gd(h)': 'dhi', @@ -113,11 +113,10 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, documentation [2]_ for more info. url: str, default: :const:`pvlib.iotools.pvgis.URL` Base url of PVGIS API. ``seriescalc`` is appended to get hourly data - endpoint. Note, a specific PVGIS version can be specified, e.g., - https://re.jrc.ec.europa.eu/api/v5_2/ + endpoint. map_variables: bool, default: True When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable :const:`VARIABLE_MAP`. + where applicable. See variable PVGIS_VARIABLE_MAP. timeout: int, default: 30 Time in seconds to wait for server response before timeout @@ -140,10 +139,10 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, Hint ---- PVGIS provides access to a number of different solar radiation datasets, - including satellite-based (SARAH, SARAH2, and NSRDB PSM3) and re-analysis - products (ERA5). Each data source has a different geographical coverage and - time stamp convention, e.g., SARAH and SARAH2 provide instantaneous values, - whereas values from ERA5 are averages for the hour. + including satellite-based (SARAH, CMSAF, and NSRDB PSM3) and re-analysis + products (ERA5 and COSMO). Each data source has a different geographical + coverage and time stamp convention, e.g., SARAH and CMSAF provide + instantaneous values, whereas values from ERA5 are averages for the hour. Notes ----- @@ -174,12 +173,6 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, -------- pvlib.iotools.read_pvgis_hourly, pvlib.iotools.get_pvgis_tmy - Examples - -------- - >>> # Retrieve two years of irradiance data from PVGIS: - >>> data, meta, inputs = pvlib.iotools.get_pvgis_hourly( # doctest: +SKIP - >>> latitude=45, longitude=8, start=2015, end=2016) # doctest: +SKIP - References ---------- .. [1] `PVGIS `_ @@ -236,7 +229,7 @@ def _parse_pvgis_hourly_json(src, map_variables): data = data.drop('time', axis=1) data = data.astype(dtype={'Int': 'int'}) # The 'Int' column to be integer if map_variables: - data = data.rename(columns=VARIABLE_MAP) + data = data.rename(columns=PVGIS_VARIABLE_MAP) return data, inputs, metadata @@ -278,7 +271,7 @@ def _parse_pvgis_hourly_csv(src, map_variables): data.index = pd.to_datetime(data['time'], format='%Y%m%d:%H%M', utc=True) data = data.drop('time', axis=1) if map_variables: - data = data.rename(columns=VARIABLE_MAP) + data = data.rename(columns=PVGIS_VARIABLE_MAP) # All columns should have the dtype=float, except 'Int' which should be # integer. It is necessary to convert to float, before converting to int data = data.astype(float).astype(dtype={'Int': 'int'}) @@ -305,7 +298,7 @@ def read_pvgis_hourly(filename, pvgis_format=None, map_variables=True): ``pvgis_format`` is required and must be in ``['csv', 'json']``. map_variables: bool, default True When true, renames columns of the DataFrame to pvlib variable names - where applicable. See variable :const:`VARIABLE_MAP`. + where applicable. See variable PVGIS_VARIABLE_MAP. Returns ------- @@ -377,9 +370,8 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, userhorizon=None, startyear=None, endyear=None, url=URL, map_variables=None, timeout=30): """ - Get TMY data from PVGIS. - - For more information see the PVGIS [1]_ TMY tool documentation [2]_. + Get TMY data from PVGIS. For more information see the PVGIS [1]_ TMY tool + documentation [2]_. Parameters ---------- @@ -405,7 +397,7 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, base url of PVGIS API, append ``tmy`` to get TMY endpoint map_variables: bool When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable const:`VARIABLE_MAP`. + where applicable. See variable PVGIS_VARIABLE_MAP. timeout : int, default 30 time in seconds to wait for server response before timeout @@ -437,12 +429,13 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, the error message in the response will be raised as an exception, otherwise raise whatever ``HTTP/1.1`` error occurred - See Also + See also -------- read_pvgis_tmy References ---------- + .. [1] `PVGIS `_ .. [2] `PVGIS TMY tool `_ .. [3] `PVGIS horizon profile tool @@ -500,7 +493,7 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, ) map_variables = False if map_variables: - data = data.rename(columns=VARIABLE_MAP) + data = data.rename(columns=PVGIS_VARIABLE_MAP) return data, months_selected, inputs, meta @@ -574,7 +567,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): be in ``['csv', 'epw', 'json', 'basic']``. map_variables: bool When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable :const:`VARIABLE_MAP`. + where applicable. See variable PVGIS_VARIABLE_MAP. Returns @@ -592,12 +585,12 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): ------ ValueError if ``pvgis_format`` is ``None`` and the file extension is neither - ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided - as input but isn't in ``['csv', 'epw', 'json', 'basic']`` + ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided as + input but isn't in ``['csv', 'epw', 'json', 'basic']`` TypeError if ``pvgis_format`` is ``None`` and ``filename`` is a buffer - See Also + See also -------- get_pvgis_tmy """ @@ -663,10 +656,11 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): ) map_variables = False if map_variables: - data = data.rename(columns=VARIABLE_MAP) + data = data.rename(columns=PVGIS_VARIABLE_MAP) return data, months_selected, inputs, meta + def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): r''' Get horizon data from PVGIS @@ -694,11 +688,11 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): # the horizon data is given in a different format then the others # numpy has an easier time decoding it array = np.genfromtxt(io.StringIO(string), - skip_header = 4, skip_footer = 7) + skip_header=4, skip_footer=7) df = pd.DataFrame(array) - + # Set the column names df.columns = ['horizon_azimuth', 'horizon_angles', - 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', - 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] + 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', + 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] return df diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 7b2565fe69..7e50e8dbde 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -21,12 +21,12 @@ testfile_radiation_csv = DATA_DIR / \ 'pvgis_hourly_Timeseries_45.000_8.000_SA_30deg_0deg_2016_2016.csv' testfile_pv_json = DATA_DIR / \ - 'pvgis_hourly_Timeseries_45.000_8.000_SA2_10kWp_CIS_5_2a_2013_2014.json' + 'pvgis_hourly_Timeseries_45.000_8.000_CM_10kWp_CIS_5_2a_2013_2014.json' index_radiation_csv = \ pd.date_range('20160101 00:10', freq='1h', periods=14, tz='UTC') index_pv_json = \ - pd.date_range('2013-01-01 00:10', freq='1h', periods=10, tz='UTC') + pd.date_range('2013-01-01 00:55', freq='1h', periods=10, tz='UTC') columns_radiation_csv = [ 'Gb(i)', 'Gd(i)', 'Gr(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] @@ -34,9 +34,10 @@ 'poa_direct', 'poa_sky_diffuse', 'poa_ground_diffuse', 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] columns_pv_json = [ - 'P', 'G(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] + 'P', 'Gb(i)', 'Gd(i)', 'Gr(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] columns_pv_json_mapped = [ - 'P', 'poa_global', 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] + 'P', 'poa_direct', 'poa_sky_diffuse', 'poa_ground_diffuse', + 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] data_radiation_csv = [ [0.0, 0.0, 0.0, 0.0, 3.44, 1.43, 0.0], @@ -54,21 +55,21 @@ [4.25, 1.88, 0.05, 21.41, 7.84, 0.4, 1.0], [0.0, 0.0, 0.0, 0.0, 7.43, 0.72, 0.0]] data_pv_json = [ - [0.0, 0.0, 0.0, -0.97, 1.52, 0.0], - [0.0, 0.0, 0.0, -1.06, 1.45, 0.0], - [0.0, 0.0, 0.0, -1.03, 1.45, 0.0], - [0.0, 0.0, 0.0, -0.48, 1.31, 0.0], - [0.0, 0.0, 0.0, -0.09, 1.24, 0.0], - [0.0, 0.0, 0.0, -0.38, 1.17, 0.0], - [0.0, 0.0, 0.0, 0.29, 1.03, 0.0], - [0.0, 0.0, 0.0, 1.0, 0.62, 0.0], - [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], - [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] - -data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, - 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, - 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, - 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, + [0.0, 0.0, 0.0, 0.0, 0.0, 3.01, 1.23, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 2.22, 1.46, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 1.43, 1.7, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.64, 1.93, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.77, 1.8, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.91, 1.66, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 1.05, 1.53, 0.0], + [3464.5, 270.35, 91.27, 6.09, 6.12, 1.92, 1.44, 0.0], + [1586.9, 80.76, 83.95, 9.04, 13.28, 2.79, 1.36, 0.0], + [713.3, 5.18, 70.57, 7.31, 18.56, 3.66, 1.27, 0.0]] + +data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, + 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, + 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, + 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1, 1.9, 3.8, 5.0, @@ -89,7 +90,7 @@ inputs_pv_json = { 'location': {'latitude': 45.0, 'longitude': 8.0, 'elevation': 250.0}, - 'meteo_data': {'radiation_db': 'PVGIS-SARAH2', 'meteo_db': 'ERA-Interim', + 'meteo_data': {'radiation_db': 'PVGIS-CMSAF', 'meteo_db': 'ERA-Interim', 'year_min': 2013, 'year_max': 2014, 'use_horizon': True, 'horizon_db': None, 'horizon_data': 'DEM-calculated'}, 'mounting_system': {'two_axis': { @@ -97,45 +98,45 @@ 'azimuth': {'value': '-', 'optimal': '-'}}}, 'pv_module': {'technology': 'CIS', 'peak_power': 10.0, 'system_loss': 5.0}} - metadata_pv_json = { 'inputs': { - 'location': - {'description': 'Selected location', 'variables': { - 'latitude': {'description': 'Latitude', 'units': 'decimal degree'}, # noqa: E501 - 'longitude': {'description': 'Longitude', 'units': 'decimal degree'}, # noqa: E501 - 'elevation': {'description': 'Elevation', 'units': 'm'}}}, - 'meteo_data': { - 'description': 'Sources of meteorological data', - 'variables': { - 'radiation_db': {'description': 'Solar radiation database'}, # noqa: E501 - 'meteo_db': {'description': 'Database used for meteorological variables other than solar radiation'}, # noqa: E501 - 'year_min': {'description': 'First year of the calculations'}, # noqa: E501 - 'year_max': {'description': 'Last year of the calculations'}, # noqa: E501 - 'use_horizon': {'description': 'Include horizon shadows'}, - 'horizon_db': {'description': 'Source of horizon data'}}}, - 'mounting_system': { - 'description': 'Mounting system', - 'choices': 'fixed, vertical_axis, inclined_axis, two_axis', - 'fields': { - 'slope': {'description': 'Inclination angle from the horizontal plane', 'units': 'degree'}, # noqa: E501 - 'azimuth': {'description': 'Orientation (azimuth) angle of the (fixed) PV system (0 = S, 90 = W, -90 = E)', 'units': 'degree'}}}, # noqa: E501 - 'pv_module': { - 'description': 'PV module parameters', - 'variables': { - 'technology': {'description': 'PV technology'}, - 'peak_power': {'description': 'Nominal (peak) power of the PV module', 'units': 'kW'}, # noqa: E501 - 'system_loss': {'description': 'Sum of system losses', 'units': '%'}}}}, # noqa: E501 - 'outputs': { - 'hourly': { - 'type': 'time series', 'timestamp': 'hourly averages', - 'variables': { - 'P': {'description': 'PV system power', 'units': 'W'}, - 'G(i)': {'description': 'Global irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 - 'H_sun': {'description': 'Sun height', 'units': 'degree'}, - 'T2m': {'description': '2-m air temperature', 'units': 'degree Celsius'}, # noqa: E501 - 'WS10m': {'description': '10-m total wind speed', 'units': 'm/s'}, # noqa: E501 - 'Int': {'description': '1 means solar radiation values are reconstructed'}}}}} # noqa: E501 + 'location': {'description': 'Selected location', 'variables': { + 'latitude': {'description': 'Latitude', 'units': 'decimal degree'}, + 'longitude': {'description': 'Longitude', 'units': 'decimal degree'}, # noqa: E501 + 'elevation': {'description': 'Elevation', 'units': 'm'}}}, + 'meteo_data': { + 'description': 'Sources of meteorological data', + 'variables': { + 'radiation_db': {'description': 'Solar radiation database'}, + 'meteo_db': {'description': 'Database used for meteorological variables other than solar radiation'}, # noqa: E501 + 'year_min': {'description': 'First year of the calculations'}, + 'year_max': {'description': 'Last year of the calculations'}, + 'use_horizon': {'description': 'Include horizon shadows'}, + 'horizon_db': {'description': 'Source of horizon data'}}}, + 'mounting_system': { + 'description': 'Mounting system', + 'choices': 'fixed, vertical_axis, inclined_axis, two_axis', + 'fields': { + 'slope': {'description': 'Inclination angle from the horizontal plane', 'units': 'degree'}, # noqa: E501 + 'azimuth': {'description': 'Orientation (azimuth) angle of the (fixed) PV system (0 = S, 90 = W, -90 = E)', 'units': 'degree'}}}, # noqa: E501 + 'pv_module': { + 'description': 'PV module parameters', + 'variables': { + 'technology': {'description': 'PV technology'}, + 'peak_power': {'description': 'Nominal (peak) power of the PV module', 'units': 'kW'}, # noqa: E501 + 'system_loss': {'description': 'Sum of system losses', 'units': '%'}}}}, # noqa: E501 + 'outputs': { + 'hourly': { + 'type': 'time series', 'timestamp': 'hourly averages', + 'variables': { + 'P': {'description': 'PV system power', 'units': 'W'}, + 'Gb(i)': {'description': 'Beam (direct) irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 + 'Gd(i)': {'description': 'Diffuse irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 + 'Gr(i)': {'description': 'Reflected irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 + 'H_sun': {'description': 'Sun height', 'units': 'degree'}, + 'T2m': {'description': '2-m air temperature', 'units': 'degree Celsius'}, # noqa: E501 + 'WS10m': {'description': '10-m total wind speed', 'units': 'm/s'}, # noqa: E501 + 'Int': {'description': '1 means solar radiation values are reconstructed'}}}}} # noqa: E501 def generate_expected_dataframe(values, columns, index): @@ -224,13 +225,12 @@ def test_read_pvgis_hourly_bad_extension(): args_pv_json = { 'surface_tilt': 30, 'surface_azimuth': 0, 'outputformat': 'json', - 'usehorizon': True, 'userhorizon': None, 'raddatabase': 'PVGIS-SARAH2', + 'usehorizon': True, 'userhorizon': None, 'raddatabase': 'PVGIS-CMSAF', 'start': pd.Timestamp(2013, 1, 1), 'end': pd.Timestamp(2014, 5, 1), 'pvcalculation': True, 'peakpower': 10, 'pvtechchoice': 'CIS', 'loss': 5, - 'trackingtype': 2, 'optimalangles': True, 'components': False, - 'url': 'https://re.jrc.ec.europa.eu/api/v5_2/'} + 'trackingtype': 2, 'optimalangles': True, 'components': True} -url_pv_json = 'https://re.jrc.ec.europa.eu/api/v5_2/seriescalc?lat=45&lon=8&outputformat=json&angle=30&aspect=0&pvtechchoice=CIS&mountingplace=free&trackingtype=2&components=0&usehorizon=1&raddatabase=PVGIS-SARAH2&startyear=2013&endyear=2014&pvcalculation=1&peakpower=10&loss=5&optimalangles=1' # noqa: E501 +url_pv_json = 'https://re.jrc.ec.europa.eu/api/seriescalc?lat=45&lon=8&outputformat=json&angle=30&aspect=0&pvtechchoice=CIS&mountingplace=free&trackingtype=2&components=1&usehorizon=1&raddatabase=PVGIS-CMSAF&startyear=2013&endyear=2014&pvcalculation=1&peakpower=10&loss=5&optimalangles=1' # noqa: E501 @pytest.mark.parametrize('testfile,expected_name,args,map_variables,url_test', [ # noqa: E501 @@ -518,11 +518,13 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): actual, _, _, _ = get_pvgis_tmy(45, 8, map_variables=True) assert all([c in pvgis_tmy_mapped_columns for c in actual.columns]) + @pytest.mark.remote_data def test_read_pvgis_horizon(): azimuth, horizon = get_pvgis_horizon(35.171051, -106.465158) assert all(np.isclose(horizon, data_horizon_abq)) - + + def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' actual, _, _, _ = read_pvgis_tmy(fn, map_variables=True) @@ -602,4 +604,3 @@ def test_read_pvgis_tmy_exception(): with pytest.raises(ValueError, match=err_msg): read_pvgis_tmy('filename', pvgis_format=bad_outputformat, map_variables=False) - \ No newline at end of file From fa14c67772446e35374cb8622310ebc2b4d9cfdf Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 13 Oct 2022 09:02:37 -0700 Subject: [PATCH 40/94] stickler2 --- pvlib/iotools/pvgis.py | 10 ++++++---- pvlib/tests/iotools/test_pvgis.py | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 005410f6c4..fd981689f0 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -585,8 +585,8 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): ------ ValueError if ``pvgis_format`` is ``None`` and the file extension is neither - ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided as - input but isn't in ``['csv', 'epw', 'json', 'basic']`` + ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided + as input but isn't in ``['csv', 'epw', 'json', 'basic']`` TypeError if ``pvgis_format`` is ``None`` and ``filename`` is a buffer @@ -693,6 +693,8 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): # Set the column names df.columns = ['horizon_azimuth', 'horizon_angles', - 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', - 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] + 'azimuth_sun_winter_solstice', + 'elevation_sun_winter_solstice', + 'azimuth_sun_summer_solstice', + 'elevation_sun_summer_solstice'] return df diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 7e50e8dbde..e6b57f6bc8 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -69,7 +69,7 @@ data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, - 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6 + 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1, 1.9, 3.8, 5.0, From 97ee479f5f899de48d750633632c13368955e448 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 13 Oct 2022 09:03:20 -0700 Subject: [PATCH 41/94] stickler3 --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index fd981689f0..f29ef5ac65 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -585,7 +585,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): ------ ValueError if ``pvgis_format`` is ``None`` and the file extension is neither - ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided + ``.csv``, ``.json``, nor ``.epw``, or if ``pvgis_format`` is provided as input but isn't in ``['csv', 'epw', 'json', 'basic']`` TypeError if ``pvgis_format`` is ``None`` and ``filename`` is a buffer From 8e01b570f5dfde67ace9d750fc96434b4907f834 Mon Sep 17 00:00:00 2001 From: Pierce Date: Thu, 13 Oct 2022 16:12:54 -0600 Subject: [PATCH 42/94] Revert "merge2" This reverts commit 6f9ad5e94f1f73045e6acc2f6a84983de197cddd, reversing changes made to 43951a228fb23249a9b8d00230276be1f0a5eea4. --- pvlib/horizon.py | 4 +- pvlib/iotools/cgiar.py | 8 +- pvlib/iotools/pvgis.py | 55 +++++++------ pvlib/tests/iotools/test_pvgis.py | 125 +++++++++++++++--------------- 4 files changed, 97 insertions(+), 95 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 8558a6d7a8..089b00bc86 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -6,7 +6,7 @@ gdal, osgeo, geoio, and scikit-image. ''' import numpy as np - +from pathlib import Path def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' @@ -46,7 +46,6 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): adjustment[mask] = 0 return adjustment - def calculate_dtf(horizon_azimuths, horizon_angles, surface_tilt, surface_azimuth): """ @@ -116,7 +115,6 @@ def calculate_dtf(horizon_azimuths, horizon_angles, dtf += 2 * (first_term + second_term) / num_points return dtf - def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ Determine the elevation angle created by the surface of a tilted plane diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 3701f85682..5de4b36e21 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -75,14 +75,14 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, org_url = 'https://srtm.csi.cgiar.org/' base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' - res = requests.get(query_URL, proxies=proxies, verify=False) + res = requests.get(query_URL, proxies = proxies, verify = False) res.raise_for_status() - + # unzip download zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [f for f in files if ext in f][0] - path = zipfile.extract(file, path=path_to_save) + file = [ f for f in files if ext in f ][0] + path = zipfile.extract(file, path = path_to_save) img = skimage.io.imread(path) return img, path diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index f29ef5ac65..c4c9575c3f 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -27,7 +27,7 @@ URL = 'https://re.jrc.ec.europa.eu/api/' # Dictionary mapping PVGIS names to pvlib names -PVGIS_VARIABLE_MAP = { +VARIABLE_MAP = { 'G(h)': 'ghi', 'Gb(n)': 'dni', 'Gd(h)': 'dhi', @@ -113,10 +113,11 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, documentation [2]_ for more info. url: str, default: :const:`pvlib.iotools.pvgis.URL` Base url of PVGIS API. ``seriescalc`` is appended to get hourly data - endpoint. + endpoint. Note, a specific PVGIS version can be specified, e.g., + https://re.jrc.ec.europa.eu/api/v5_2/ map_variables: bool, default: True When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable PVGIS_VARIABLE_MAP. + where applicable. See variable :const:`VARIABLE_MAP`. timeout: int, default: 30 Time in seconds to wait for server response before timeout @@ -139,10 +140,10 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, Hint ---- PVGIS provides access to a number of different solar radiation datasets, - including satellite-based (SARAH, CMSAF, and NSRDB PSM3) and re-analysis - products (ERA5 and COSMO). Each data source has a different geographical - coverage and time stamp convention, e.g., SARAH and CMSAF provide - instantaneous values, whereas values from ERA5 are averages for the hour. + including satellite-based (SARAH, SARAH2, and NSRDB PSM3) and re-analysis + products (ERA5). Each data source has a different geographical coverage and + time stamp convention, e.g., SARAH and SARAH2 provide instantaneous values, + whereas values from ERA5 are averages for the hour. Notes ----- @@ -173,6 +174,12 @@ def get_pvgis_hourly(latitude, longitude, start=None, end=None, -------- pvlib.iotools.read_pvgis_hourly, pvlib.iotools.get_pvgis_tmy + Examples + -------- + >>> # Retrieve two years of irradiance data from PVGIS: + >>> data, meta, inputs = pvlib.iotools.get_pvgis_hourly( # doctest: +SKIP + >>> latitude=45, longitude=8, start=2015, end=2016) # doctest: +SKIP + References ---------- .. [1] `PVGIS `_ @@ -229,7 +236,7 @@ def _parse_pvgis_hourly_json(src, map_variables): data = data.drop('time', axis=1) data = data.astype(dtype={'Int': 'int'}) # The 'Int' column to be integer if map_variables: - data = data.rename(columns=PVGIS_VARIABLE_MAP) + data = data.rename(columns=VARIABLE_MAP) return data, inputs, metadata @@ -271,7 +278,7 @@ def _parse_pvgis_hourly_csv(src, map_variables): data.index = pd.to_datetime(data['time'], format='%Y%m%d:%H%M', utc=True) data = data.drop('time', axis=1) if map_variables: - data = data.rename(columns=PVGIS_VARIABLE_MAP) + data = data.rename(columns=VARIABLE_MAP) # All columns should have the dtype=float, except 'Int' which should be # integer. It is necessary to convert to float, before converting to int data = data.astype(float).astype(dtype={'Int': 'int'}) @@ -298,7 +305,7 @@ def read_pvgis_hourly(filename, pvgis_format=None, map_variables=True): ``pvgis_format`` is required and must be in ``['csv', 'json']``. map_variables: bool, default True When true, renames columns of the DataFrame to pvlib variable names - where applicable. See variable PVGIS_VARIABLE_MAP. + where applicable. See variable :const:`VARIABLE_MAP`. Returns ------- @@ -370,8 +377,9 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, userhorizon=None, startyear=None, endyear=None, url=URL, map_variables=None, timeout=30): """ - Get TMY data from PVGIS. For more information see the PVGIS [1]_ TMY tool - documentation [2]_. + Get TMY data from PVGIS. + + For more information see the PVGIS [1]_ TMY tool documentation [2]_. Parameters ---------- @@ -397,7 +405,7 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, base url of PVGIS API, append ``tmy`` to get TMY endpoint map_variables: bool When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable PVGIS_VARIABLE_MAP. + where applicable. See variable const:`VARIABLE_MAP`. timeout : int, default 30 time in seconds to wait for server response before timeout @@ -429,13 +437,12 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, the error message in the response will be raised as an exception, otherwise raise whatever ``HTTP/1.1`` error occurred - See also + See Also -------- read_pvgis_tmy References ---------- - .. [1] `PVGIS `_ .. [2] `PVGIS TMY tool `_ .. [3] `PVGIS horizon profile tool @@ -493,7 +500,7 @@ def get_pvgis_tmy(latitude, longitude, outputformat='json', usehorizon=True, ) map_variables = False if map_variables: - data = data.rename(columns=PVGIS_VARIABLE_MAP) + data = data.rename(columns=VARIABLE_MAP) return data, months_selected, inputs, meta @@ -567,7 +574,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): be in ``['csv', 'epw', 'json', 'basic']``. map_variables: bool When true, renames columns of the Dataframe to pvlib variable names - where applicable. See variable PVGIS_VARIABLE_MAP. + where applicable. See variable :const:`VARIABLE_MAP`. Returns @@ -590,7 +597,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): TypeError if ``pvgis_format`` is ``None`` and ``filename`` is a buffer - See also + See Also -------- get_pvgis_tmy """ @@ -656,7 +663,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): ) map_variables = False if map_variables: - data = data.rename(columns=PVGIS_VARIABLE_MAP) + data = data.rename(columns=VARIABLE_MAP) return data, months_selected, inputs, meta @@ -688,13 +695,11 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): # the horizon data is given in a different format then the others # numpy has an easier time decoding it array = np.genfromtxt(io.StringIO(string), - skip_header=4, skip_footer=7) + skip_header = 4, skip_footer = 7) df = pd.DataFrame(array) - + # Set the column names df.columns = ['horizon_azimuth', 'horizon_angles', - 'azimuth_sun_winter_solstice', - 'elevation_sun_winter_solstice', - 'azimuth_sun_summer_solstice', - 'elevation_sun_summer_solstice'] + 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', + 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] return df diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index e6b57f6bc8..7b2565fe69 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -21,12 +21,12 @@ testfile_radiation_csv = DATA_DIR / \ 'pvgis_hourly_Timeseries_45.000_8.000_SA_30deg_0deg_2016_2016.csv' testfile_pv_json = DATA_DIR / \ - 'pvgis_hourly_Timeseries_45.000_8.000_CM_10kWp_CIS_5_2a_2013_2014.json' + 'pvgis_hourly_Timeseries_45.000_8.000_SA2_10kWp_CIS_5_2a_2013_2014.json' index_radiation_csv = \ pd.date_range('20160101 00:10', freq='1h', periods=14, tz='UTC') index_pv_json = \ - pd.date_range('2013-01-01 00:55', freq='1h', periods=10, tz='UTC') + pd.date_range('2013-01-01 00:10', freq='1h', periods=10, tz='UTC') columns_radiation_csv = [ 'Gb(i)', 'Gd(i)', 'Gr(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] @@ -34,10 +34,9 @@ 'poa_direct', 'poa_sky_diffuse', 'poa_ground_diffuse', 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] columns_pv_json = [ - 'P', 'Gb(i)', 'Gd(i)', 'Gr(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] + 'P', 'G(i)', 'H_sun', 'T2m', 'WS10m', 'Int'] columns_pv_json_mapped = [ - 'P', 'poa_direct', 'poa_sky_diffuse', 'poa_ground_diffuse', - 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] + 'P', 'poa_global', 'solar_elevation', 'temp_air', 'wind_speed', 'Int'] data_radiation_csv = [ [0.0, 0.0, 0.0, 0.0, 3.44, 1.43, 0.0], @@ -55,20 +54,20 @@ [4.25, 1.88, 0.05, 21.41, 7.84, 0.4, 1.0], [0.0, 0.0, 0.0, 0.0, 7.43, 0.72, 0.0]] data_pv_json = [ - [0.0, 0.0, 0.0, 0.0, 0.0, 3.01, 1.23, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 2.22, 1.46, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 1.43, 1.7, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.64, 1.93, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.77, 1.8, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.91, 1.66, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 1.05, 1.53, 0.0], - [3464.5, 270.35, 91.27, 6.09, 6.12, 1.92, 1.44, 0.0], - [1586.9, 80.76, 83.95, 9.04, 13.28, 2.79, 1.36, 0.0], - [713.3, 5.18, 70.57, 7.31, 18.56, 3.66, 1.27, 0.0]] - -data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, - 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, - 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, + [0.0, 0.0, 0.0, -0.97, 1.52, 0.0], + [0.0, 0.0, 0.0, -1.06, 1.45, 0.0], + [0.0, 0.0, 0.0, -1.03, 1.45, 0.0], + [0.0, 0.0, 0.0, -0.48, 1.31, 0.0], + [0.0, 0.0, 0.0, -0.09, 1.24, 0.0], + [0.0, 0.0, 0.0, -0.38, 1.17, 0.0], + [0.0, 0.0, 0.0, 0.29, 1.03, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.62, 0.0], + [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], + [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] + +data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, + 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, + 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -90,7 +89,7 @@ inputs_pv_json = { 'location': {'latitude': 45.0, 'longitude': 8.0, 'elevation': 250.0}, - 'meteo_data': {'radiation_db': 'PVGIS-CMSAF', 'meteo_db': 'ERA-Interim', + 'meteo_data': {'radiation_db': 'PVGIS-SARAH2', 'meteo_db': 'ERA-Interim', 'year_min': 2013, 'year_max': 2014, 'use_horizon': True, 'horizon_db': None, 'horizon_data': 'DEM-calculated'}, 'mounting_system': {'two_axis': { @@ -98,45 +97,45 @@ 'azimuth': {'value': '-', 'optimal': '-'}}}, 'pv_module': {'technology': 'CIS', 'peak_power': 10.0, 'system_loss': 5.0}} + metadata_pv_json = { 'inputs': { - 'location': {'description': 'Selected location', 'variables': { - 'latitude': {'description': 'Latitude', 'units': 'decimal degree'}, - 'longitude': {'description': 'Longitude', 'units': 'decimal degree'}, # noqa: E501 - 'elevation': {'description': 'Elevation', 'units': 'm'}}}, - 'meteo_data': { - 'description': 'Sources of meteorological data', - 'variables': { - 'radiation_db': {'description': 'Solar radiation database'}, - 'meteo_db': {'description': 'Database used for meteorological variables other than solar radiation'}, # noqa: E501 - 'year_min': {'description': 'First year of the calculations'}, - 'year_max': {'description': 'Last year of the calculations'}, - 'use_horizon': {'description': 'Include horizon shadows'}, - 'horizon_db': {'description': 'Source of horizon data'}}}, - 'mounting_system': { - 'description': 'Mounting system', - 'choices': 'fixed, vertical_axis, inclined_axis, two_axis', - 'fields': { - 'slope': {'description': 'Inclination angle from the horizontal plane', 'units': 'degree'}, # noqa: E501 - 'azimuth': {'description': 'Orientation (azimuth) angle of the (fixed) PV system (0 = S, 90 = W, -90 = E)', 'units': 'degree'}}}, # noqa: E501 - 'pv_module': { - 'description': 'PV module parameters', - 'variables': { - 'technology': {'description': 'PV technology'}, - 'peak_power': {'description': 'Nominal (peak) power of the PV module', 'units': 'kW'}, # noqa: E501 - 'system_loss': {'description': 'Sum of system losses', 'units': '%'}}}}, # noqa: E501 - 'outputs': { - 'hourly': { - 'type': 'time series', 'timestamp': 'hourly averages', - 'variables': { - 'P': {'description': 'PV system power', 'units': 'W'}, - 'Gb(i)': {'description': 'Beam (direct) irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 - 'Gd(i)': {'description': 'Diffuse irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 - 'Gr(i)': {'description': 'Reflected irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 - 'H_sun': {'description': 'Sun height', 'units': 'degree'}, - 'T2m': {'description': '2-m air temperature', 'units': 'degree Celsius'}, # noqa: E501 - 'WS10m': {'description': '10-m total wind speed', 'units': 'm/s'}, # noqa: E501 - 'Int': {'description': '1 means solar radiation values are reconstructed'}}}}} # noqa: E501 + 'location': + {'description': 'Selected location', 'variables': { + 'latitude': {'description': 'Latitude', 'units': 'decimal degree'}, # noqa: E501 + 'longitude': {'description': 'Longitude', 'units': 'decimal degree'}, # noqa: E501 + 'elevation': {'description': 'Elevation', 'units': 'm'}}}, + 'meteo_data': { + 'description': 'Sources of meteorological data', + 'variables': { + 'radiation_db': {'description': 'Solar radiation database'}, # noqa: E501 + 'meteo_db': {'description': 'Database used for meteorological variables other than solar radiation'}, # noqa: E501 + 'year_min': {'description': 'First year of the calculations'}, # noqa: E501 + 'year_max': {'description': 'Last year of the calculations'}, # noqa: E501 + 'use_horizon': {'description': 'Include horizon shadows'}, + 'horizon_db': {'description': 'Source of horizon data'}}}, + 'mounting_system': { + 'description': 'Mounting system', + 'choices': 'fixed, vertical_axis, inclined_axis, two_axis', + 'fields': { + 'slope': {'description': 'Inclination angle from the horizontal plane', 'units': 'degree'}, # noqa: E501 + 'azimuth': {'description': 'Orientation (azimuth) angle of the (fixed) PV system (0 = S, 90 = W, -90 = E)', 'units': 'degree'}}}, # noqa: E501 + 'pv_module': { + 'description': 'PV module parameters', + 'variables': { + 'technology': {'description': 'PV technology'}, + 'peak_power': {'description': 'Nominal (peak) power of the PV module', 'units': 'kW'}, # noqa: E501 + 'system_loss': {'description': 'Sum of system losses', 'units': '%'}}}}, # noqa: E501 + 'outputs': { + 'hourly': { + 'type': 'time series', 'timestamp': 'hourly averages', + 'variables': { + 'P': {'description': 'PV system power', 'units': 'W'}, + 'G(i)': {'description': 'Global irradiance on the inclined plane (plane of the array)', 'units': 'W/m2'}, # noqa: E501 + 'H_sun': {'description': 'Sun height', 'units': 'degree'}, + 'T2m': {'description': '2-m air temperature', 'units': 'degree Celsius'}, # noqa: E501 + 'WS10m': {'description': '10-m total wind speed', 'units': 'm/s'}, # noqa: E501 + 'Int': {'description': '1 means solar radiation values are reconstructed'}}}}} # noqa: E501 def generate_expected_dataframe(values, columns, index): @@ -225,12 +224,13 @@ def test_read_pvgis_hourly_bad_extension(): args_pv_json = { 'surface_tilt': 30, 'surface_azimuth': 0, 'outputformat': 'json', - 'usehorizon': True, 'userhorizon': None, 'raddatabase': 'PVGIS-CMSAF', + 'usehorizon': True, 'userhorizon': None, 'raddatabase': 'PVGIS-SARAH2', 'start': pd.Timestamp(2013, 1, 1), 'end': pd.Timestamp(2014, 5, 1), 'pvcalculation': True, 'peakpower': 10, 'pvtechchoice': 'CIS', 'loss': 5, - 'trackingtype': 2, 'optimalangles': True, 'components': True} + 'trackingtype': 2, 'optimalangles': True, 'components': False, + 'url': 'https://re.jrc.ec.europa.eu/api/v5_2/'} -url_pv_json = 'https://re.jrc.ec.europa.eu/api/seriescalc?lat=45&lon=8&outputformat=json&angle=30&aspect=0&pvtechchoice=CIS&mountingplace=free&trackingtype=2&components=1&usehorizon=1&raddatabase=PVGIS-CMSAF&startyear=2013&endyear=2014&pvcalculation=1&peakpower=10&loss=5&optimalangles=1' # noqa: E501 +url_pv_json = 'https://re.jrc.ec.europa.eu/api/v5_2/seriescalc?lat=45&lon=8&outputformat=json&angle=30&aspect=0&pvtechchoice=CIS&mountingplace=free&trackingtype=2&components=0&usehorizon=1&raddatabase=PVGIS-SARAH2&startyear=2013&endyear=2014&pvcalculation=1&peakpower=10&loss=5&optimalangles=1' # noqa: E501 @pytest.mark.parametrize('testfile,expected_name,args,map_variables,url_test', [ # noqa: E501 @@ -518,13 +518,11 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): actual, _, _, _ = get_pvgis_tmy(45, 8, map_variables=True) assert all([c in pvgis_tmy_mapped_columns for c in actual.columns]) - @pytest.mark.remote_data def test_read_pvgis_horizon(): azimuth, horizon = get_pvgis_horizon(35.171051, -106.465158) assert all(np.isclose(horizon, data_horizon_abq)) - - + def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' actual, _, _, _ = read_pvgis_tmy(fn, map_variables=True) @@ -604,3 +602,4 @@ def test_read_pvgis_tmy_exception(): with pytest.raises(ValueError, match=err_msg): read_pvgis_tmy('filename', pvgis_format=bad_outputformat, map_variables=False) + \ No newline at end of file From 7b99b8d03e50048decf7908585cb798fc5b5a9c4 Mon Sep 17 00:00:00 2001 From: Pierce Date: Thu, 13 Oct 2022 16:33:28 -0600 Subject: [PATCH 43/94] fix formatting --- pvlib/horizon.py | 4 +++- pvlib/iotools/cgiar.py | 8 ++++---- pvlib/iotools/pvgis.py | 10 ++++++---- pvlib/tests/iotools/test_pvgis.py | 9 +++++---- 4 files changed, 18 insertions(+), 13 deletions(-) diff --git a/pvlib/horizon.py b/pvlib/horizon.py index 089b00bc86..8558a6d7a8 100644 --- a/pvlib/horizon.py +++ b/pvlib/horizon.py @@ -6,7 +6,7 @@ gdal, osgeo, geoio, and scikit-image. ''' import numpy as np -from pathlib import Path + def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): ''' @@ -46,6 +46,7 @@ def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): adjustment[mask] = 0 return adjustment + def calculate_dtf(horizon_azimuths, horizon_angles, surface_tilt, surface_azimuth): """ @@ -115,6 +116,7 @@ def calculate_dtf(horizon_azimuths, horizon_angles, dtf += 2 * (first_term + second_term) / num_points return dtf + def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): """ Determine the elevation angle created by the surface of a tilted plane diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 5de4b36e21..3701f85682 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -75,14 +75,14 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, org_url = 'https://srtm.csi.cgiar.org/' base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' - res = requests.get(query_URL, proxies = proxies, verify = False) + res = requests.get(query_URL, proxies=proxies, verify=False) res.raise_for_status() - + # unzip download zipfile = ZipFile(io.BytesIO(res.content)) ext = '.tif' files = zipfile.namelist() - file = [ f for f in files if ext in f ][0] - path = zipfile.extract(file, path = path_to_save) + file = [f for f in files if ext in f][0] + path = zipfile.extract(file, path=path_to_save) img = skimage.io.imread(path) return img, path diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index c4c9575c3f..50608d701f 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -695,11 +695,13 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): # the horizon data is given in a different format then the others # numpy has an easier time decoding it array = np.genfromtxt(io.StringIO(string), - skip_header = 4, skip_footer = 7) + skip_header=4, skip_footer=7) df = pd.DataFrame(array) - + # Set the column names df.columns = ['horizon_azimuth', 'horizon_angles', - 'azimuth_sun_winter_solstice', 'elevation_sun_winter_solstice', - 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] + 'azimuth_sun_winter_solstice', + 'elevation_sun_winter_solstice', + 'azimuth_sun_summer_solstice', + 'elevation_sun_summer_solstice'] return df diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 7b2565fe69..59d07dd652 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -65,9 +65,9 @@ [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] -data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, - 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, - 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, +data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, + 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, + 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -518,11 +518,12 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): actual, _, _, _ = get_pvgis_tmy(45, 8, map_variables=True) assert all([c in pvgis_tmy_mapped_columns for c in actual.columns]) + @pytest.mark.remote_data def test_read_pvgis_horizon(): azimuth, horizon = get_pvgis_horizon(35.171051, -106.465158) assert all(np.isclose(horizon, data_horizon_abq)) - + def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' actual, _, _, _ = read_pvgis_tmy(fn, map_variables=True) From 9ea61e2a1643ef9d56ea4a523939bc2a696796a6 Mon Sep 17 00:00:00 2001 From: Pierce Date: Thu, 13 Oct 2022 16:34:32 -0600 Subject: [PATCH 44/94] fix formatting, 2 --- pvlib/tests/iotools/test_pvgis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 59d07dd652..8c0026e99a 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -603,4 +603,3 @@ def test_read_pvgis_tmy_exception(): with pytest.raises(ValueError, match=err_msg): read_pvgis_tmy('filename', pvgis_format=bad_outputformat, map_variables=False) - \ No newline at end of file From e23d58c06ad2943c3f9c9d1f234cccde383f7bfb Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 31 Oct 2022 09:49:42 -0700 Subject: [PATCH 45/94] removing horizon.py --- pvlib/horizon.py | 155 ----------------------------------------------- 1 file changed, 155 deletions(-) delete mode 100644 pvlib/horizon.py diff --git a/pvlib/horizon.py b/pvlib/horizon.py deleted file mode 100644 index 8558a6d7a8..0000000000 --- a/pvlib/horizon.py +++ /dev/null @@ -1,155 +0,0 @@ -''' -The 'horizon' module contains function definitions that -retrive & calculate the surrounding horizon using DEM -elevation data and its impact on performance. -Optional dependencies for this module include -gdal, osgeo, geoio, and scikit-image. -''' -import numpy as np - - -def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): - ''' - 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 - ---------- - 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] != 360): - 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 - horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] - mask = solar_zenith > horizon_zenith - adjustment[mask] = 0 - return adjustment - - -def calculate_dtf(horizon_azimuths, horizon_angles, - surface_tilt, surface_azimuth): - """ - 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 - ---------- - 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 - 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 - 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. - 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. - A similar method is used in section II of [1] although it is looking at - both ground and sky diffuse irradiation. - [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 - """ - 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) - b = np.sin(tilt_rad) * np.sin(plane_az_rad) - c = np.cos(tilt_rad) - - # this gets either a float or an array of zeros - dtf = np.multiply(0.0, surface_tilt) - num_points = horizon_azimuths.shape[0] - for i in range(horizon_azimuths.shape[0]): - az = np.radians(horizon_azimuths[i]) - 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 - 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 - - -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. This comes from Equation 10 in [1] - 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 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 - -------- - 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. - [1] doi.org/10.1016/j.solener.2014.09.037 - """ - 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) - elevation_angle = 90.0 - declination - elevation_angle[mask] = 0.0 - - return elevation_angle From 801f4d31cece233d1b033db4b6060641f11e5b8d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 31 Oct 2022 13:16:09 -0700 Subject: [PATCH 46/94] tests added --- pvlib/iotools/__init__.py | 2 ++ pvlib/tests/iotools/test_cgiar.py | 18 ++++++++++++++++++ setup.py | 2 +- 3 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 pvlib/tests/iotools/test_cgiar.py diff --git a/pvlib/iotools/__init__.py b/pvlib/iotools/__init__.py index b02ce243ae..4cca756283 100644 --- a/pvlib/iotools/__init__.py +++ b/pvlib/iotools/__init__.py @@ -14,6 +14,7 @@ from pvlib.iotools.psm3 import parse_psm3 # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_tmy, read_pvgis_tmy # noqa: F401 from pvlib.iotools.pvgis import read_pvgis_hourly # noqa: F401 +from pvlib.iotools.pvgis import get_pvgis_horizon from pvlib.iotools.pvgis import get_pvgis_hourly # noqa: F401 from pvlib.iotools.bsrn import get_bsrn # noqa: F401 from pvlib.iotools.bsrn import read_bsrn # noqa: F401 @@ -21,3 +22,4 @@ from pvlib.iotools.sodapro import get_cams # noqa: F401 from pvlib.iotools.sodapro import read_cams # noqa: F401 from pvlib.iotools.sodapro import parse_cams # noqa: F401 +from pvlib.iotools.cgiar import download_SRTM diff --git a/pvlib/tests/iotools/test_cgiar.py b/pvlib/tests/iotools/test_cgiar.py new file mode 100644 index 0000000000..c025798df4 --- /dev/null +++ b/pvlib/tests/iotools/test_cgiar.py @@ -0,0 +1,18 @@ +import pytest +import requests +import numpy as np +from pvlib.iotools import download_SRTM +from ..conftest import (DATA_DIR, RERUNS, RERUNS_DELAY, assert_frame_equal, + fail_on_pvlib_version) +from pvlib._deprecation import pvlibDeprecationWarning + +NM_lat, NM_lon = 35, -107 +NM_elevations = [1418, 2647, 3351] +NM_points = ((1000, 2000, 3000), + (1000, 2000, 3000)) + + +@pytest.mark.remote_data +def test_cgiar_download(): + DEM, file_path = download_SRTM(35, -107, path_to_save=DATA_DIR) + assert np.allclose(DEM[NM_points], NM_elevations) diff --git a/setup.py b/setup.py index 134ec1d88c..bddc953731 100755 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ EXTRAS_REQUIRE = { 'optional': ['cython', 'ephem', 'netcdf4', 'nrel-pysam', 'numba', 'pvfactors', 'siphon', 'statsmodels', - 'cftime >= 1.1.1'], + 'cftime >= 1.1.1', 'scikit-image'], 'doc': ['ipython', 'matplotlib', 'sphinx == 4.5.0', 'pydata-sphinx-theme == 0.8.1', 'sphinx-gallery', 'docutils == 0.15.2', 'pillow', 'netcdf4', 'siphon', From bcca59b311821d7e8be181dcb5430f738dbe4151 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 24 Apr 2023 09:03:54 -0700 Subject: [PATCH 47/94] updated test --- pvlib/tests/iotools/test_pvgis.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 8c0026e99a..7ccfe061ff 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -521,8 +521,9 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data def test_read_pvgis_horizon(): - azimuth, horizon = get_pvgis_horizon(35.171051, -106.465158) - assert all(np.isclose(horizon, data_horizon_abq)) + df = get_pvgis_horizon(35.171051, -106.465158) + az, elv = df.horizon_azimuth, df.horizon_angles + assert all(np.isclose(elv, data_horizon_abq)) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' From 6cfb0c49ccc3b934be2496bf2a66e1fef204cceb Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 24 Apr 2023 09:32:09 -0700 Subject: [PATCH 48/94] fixed stickler --- pvlib/iotools/__init__.py | 4 ++-- pvlib/tests/iotools/test_cgiar.py | 5 ++--- pvlib/tests/iotools/test_pvgis.py | 2 +- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/pvlib/iotools/__init__.py b/pvlib/iotools/__init__.py index 4cca756283..c1b1372faa 100644 --- a/pvlib/iotools/__init__.py +++ b/pvlib/iotools/__init__.py @@ -14,7 +14,7 @@ from pvlib.iotools.psm3 import parse_psm3 # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_tmy, read_pvgis_tmy # noqa: F401 from pvlib.iotools.pvgis import read_pvgis_hourly # noqa: F401 -from pvlib.iotools.pvgis import get_pvgis_horizon +from pvlib.iotools.pvgis import get_pvgis_horizon # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_hourly # noqa: F401 from pvlib.iotools.bsrn import get_bsrn # noqa: F401 from pvlib.iotools.bsrn import read_bsrn # noqa: F401 @@ -22,4 +22,4 @@ from pvlib.iotools.sodapro import get_cams # noqa: F401 from pvlib.iotools.sodapro import read_cams # noqa: F401 from pvlib.iotools.sodapro import parse_cams # noqa: F401 -from pvlib.iotools.cgiar import download_SRTM +from pvlib.iotools.cgiar import download_SRTM # noqa: F401 diff --git a/pvlib/tests/iotools/test_cgiar.py b/pvlib/tests/iotools/test_cgiar.py index c025798df4..8ee1ef7eca 100644 --- a/pvlib/tests/iotools/test_cgiar.py +++ b/pvlib/tests/iotools/test_cgiar.py @@ -1,9 +1,8 @@ import pytest -import requests import numpy as np from pvlib.iotools import download_SRTM -from ..conftest import (DATA_DIR, RERUNS, RERUNS_DELAY, assert_frame_equal, - fail_on_pvlib_version) +from ..conftest import DATA_DIR + from pvlib._deprecation import pvlibDeprecationWarning NM_lat, NM_lon = 35, -107 diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 7ccfe061ff..29a9fd4ecf 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -522,7 +522,7 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data def test_read_pvgis_horizon(): df = get_pvgis_horizon(35.171051, -106.465158) - az, elv = df.horizon_azimuth, df.horizon_angles + elv = df.horizon_angles assert all(np.isclose(elv, data_horizon_abq)) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): From d5cd1aff93ed62976c50bf5a56de8a77d9fc6a8b Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 24 Apr 2023 09:45:39 -0700 Subject: [PATCH 49/94] minor formatting/stickler --- pvlib/iotools/__init__.py | 4 ++-- pvlib/tests/iotools/test_cgiar.py | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/pvlib/iotools/__init__.py b/pvlib/iotools/__init__.py index c1b1372faa..7fc57bfd2c 100644 --- a/pvlib/iotools/__init__.py +++ b/pvlib/iotools/__init__.py @@ -14,7 +14,7 @@ from pvlib.iotools.psm3 import parse_psm3 # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_tmy, read_pvgis_tmy # noqa: F401 from pvlib.iotools.pvgis import read_pvgis_hourly # noqa: F401 -from pvlib.iotools.pvgis import get_pvgis_horizon # noqa: F401 +from pvlib.iotools.pvgis import get_pvgis_horizon # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_hourly # noqa: F401 from pvlib.iotools.bsrn import get_bsrn # noqa: F401 from pvlib.iotools.bsrn import read_bsrn # noqa: F401 @@ -22,4 +22,4 @@ from pvlib.iotools.sodapro import get_cams # noqa: F401 from pvlib.iotools.sodapro import read_cams # noqa: F401 from pvlib.iotools.sodapro import parse_cams # noqa: F401 -from pvlib.iotools.cgiar import download_SRTM # noqa: F401 +from pvlib.iotools.cgiar import download_SRTM # noqa: F401 diff --git a/pvlib/tests/iotools/test_cgiar.py b/pvlib/tests/iotools/test_cgiar.py index 8ee1ef7eca..3943c7807d 100644 --- a/pvlib/tests/iotools/test_cgiar.py +++ b/pvlib/tests/iotools/test_cgiar.py @@ -3,7 +3,6 @@ from pvlib.iotools import download_SRTM from ..conftest import DATA_DIR -from pvlib._deprecation import pvlibDeprecationWarning NM_lat, NM_lon = 35, -107 NM_elevations = [1418, 2647, 3351] From 53d9e6c052c32331006d0455c648f41cd08e1051 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 24 Apr 2023 11:19:14 -0600 Subject: [PATCH 50/94] Update v0.9.6.rst --- docs/sphinx/source/whatsnew/v0.9.6.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.6.rst b/docs/sphinx/source/whatsnew/v0.9.6.rst index f21fe80ac6..6c01bc46a0 100644 --- a/docs/sphinx/source/whatsnew/v0.9.6.rst +++ b/docs/sphinx/source/whatsnew/v0.9.6.rst @@ -11,7 +11,10 @@ Deprecations Enhancements ~~~~~~~~~~~~ - +* Added horizon retrieval from PVGIS function in + :py:func:`pvlib.iotools.pvgis.get_pvgis_horizon`. (:issue:`1290`, :pull:`1395`) +* Added SRTM Digital Elevation Model retrieval from CGIAR in + :py:func:`pvlib.iotools.cgiar.download_SRTM`. (:issue:`1290`, :pull:`1395`) Bug fixes ~~~~~~~~~ From c8311072966d968062b106a7f17723411a403789 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 24 Apr 2023 12:54:49 -0600 Subject: [PATCH 51/94] Update iotools.rst --- docs/sphinx/source/reference/iotools.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/sphinx/source/reference/iotools.rst b/docs/sphinx/source/reference/iotools.rst index 514aeac2f5..a63bf52284 100644 --- a/docs/sphinx/source/reference/iotools.rst +++ b/docs/sphinx/source/reference/iotools.rst @@ -37,6 +37,8 @@ of sources and file formats relevant to solar energy modeling. iotools.get_cams iotools.read_cams iotools.parse_cams + iotools.get_pvgis_horizon + iotools.download_SRTM A :py:class:`~pvlib.location.Location` object may be created from metadata in some files. From 17e4f92d7d43e83c1a4ace38c996e69766740ada Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 14:47:48 -0600 Subject: [PATCH 52/94] Apply suggestions from code review Co-authored-by: Kevin Anderson --- docs/sphinx/source/whatsnew/v0.9.6.rst | 4 ++-- pvlib/iotools/cgiar.py | 4 +++- pvlib/iotools/pvgis.py | 1 + 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.6.rst b/docs/sphinx/source/whatsnew/v0.9.6.rst index 6c01bc46a0..c8156ec0ba 100644 --- a/docs/sphinx/source/whatsnew/v0.9.6.rst +++ b/docs/sphinx/source/whatsnew/v0.9.6.rst @@ -12,9 +12,9 @@ Deprecations Enhancements ~~~~~~~~~~~~ * Added horizon retrieval from PVGIS function in - :py:func:`pvlib.iotools.pvgis.get_pvgis_horizon`. (:issue:`1290`, :pull:`1395`) + :py:func:`pvlib.iotools.get_pvgis_horizon`. (:issue:`1290`, :pull:`1395`) * Added SRTM Digital Elevation Model retrieval from CGIAR in - :py:func:`pvlib.iotools.cgiar.download_SRTM`. (:issue:`1290`, :pull:`1395`) + :py:func:`pvlib.iotools.download_SRTM`. (:issue:`1290`, :pull:`1395`) Bug fixes ~~~~~~~~~ diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 3701f85682..526699286b 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -47,13 +47,15 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, path_to_save="./", proxies=None): r'''Downloads a SRTM DEM tile from CGIAR, saves it to a path, and loads it as an array + + Parameters ---------- latitude: numeric latitude value to be included in the DEM longitude : numeric longitude value to be included in the DEM, negative W of prime meridian - srtm_arc_sec: numeric, {1,3} + srtm_arc_sec : int, {3, 1} The resolution (arc-seconds) of the desired DEM. Either 1 or 3 path_to_save: string diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 50608d701f..5251c06114 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -671,6 +671,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): r''' Get horizon data from PVGIS + Parameters ---------- latitude : float From c176d9bdee731b4b8b79f849f8d25d86ae745428 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 14:51:58 -0600 Subject: [PATCH 53/94] Change assert to use np builtin --- pvlib/tests/iotools/test_cgiar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/tests/iotools/test_cgiar.py b/pvlib/tests/iotools/test_cgiar.py index 3943c7807d..b78d65ffcc 100644 --- a/pvlib/tests/iotools/test_cgiar.py +++ b/pvlib/tests/iotools/test_cgiar.py @@ -13,4 +13,4 @@ @pytest.mark.remote_data def test_cgiar_download(): DEM, file_path = download_SRTM(35, -107, path_to_save=DATA_DIR) - assert np.allclose(DEM[NM_points], NM_elevations) + np.testing.assert_allclose(DEM[NM_points], NM_elevations) From 0d2a061f9aa98fb95b03be713676d9b4c7291aa1 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 14:57:39 -0600 Subject: [PATCH 54/94] correct docstring for hidden function --- pvlib/iotools/cgiar.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index 526699286b..a301c6f19a 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -17,19 +17,19 @@ def _lat_lon_to_query(longitude, latitude, srtm_arc_sec=3): r'''Converts latitude, longitude from the format used as input to the other functions to format used by CGIAR ---------- - longitude : numeric + longitude : float longitude value, negative W of prime meridian - latitude: numeric + latitude: float latitude value - srtm_arc_sec: numeric {1,3} + srtm_arc_sec: int, {3,1} The resolution (arc-seconds) of the desired DEM. Either 1 or 3 Returns ------- - rounded_lon : numeric + rounded_lon : int Rounded/adjusted longitude value - rounded_lat : numeric + rounded_lat : int Rounded/adjusted latitude value ''' From 6980a4d8e6a621522cf3b5b4d17ad92fb129752c Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 15:02:21 -0600 Subject: [PATCH 55/94] Chnage proxy settings to passing through kwargs --- pvlib/iotools/cgiar.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py index a301c6f19a..94e04c6859 100644 --- a/pvlib/iotools/cgiar.py +++ b/pvlib/iotools/cgiar.py @@ -44,7 +44,7 @@ def _lat_lon_to_query(longitude, latitude, srtm_arc_sec=3): def download_SRTM(latitude, longitude, srtm_arc_sec=3, - path_to_save="./", proxies=None): + path_to_save="./", **kwargs): r'''Downloads a SRTM DEM tile from CGIAR, saves it to a path, and loads it as an array @@ -62,6 +62,8 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, The base path to save the DEM as a .tif file proxies: dict Proxy table for a corporate network + kwargs: + Passed to requests.get Returns ------- @@ -77,7 +79,7 @@ def download_SRTM(latitude, longitude, srtm_arc_sec=3, org_url = 'https://srtm.csi.cgiar.org/' base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' - res = requests.get(query_URL, proxies=proxies, verify=False) + res = requests.get(query_URL, **kwargs) res.raise_for_status() # unzip download From cf2197496e87fa4a67d0fa142458835ca2dea07a Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 15:04:52 -0600 Subject: [PATCH 56/94] Changed proxy settings to pass kwargs --- pvlib/iotools/pvgis.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 5251c06114..616d1cc6b9 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -668,7 +668,7 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta -def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): +def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): r''' Get horizon data from PVGIS @@ -678,10 +678,10 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): Latitude in degrees north longitude : float Longitude in degrees east - proxies: dict - Dictionary of proxies to access through a corporate network url: string Base URL for PVGIS + kwargs: + Passed to requests.get Returns ------- @@ -690,7 +690,7 @@ def get_pvgis_horizon(latitude, longitude, proxies=None, url=URL): ''' res = requests.get(url + f'printhorizon?lat={latitude}&lon={longitude}', - proxies=proxies, verify=False) + **kwargs) res.raise_for_status() string = str(io.BytesIO(res.content).read().decode('UTF-8')) # the horizon data is given in a different format then the others From 0fb9e05751f496236f2542c7987eee41baf89d2d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 15:15:38 -0600 Subject: [PATCH 57/94] Change assertion to use assert_series_equal from conftest --- pvlib/tests/iotools/test_pvgis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 29a9fd4ecf..a1993fa9e6 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -11,7 +11,7 @@ from pvlib.iotools import get_pvgis_hourly, read_pvgis_hourly from pvlib.iotools import get_pvgis_horizon from ..conftest import (DATA_DIR, RERUNS, RERUNS_DELAY, assert_frame_equal, - fail_on_pvlib_version) + fail_on_pvlib_version, assert_series_equal) from pvlib._deprecation import pvlibDeprecationWarning @@ -523,7 +523,7 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): def test_read_pvgis_horizon(): df = get_pvgis_horizon(35.171051, -106.465158) elv = df.horizon_angles - assert all(np.isclose(elv, data_horizon_abq)) + assert_series_equal(elv, data_horizon_abq) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' From 7736e35fb0d009ab5363c0d5d912194006b7a659 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 25 Apr 2023 15:25:28 -0600 Subject: [PATCH 58/94] Added horizon profile contributions --- docs/sphinx/source/whatsnew/v0.9.6.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/sphinx/source/whatsnew/v0.9.6.rst b/docs/sphinx/source/whatsnew/v0.9.6.rst index c8156ec0ba..df4d9bca6d 100644 --- a/docs/sphinx/source/whatsnew/v0.9.6.rst +++ b/docs/sphinx/source/whatsnew/v0.9.6.rst @@ -40,3 +40,12 @@ Requirements Contributors ~~~~~~~~~~~~ * Adam R. Jensen (:ghuser:`adamrjensen`) +* Ben Pierce (:ghuser:`bgpierc`) +* Joseph Palakapilly (:ghuser:`JPalakapillyKWH`) +* Cliff Hansen (:ghuser:`cwhanse`) +* Anton Driesse (:ghuser:`adriesse`) +* Will Holmgren (:ghuser:`wholmgren`) +* Mark Mikofski (:ghuser:`mikofski`) +* Karel De Brabandere (:ghuser:`kdebrab`) +* Josh Stein (:ghuser:`jsstein`) +* Kevin Anderson (:ghuser:`kandersolar`) From 7f8453156bd0cb096fd513669353991b43ca814a Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:21:23 -0600 Subject: [PATCH 59/94] Update v0.9.6.rst --- docs/sphinx/source/whatsnew/v0.9.6.rst | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.6.rst b/docs/sphinx/source/whatsnew/v0.9.6.rst index df4d9bca6d..06359e3959 100644 --- a/docs/sphinx/source/whatsnew/v0.9.6.rst +++ b/docs/sphinx/source/whatsnew/v0.9.6.rst @@ -13,8 +13,6 @@ Enhancements ~~~~~~~~~~~~ * Added horizon retrieval from PVGIS function in :py:func:`pvlib.iotools.get_pvgis_horizon`. (:issue:`1290`, :pull:`1395`) -* Added SRTM Digital Elevation Model retrieval from CGIAR in - :py:func:`pvlib.iotools.download_SRTM`. (:issue:`1290`, :pull:`1395`) Bug fixes ~~~~~~~~~ From cb7eb41b1e4acb37a5ff8032fa34751101f648c2 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:22:02 -0600 Subject: [PATCH 60/94] Delete cgiar.py --- pvlib/iotools/cgiar.py | 92 ------------------------------------------ 1 file changed, 92 deletions(-) delete mode 100644 pvlib/iotools/cgiar.py diff --git a/pvlib/iotools/cgiar.py b/pvlib/iotools/cgiar.py deleted file mode 100644 index 94e04c6859..0000000000 --- a/pvlib/iotools/cgiar.py +++ /dev/null @@ -1,92 +0,0 @@ -""" -Get, read, and parse data from -`CGIAR GeoPortal `_. - -For more information, see the following links: -* `FAQ `_ - - -""" -import math -from zipfile import ZipFile -import requests -import io - - -def _lat_lon_to_query(longitude, latitude, srtm_arc_sec=3): - r'''Converts latitude, longitude from the format used as - input to the other functions to format used by CGIAR - ---------- - longitude : float - longitude value, negative W of prime meridian - latitude: float - latitude value - srtm_arc_sec: int, {3,1} - The resolution (arc-seconds) of the desired DEM. - Either 1 or 3 - - Returns - ------- - rounded_lon : int - Rounded/adjusted longitude value - rounded_lat : int - Rounded/adjusted latitude value - - ''' - rounded = (int(math.floor(longitude)), int(math.floor(latitude))) - if srtm_arc_sec == 1: - return rounded - elif srtm_arc_sec == 3: - rounded_lon, rounded_lat = rounded - return (rounded_lon + 180) // 5 + 1, (64 - rounded_lat) // 5 - else: - raise Exception("Please use SRTM 1 Arc Second or SRTM 3 Arc Second") - - -def download_SRTM(latitude, longitude, srtm_arc_sec=3, - path_to_save="./", **kwargs): - r'''Downloads a SRTM DEM tile from CGIAR, - saves it to a path, and loads it as an array - - Parameters - ---------- - latitude: numeric - latitude value to be included in the DEM - longitude : numeric - longitude value to be included in the DEM, - negative W of prime meridian - srtm_arc_sec : int, {3, 1} - The resolution (arc-seconds) of the desired DEM. - Either 1 or 3 - path_to_save: string - The base path to save the DEM as a .tif file - proxies: dict - Proxy table for a corporate network - kwargs: - Passed to requests.get - - Returns - ------- - img : np.array - Numpy array filled with elevation values in [m] - path : string - Path and filename where the DEM .tif filed was saved - - ''' - import skimage - # convert lat/lon to format cgiar uses - long, lat = _lat_lon_to_query(longitude, latitude, srtm_arc_sec) - org_url = 'https://srtm.csi.cgiar.org/' - base_url = org_url + 'wp-content/uploads/files/srtm_5x5/TIFF/' - query_URL = base_url + f'srtm_{long:02d}_{lat:02d}.zip' - res = requests.get(query_URL, **kwargs) - res.raise_for_status() - - # unzip download - zipfile = ZipFile(io.BytesIO(res.content)) - ext = '.tif' - files = zipfile.namelist() - file = [f for f in files if ext in f][0] - path = zipfile.extract(file, path=path_to_save) - img = skimage.io.imread(path) - return img, path From 25baf2bc492d290c69e1ad6bbb4ee8a266bcc91d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:22:30 -0600 Subject: [PATCH 61/94] Delete test_cgiar.py --- pvlib/tests/iotools/test_cgiar.py | 16 ---------------- 1 file changed, 16 deletions(-) delete mode 100644 pvlib/tests/iotools/test_cgiar.py diff --git a/pvlib/tests/iotools/test_cgiar.py b/pvlib/tests/iotools/test_cgiar.py deleted file mode 100644 index b78d65ffcc..0000000000 --- a/pvlib/tests/iotools/test_cgiar.py +++ /dev/null @@ -1,16 +0,0 @@ -import pytest -import numpy as np -from pvlib.iotools import download_SRTM -from ..conftest import DATA_DIR - - -NM_lat, NM_lon = 35, -107 -NM_elevations = [1418, 2647, 3351] -NM_points = ((1000, 2000, 3000), - (1000, 2000, 3000)) - - -@pytest.mark.remote_data -def test_cgiar_download(): - DEM, file_path = download_SRTM(35, -107, path_to_save=DATA_DIR) - np.testing.assert_allclose(DEM[NM_points], NM_elevations) From b9edf2504050c2741b4c234e5ffc7b4398290eb7 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:24:28 -0600 Subject: [PATCH 62/94] Update shading.py --- pvlib/shading.py | 161 ----------------------------------------------- 1 file changed, 161 deletions(-) diff --git a/pvlib/shading.py b/pvlib/shading.py index 7c00107a7d..01d725207c 100644 --- a/pvlib/shading.py +++ b/pvlib/shading.py @@ -191,164 +191,3 @@ def sky_diffuse_passias(masking_angle): Available at https://www.nrel.gov/docs/fy18osti/67399.pdf """ return 1 - cosd(masking_angle/2)**2 - - -def calculate_dtf(horizon_azimuths, horizon_angles, - surface_tilt, surface_azimuth): - r"""Author: JPalakapillyKWH - 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 - ---------- - 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 - 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 - 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. - - 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. - A similar method is used in section II of [1] although it is looking at - both ground and sky diffuse irradiation. - [2] Wright D. (2019) IEEE Journal of Photovoltaics 9(2), 391-396 - - This function was written by @JPalakapillyKWH - in an uncompleted pvlib-python pull request #758. - """ - 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) - b = np.sin(tilt_rad) * np.sin(plane_az_rad) - c = np.cos(tilt_rad) - - # this gets either a float or an array of zeros - dtf = np.multiply(0.0, surface_tilt) - num_points = horizon_azimuths.shape[0] - for i in range(horizon_azimuths.shape[0]): - az = np.radians(horizon_azimuths[i]) - 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 - 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 - - -def collection_plane_elev_angle(surface_tilt, surface_azimuth, direction): - r""" - 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. This comes from Equation 10 in [1] - 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 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 - -------- - 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. - [1] doi.org/10.1016/j.solener.2014.09.037 - - Notes - ------ - This function was written by @JPalakapillyKWH - in an uncompleted pvlib-python pull request #758. - """ - 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) - elevation_angle = 90.0 - declination - elevation_angle[mask] = 0.0 - - return elevation_angle - - -def dni_horizon_adjustment(horizon_angles, solar_zenith, solar_azimuth): - r''' - 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 - ---------- - 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. - - Notes - ------ - This function was written by @JPalakapillyKWH - in an uncompleted pvlib-python pull request #758. - ''' - adjustment = np.ones(solar_zenith.shape) - - if (horizon_angles.shape[0] != 360): - 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 - horizon_zenith = 90 - horizon_angles[rounded_solar_azimuth] - mask = solar_zenith > horizon_zenith - adjustment[mask] = 0 - return adjustment From d15c7ff93f82eb69a279c2a4bf06911b9c886672 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:25:34 -0600 Subject: [PATCH 63/94] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 1cef5be2b7..6e25d097da 100755 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ EXTRAS_REQUIRE = { 'optional': ['cython', 'ephem', 'netcdf4', 'nrel-pysam', 'numba', 'pvfactors', 'siphon', 'statsmodels', - 'cftime >= 1.1.1', 'scikit-image'], + 'cftime >= 1.1.1'], 'doc': ['ipython', 'matplotlib', 'sphinx == 4.5.0', 'pydata-sphinx-theme == 0.8.1', 'sphinx-gallery', 'docutils == 0.15.2', 'pillow', 'netcdf4', 'siphon', From 2bf5cb28224e3b43d7ae33f06c546a711699073a Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:26:49 -0600 Subject: [PATCH 64/94] Update __init__.py --- pvlib/iotools/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pvlib/iotools/__init__.py b/pvlib/iotools/__init__.py index 7fc57bfd2c..4358a208a5 100644 --- a/pvlib/iotools/__init__.py +++ b/pvlib/iotools/__init__.py @@ -22,4 +22,3 @@ from pvlib.iotools.sodapro import get_cams # noqa: F401 from pvlib.iotools.sodapro import read_cams # noqa: F401 from pvlib.iotools.sodapro import parse_cams # noqa: F401 -from pvlib.iotools.cgiar import download_SRTM # noqa: F401 From 89e53d92dce63f8255bcba65205304707e8be8f9 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 12:30:52 -0600 Subject: [PATCH 65/94] Update iotools.rst --- docs/sphinx/source/reference/iotools.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/sphinx/source/reference/iotools.rst b/docs/sphinx/source/reference/iotools.rst index a63bf52284..e6bd33907e 100644 --- a/docs/sphinx/source/reference/iotools.rst +++ b/docs/sphinx/source/reference/iotools.rst @@ -38,7 +38,6 @@ of sources and file formats relevant to solar energy modeling. iotools.read_cams iotools.parse_cams iotools.get_pvgis_horizon - iotools.download_SRTM A :py:class:`~pvlib.location.Location` object may be created from metadata in some files. From e13e322c0409ab92f7030f23f8a616d320a7bf03 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 14:14:58 -0600 Subject: [PATCH 66/94] Update test_pvgis.py --- pvlib/tests/iotools/test_pvgis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index a1993fa9e6..e457f3021b 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -65,14 +65,14 @@ [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] -data_horizon_abq = [9.9, 13.0, 14.5, 15.7, 14.9, 15.3, +data_horizon_abq = pd.Series([9.9, 13.0, 14.5, 15.7, 14.9, 15.3, 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1, 1.9, 3.8, 5.0, - 6.5, 9.2, 9.9] + 6.5, 9.2, 9.9], name='horizon_angles') inputs_radiation_csv = {'latitude': 45.0, 'longitude': 8.0, 'elevation': 250.0, 'radiation_database': 'PVGIS-SARAH', From 50e79cd884f26bc34ef4102454d8fe2332d9e5de Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 14:36:35 -0600 Subject: [PATCH 67/94] return only two columns and makde bytes->utf-8 implicit --- pvlib/iotools/pvgis.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 616d1cc6b9..349c642d06 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -692,10 +692,9 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): f'printhorizon?lat={latitude}&lon={longitude}', **kwargs) res.raise_for_status() - string = str(io.BytesIO(res.content).read().decode('UTF-8')) # the horizon data is given in a different format then the others # numpy has an easier time decoding it - array = np.genfromtxt(io.StringIO(string), + array = np.genfromtxt(io.StringIO(res.text), skip_header=4, skip_footer=7) df = pd.DataFrame(array) @@ -705,4 +704,4 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): 'elevation_sun_winter_solstice', 'azimuth_sun_summer_solstice', 'elevation_sun_summer_solstice'] - return df + return df[['horizon_azimuth', 'horizon_angles']] From 75a94077f58f9614632f6e7983b4aa6202507ec2 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 15:06:57 -0600 Subject: [PATCH 68/94] Use JSON to streamline processing --- pvlib/iotools/pvgis.py | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 349c642d06..8a6459b9cb 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -668,7 +668,8 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta -def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): +def get_pvgis_horizon(latitude, longitude, + url=URL, **kwargs): r''' Get horizon data from PVGIS @@ -678,7 +679,7 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): Latitude in degrees north longitude : float Longitude in degrees east - url: string + url: str, default Base URL for PVGIS kwargs: Passed to requests.get @@ -687,21 +688,24 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): ------- df : pd.DataFrame Pandas dataframe of the retrived horizon + + References + ---------- + .. [1] `PVGIS horizon profile tool + `_ ''' - res = requests.get(url + - f'printhorizon?lat={latitude}&lon={longitude}', + params = {'lat': latitude, 'lon': longitude, + 'outputformat': outputformat} + res = requests.get(url + 'printhorizon', params=params, **kwargs) - res.raise_for_status() - # the horizon data is given in a different format then the others - # numpy has an easier time decoding it - array = np.genfromtxt(io.StringIO(res.text), - skip_header=4, skip_footer=7) - df = pd.DataFrame(array) - - # Set the column names - df.columns = ['horizon_azimuth', 'horizon_angles', - 'azimuth_sun_winter_solstice', - 'elevation_sun_winter_solstice', - 'azimuth_sun_summer_solstice', - 'elevation_sun_summer_solstice'] - return df[['horizon_azimuth', 'horizon_angles']] + if not res.ok: + try: + err_msg = res.json() + except Exception: + res.raise_for_status() + else: + raise requests.HTTPError(err_msg['message']) + json_output = res.json() + df = pd.DataFrame(json_output['outputs']['horizon_profile']) + df.columns = ['horizon_azimuth', 'horizon_angles'] + return df From d860570eadb938590274ac291239dbd37b22ae5c Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 15:07:39 -0600 Subject: [PATCH 69/94] Update pvgis.py --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 8a6459b9cb..cfbf8a48d7 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -695,7 +695,7 @@ def get_pvgis_horizon(latitude, longitude, `_ ''' params = {'lat': latitude, 'lon': longitude, - 'outputformat': outputformat} + 'outputformat': 'json'} res = requests.get(url + 'printhorizon', params=params, **kwargs) if not res.ok: From 8aedbe66e53799782a30d2fc7a0a5cc03987b00b Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 15:13:21 -0600 Subject: [PATCH 70/94] added url default info --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index cfbf8a48d7..9c30fc07ea 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -679,7 +679,7 @@ def get_pvgis_horizon(latitude, longitude, Latitude in degrees north longitude : float Longitude in degrees east - url: str, default + url: str, default: :const:`pvlib.iotools.pvgis.URL` Base URL for PVGIS kwargs: Passed to requests.get From bef7c8145d61710c1be413ec010654795220c0b3 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Mon, 1 May 2023 15:19:32 -0600 Subject: [PATCH 71/94] remove np dependency --- pvlib/iotools/pvgis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 9c30fc07ea..79d250c99a 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -19,7 +19,6 @@ from pathlib import Path import requests import pandas as pd -import numpy as np from pvlib.iotools import read_epw, parse_epw import warnings from pvlib._deprecation import pvlibDeprecationWarning From d9689d3e5d8c0907d41e9d78906e430f1fa39083 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 11:43:30 -0600 Subject: [PATCH 72/94] Apply suggestions from code review Co-authored-by: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> --- docs/sphinx/source/whatsnew/v0.9.6.rst | 2 +- pvlib/iotools/__init__.py | 2 +- pvlib/iotools/pvgis.py | 16 ++++++---------- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.6.rst b/docs/sphinx/source/whatsnew/v0.9.6.rst index 161792ffda..4aa054a51a 100644 --- a/docs/sphinx/source/whatsnew/v0.9.6.rst +++ b/docs/sphinx/source/whatsnew/v0.9.6.rst @@ -11,7 +11,7 @@ Deprecations Enhancements ~~~~~~~~~~~~ -* Added horizon retrieval from PVGIS function in +* Added function to retrieve horizon data from PVGIS :py:func:`pvlib.iotools.get_pvgis_horizon`. (:issue:`1290`, :pull:`1395`) Bug fixes diff --git a/pvlib/iotools/__init__.py b/pvlib/iotools/__init__.py index 4358a208a5..6f6a254a60 100644 --- a/pvlib/iotools/__init__.py +++ b/pvlib/iotools/__init__.py @@ -14,8 +14,8 @@ from pvlib.iotools.psm3 import parse_psm3 # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_tmy, read_pvgis_tmy # noqa: F401 from pvlib.iotools.pvgis import read_pvgis_hourly # noqa: F401 -from pvlib.iotools.pvgis import get_pvgis_horizon # noqa: F401 from pvlib.iotools.pvgis import get_pvgis_hourly # noqa: F401 +from pvlib.iotools.pvgis import get_pvgis_horizon # noqa: F401 from pvlib.iotools.bsrn import get_bsrn # noqa: F401 from pvlib.iotools.bsrn import read_bsrn # noqa: F401 from pvlib.iotools.bsrn import parse_bsrn # noqa: F401 diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 79d250c99a..82355bfb4c 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -667,10 +667,8 @@ def read_pvgis_tmy(filename, pvgis_format=None, map_variables=None): return data, months_selected, inputs, meta -def get_pvgis_horizon(latitude, longitude, - url=URL, **kwargs): - r''' - Get horizon data from PVGIS +def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): + """Get horizon data from PVGIS. Parameters ---------- @@ -692,11 +690,9 @@ def get_pvgis_horizon(latitude, longitude, ---------- .. [1] `PVGIS horizon profile tool `_ - ''' - params = {'lat': latitude, 'lon': longitude, - 'outputformat': 'json'} - res = requests.get(url + 'printhorizon', params=params, - **kwargs) + """ + params = {'lat': latitude, 'lon': longitude, 'outputformat': 'json'} + res = requests.get(url + 'printhorizon', params=params, **kwargs) if not res.ok: try: err_msg = res.json() @@ -704,7 +700,7 @@ def get_pvgis_horizon(latitude, longitude, res.raise_for_status() else: raise requests.HTTPError(err_msg['message']) - json_output = res.json() + json_output = res.json() df = pd.DataFrame(json_output['outputs']['horizon_profile']) df.columns = ['horizon_azimuth', 'horizon_angles'] return df From 302abb0426c12dae53377658acde884bfb903803 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 11:44:47 -0600 Subject: [PATCH 73/94] Moved ref for consistency --- docs/sphinx/source/reference/iotools.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/sphinx/source/reference/iotools.rst b/docs/sphinx/source/reference/iotools.rst index e6bd33907e..14271cf3ee 100644 --- a/docs/sphinx/source/reference/iotools.rst +++ b/docs/sphinx/source/reference/iotools.rst @@ -31,13 +31,13 @@ of sources and file formats relevant to solar energy modeling. iotools.read_pvgis_tmy iotools.get_pvgis_hourly iotools.read_pvgis_hourly + iotools.get_pvgis_horizon iotools.get_bsrn iotools.read_bsrn iotools.parse_bsrn iotools.get_cams iotools.read_cams iotools.parse_cams - iotools.get_pvgis_horizon A :py:class:`~pvlib.location.Location` object may be created from metadata in some files. From 4ad491bb7f77ce245cec0f1789c6c76515ffa23e Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 11:46:01 -0600 Subject: [PATCH 74/94] Updated name of horizon data variable --- pvlib/tests/iotools/test_pvgis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index e457f3021b..db1c1ea3d5 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -65,7 +65,7 @@ [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] -data_horizon_abq = pd.Series([9.9, 13.0, 14.5, 15.7, 14.9, 15.3, +data_horizon = pd.Series([9.9, 13.0, 14.5, 15.7, 14.9, 15.3, 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, @@ -523,7 +523,7 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): def test_read_pvgis_horizon(): df = get_pvgis_horizon(35.171051, -106.465158) elv = df.horizon_angles - assert_series_equal(elv, data_horizon_abq) + assert_series_equal(elv, data_horizon) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' From 7901797b431267e012f249ac4c8fbc0b02171719 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 14:31:53 -0400 Subject: [PATCH 75/94] Add CSV for horizon tests --- pvlib/data/test_read_pvgis_horizon.csv | 50 ++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 pvlib/data/test_read_pvgis_horizon.csv diff --git a/pvlib/data/test_read_pvgis_horizon.csv b/pvlib/data/test_read_pvgis_horizon.csv new file mode 100644 index 0000000000..44772b37dc --- /dev/null +++ b/pvlib/data/test_read_pvgis_horizon.csv @@ -0,0 +1,50 @@ +,horizon_azimuth,horizon_angles +0,0.0,9.9 +1,7.5,13.0 +2,15.0,14.5 +3,22.5,15.7 +4,30.0,14.9 +5,37.5,15.3 +6,45.0,15.7 +7,52.5,15.7 +8,60.0,13.0 +9,67.5,11.5 +10,75.0,11.1 +11,82.5,11.5 +12,90.0,10.3 +13,97.5,11.5 +14,105.0,10.3 +15,112.5,9.5 +16,120.0,10.7 +17,127.5,11.8 +18,135.0,11.8 +19,142.5,8.8 +20,150.0,8.4 +21,157.5,7.3 +22,165.0,5.7 +23,172.5,5.7 +24,180.0,4.6 +25,187.5,3.4 +26,195.0,0.8 +27,202.5,0.0 +28,210.0,0.0 +29,217.5,0.0 +30,225.0,0.0 +31,232.5,0.0 +32,240.0,0.0 +33,247.5,0.0 +34,255.0,0.0 +35,262.5,0.0 +36,270.0,0.0 +37,277.5,0.0 +38,285.0,0.0 +39,292.5,0.0 +40,300.0,0.0 +41,307.5,0.0 +42,315.0,1.1 +43,322.5,1.9 +44,330.0,3.8 +45,337.5,5.0 +46,345.0,6.5 +47,352.5,9.2 +48,360.0,9.9 From af8988f233af69efef0da61f284466b6d5f1523f Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 12:37:42 -0600 Subject: [PATCH 76/94] Test whole dataframe --- pvlib/tests/iotools/test_pvgis.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index db1c1ea3d5..1603e33e2c 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -65,15 +65,6 @@ [1187.2, 129.59, 8.06, 0.97, 0.97, 0.0], [3950.1, 423.28, 14.8, 1.89, 0.69, 0.0]] -data_horizon = pd.Series([9.9, 13.0, 14.5, 15.7, 14.9, 15.3, - 15.7, 15.7, 13.0, 11.5, 11.1, 11.5, - 10.3, 11.5, 10.3, 9.5, 10.7, 11.8, - 11.8, 8.8, 8.4, 7.3, 5.7, 5.7, 4.6, - 3.4, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 1.1, 1.9, 3.8, 5.0, - 6.5, 9.2, 9.9], name='horizon_angles') - inputs_radiation_csv = {'latitude': 45.0, 'longitude': 8.0, 'elevation': 250.0, 'radiation_database': 'PVGIS-SARAH', 'Slope': '30 deg.', 'Azimuth': '0 deg.'} @@ -521,9 +512,9 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data def test_read_pvgis_horizon(): - df = get_pvgis_horizon(35.171051, -106.465158) - elv = df.horizon_angles - assert_series_equal(elv, data_horizon) + pvgis_data = get_pvgis_horizon(35.171051, -106.465158) + horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv') + assert_frame_equal(pvgis_data, horizon_data) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' From f58f8d7eb291f74b6e613d8db0bc0a4c7dd5cac1 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 14:10:16 -0600 Subject: [PATCH 77/94] Change coordinate system to pvlib --- pvlib/iotools/pvgis.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 82355bfb4c..678d018e00 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -703,4 +703,5 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): json_output = res.json() df = pd.DataFrame(json_output['outputs']['horizon_profile']) df.columns = ['horizon_azimuth', 'horizon_angles'] + df['horizon_azimuth'] += 180 return df From a305ab6736ecc58814c83c3f5f5367dfc614431c Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Wed, 3 May 2023 14:12:57 -0600 Subject: [PATCH 78/94] fix test to read in csv properly --- pvlib/tests/iotools/test_pvgis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 1603e33e2c..1545642231 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -513,7 +513,8 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data def test_read_pvgis_horizon(): pvgis_data = get_pvgis_horizon(35.171051, -106.465158) - horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv') + horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv', + index_col = 0) assert_frame_equal(pvgis_data, horizon_data) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): From 15e5c9ceee2a6521f78db5d88f0105af14c97d5d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 4 May 2023 09:46:08 -0600 Subject: [PATCH 79/94] Apply suggestions from code review Co-authored-by: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com> --- pvlib/iotools/pvgis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 678d018e00..a0fdb9fbc6 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -683,7 +683,7 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): Returns ------- - df : pd.DataFrame + data : pd.DataFrame Pandas dataframe of the retrived horizon References @@ -703,5 +703,6 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): json_output = res.json() df = pd.DataFrame(json_output['outputs']['horizon_profile']) df.columns = ['horizon_azimuth', 'horizon_angles'] + # Convert azimuth to pvlib convention (north=0, south=180) df['horizon_azimuth'] += 180 return df From 482bc23944ca43dd93e9f56096ed5325e002d88b Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 4 May 2023 09:48:32 -0600 Subject: [PATCH 80/94] Changed names for consistency --- pvlib/iotools/pvgis.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index a0fdb9fbc6..153be2b7b6 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -701,8 +701,8 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): else: raise requests.HTTPError(err_msg['message']) json_output = res.json() - df = pd.DataFrame(json_output['outputs']['horizon_profile']) - df.columns = ['horizon_azimuth', 'horizon_angles'] + data = pd.DataFrame(json_output['outputs']['horizon_profile']) + data.columns = ['horizon_azimuth', 'horizon_elevation'] # Convert azimuth to pvlib convention (north=0, south=180) - df['horizon_azimuth'] += 180 - return df + data['horizon_azimuth'] += 180 + return data From 94fee4c2832aec0db6e2a198dda48495861b79d1 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Thu, 4 May 2023 09:49:22 -0600 Subject: [PATCH 81/94] Match column names in test dataframe --- pvlib/data/test_read_pvgis_horizon.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/data/test_read_pvgis_horizon.csv b/pvlib/data/test_read_pvgis_horizon.csv index 44772b37dc..d54abcb242 100644 --- a/pvlib/data/test_read_pvgis_horizon.csv +++ b/pvlib/data/test_read_pvgis_horizon.csv @@ -1,4 +1,4 @@ -,horizon_azimuth,horizon_angles +,horizon_azimuth,horizon_elevation 0,0.0,9.9 1,7.5,13.0 2,15.0,14.5 From 2eea3845592187d319587fdf4205c9a316102683 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 9 May 2023 10:54:52 -0600 Subject: [PATCH 82/94] stickler --- pvlib/tests/iotools/test_pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 1545642231..66c8122f36 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -514,7 +514,7 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): def test_read_pvgis_horizon(): pvgis_data = get_pvgis_horizon(35.171051, -106.465158) horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv', - index_col = 0) + index_col=0) assert_frame_equal(pvgis_data, horizon_data) def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): From bf342561b5b88f59c3b9e9c6e436ee8a36ef3d8d Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 9 May 2023 11:03:42 -0600 Subject: [PATCH 83/94] Set horizon_azimuth to the index --- pvlib/iotools/pvgis.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 153be2b7b6..5dd6d713f1 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -705,4 +705,5 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): data.columns = ['horizon_azimuth', 'horizon_elevation'] # Convert azimuth to pvlib convention (north=0, south=180) data['horizon_azimuth'] += 180 + data.set_index('horizon_azimuth', inplace=True) return data From 01893d871bbb07e26096924d6fcf2d8420539572 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 9 May 2023 11:36:45 -0600 Subject: [PATCH 84/94] Added invalid coords test --- pvlib/tests/iotools/test_pvgis.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 66c8122f36..ee7d0908a2 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -517,6 +517,13 @@ def test_read_pvgis_horizon(): index_col=0) assert_frame_equal(pvgis_data, horizon_data) + +@pytest.mark.remote_data +def test_read_pvgis_horizon_invalid_coords(): + with pytest.warns(UserWarning, match='lat: Incorrect value'): + _ = get_pvgis_horizon(100, 50) # unfeasible latitude + + def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): fn = DATA_DIR / 'tmy_45.000_8.000_2005_2016.json' actual, _, _, _ = read_pvgis_tmy(fn, map_variables=True) From d978b2c490cee4ea01ffc3b17655cc08c6105ba7 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 9 May 2023 12:13:19 -0600 Subject: [PATCH 85/94] update test to reflect change to index --- pvlib/tests/iotools/test_pvgis.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index ee7d0908a2..1194c2a12f 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -515,6 +515,7 @@ def test_read_pvgis_horizon(): pvgis_data = get_pvgis_horizon(35.171051, -106.465158) horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv', index_col=0) + horizon_data.set_index('horizon_azimuth', inplace=True) assert_frame_equal(pvgis_data, horizon_data) From b69381312e83b1ba716c1cc451622ab5eed40f07 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 9 May 2023 14:09:58 -0600 Subject: [PATCH 86/94] forgot xfail flag --- pvlib/tests/iotools/test_pvgis.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 1194c2a12f..792a17ccbc 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -520,6 +520,7 @@ def test_read_pvgis_horizon(): @pytest.mark.remote_data +@pytest.mark.xfail(strict=True, reason='Ensure PVGIS raises exception') def test_read_pvgis_horizon_invalid_coords(): with pytest.warns(UserWarning, match='lat: Incorrect value'): _ = get_pvgis_horizon(100, 50) # unfeasible latitude From 8de25069c1cdb5e9e321c6ee94a93c318a73b598 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Thu, 11 May 2023 18:20:28 +0200 Subject: [PATCH 87/94] Update test_read_pvgis_horzion_invalid_coords --- pvlib/tests/iotools/test_pvgis.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 792a17ccbc..6c6c39b67f 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -520,9 +520,8 @@ def test_read_pvgis_horizon(): @pytest.mark.remote_data -@pytest.mark.xfail(strict=True, reason='Ensure PVGIS raises exception') def test_read_pvgis_horizon_invalid_coords(): - with pytest.warns(UserWarning, match='lat: Incorrect value'): + with pytest.raises(requests.HTTPError, match='lat: Incorrect value'): _ = get_pvgis_horizon(100, 50) # unfeasible latitude From 568c9cdd6a53b432ebb1e0b9f89b5a7397a48863 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 23 May 2023 15:15:58 -0600 Subject: [PATCH 88/94] return metadata object --- pvlib/iotools/pvgis.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 5dd6d713f1..91784f87e1 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -685,6 +685,8 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): ------- data : pd.DataFrame Pandas dataframe of the retrived horizon + metadata : dict + Metadata returned by PVGIS References ---------- @@ -701,9 +703,10 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): else: raise requests.HTTPError(err_msg['message']) json_output = res.json() + metadata = json_output['meta'] data = pd.DataFrame(json_output['outputs']['horizon_profile']) data.columns = ['horizon_azimuth', 'horizon_elevation'] # Convert azimuth to pvlib convention (north=0, south=180) data['horizon_azimuth'] += 180 data.set_index('horizon_azimuth', inplace=True) - return data + return data, meta From fb3d5712015e474d01a3eba829ef69828990bc63 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 23 May 2023 15:21:28 -0600 Subject: [PATCH 89/94] Update pvgis.py --- pvlib/iotools/pvgis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 91784f87e1..e6920a0baa 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -709,4 +709,4 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): # Convert azimuth to pvlib convention (north=0, south=180) data['horizon_azimuth'] += 180 data.set_index('horizon_azimuth', inplace=True) - return data, meta + return data, metadata From d3bc17b2f33b56b50cb2c18edbd4a872ca638931 Mon Sep 17 00:00:00 2001 From: Ben Pierce Date: Tue, 23 May 2023 15:23:40 -0600 Subject: [PATCH 90/94] fix tests for metadata update --- pvlib/tests/iotools/test_pvgis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 6c6c39b67f..2c9a9deeb1 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -512,7 +512,7 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data def test_read_pvgis_horizon(): - pvgis_data = get_pvgis_horizon(35.171051, -106.465158) + pvgis_data, _ = get_pvgis_horizon(35.171051, -106.465158) horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv', index_col=0) horizon_data.set_index('horizon_azimuth', inplace=True) @@ -522,7 +522,7 @@ def test_read_pvgis_horizon(): @pytest.mark.remote_data def test_read_pvgis_horizon_invalid_coords(): with pytest.raises(requests.HTTPError, match='lat: Incorrect value'): - _ = get_pvgis_horizon(100, 50) # unfeasible latitude + _, _ = get_pvgis_horizon(100, 50) # unfeasible latitude def test_read_pvgis_tmy_map_variables(pvgis_tmy_mapped_columns): From a696fafea04eab051ff9e9a2f448087ade3e509b Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Thu, 25 May 2023 16:04:16 +0200 Subject: [PATCH 91/94] Remove index column from test file --- pvlib/data/test_read_pvgis_horizon.csv | 100 ++++++++++++------------- 1 file changed, 50 insertions(+), 50 deletions(-) diff --git a/pvlib/data/test_read_pvgis_horizon.csv b/pvlib/data/test_read_pvgis_horizon.csv index d54abcb242..d85dbb6e6f 100644 --- a/pvlib/data/test_read_pvgis_horizon.csv +++ b/pvlib/data/test_read_pvgis_horizon.csv @@ -1,50 +1,50 @@ -,horizon_azimuth,horizon_elevation -0,0.0,9.9 -1,7.5,13.0 -2,15.0,14.5 -3,22.5,15.7 -4,30.0,14.9 -5,37.5,15.3 -6,45.0,15.7 -7,52.5,15.7 -8,60.0,13.0 -9,67.5,11.5 -10,75.0,11.1 -11,82.5,11.5 -12,90.0,10.3 -13,97.5,11.5 -14,105.0,10.3 -15,112.5,9.5 -16,120.0,10.7 -17,127.5,11.8 -18,135.0,11.8 -19,142.5,8.8 -20,150.0,8.4 -21,157.5,7.3 -22,165.0,5.7 -23,172.5,5.7 -24,180.0,4.6 -25,187.5,3.4 -26,195.0,0.8 -27,202.5,0.0 -28,210.0,0.0 -29,217.5,0.0 -30,225.0,0.0 -31,232.5,0.0 -32,240.0,0.0 -33,247.5,0.0 -34,255.0,0.0 -35,262.5,0.0 -36,270.0,0.0 -37,277.5,0.0 -38,285.0,0.0 -39,292.5,0.0 -40,300.0,0.0 -41,307.5,0.0 -42,315.0,1.1 -43,322.5,1.9 -44,330.0,3.8 -45,337.5,5.0 -46,345.0,6.5 -47,352.5,9.2 -48,360.0,9.9 +horizon_azimuth,horizon_elevation +0,9.9 +7.5,13 +15,14.5 +22.5,15.7 +30,14.9 +37.5,15.3 +45,15.7 +52.5,15.7 +60,13 +67.5,11.5 +75,11.1 +82.5,11.5 +90,10.3 +97.5,11.5 +105,10.3 +112.5,9.5 +120,10.7 +127.5,11.8 +135,11.8 +142.5,8.8 +150,8.4 +157.5,7.3 +165,5.7 +172.5,5.7 +180,4.6 +187.5,3.4 +195,0.8 +202.5,0 +210,0 +217.5,0 +225,0 +232.5,0 +240,0 +247.5,0 +255,0 +262.5,0 +270,0 +277.5,0 +285,0 +292.5,0 +300,0 +307.5,0 +315,1.1 +322.5,1.9 +330,3.8 +337.5,5 +345,6.5 +352.5,9.2 +360,9.9 From b65d58b853a35e4569796456b7ab46b6f8009518 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Thu, 25 May 2023 16:04:38 +0200 Subject: [PATCH 92/94] Convert output to pd.Series --- pvlib/iotools/pvgis.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index e6920a0baa..3a19262c8c 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -683,10 +683,17 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): Returns ------- - data : pd.DataFrame - Pandas dataframe of the retrived horizon + data : pd.Series + Pandas Series of the retrived horizon elevation angles. Index is the + corresponding horizon azimuth angles. metadata : dict - Metadata returned by PVGIS + Metadata returned by PVGIS. + + Notes + ----- + The horizon azimuths are specified clockwise from north, e.g., south=180. + This is the standard pvlib convention, although the PVGIS website specifies + south=0. References ---------- @@ -709,4 +716,5 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): # Convert azimuth to pvlib convention (north=0, south=180) data['horizon_azimuth'] += 180 data.set_index('horizon_azimuth', inplace=True) + data = data['horizon_elevation'] # convert to pd.Series return data, metadata From d008dab80832d442b1f7aa992bc174f44532de40 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Thu, 25 May 2023 16:04:52 +0200 Subject: [PATCH 93/94] Update tests to assert pd.Series --- pvlib/tests/iotools/test_pvgis.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py index 2c9a9deeb1..a5a5e3fbd7 100644 --- a/pvlib/tests/iotools/test_pvgis.py +++ b/pvlib/tests/iotools/test_pvgis.py @@ -511,15 +511,17 @@ def test_get_pvgis_map_variables(pvgis_tmy_mapped_columns): @pytest.mark.remote_data +@pytest.mark.flaky(reruns=RERUNS, reruns_delay=RERUNS_DELAY) def test_read_pvgis_horizon(): pvgis_data, _ = get_pvgis_horizon(35.171051, -106.465158) horizon_data = pd.read_csv(DATA_DIR / 'test_read_pvgis_horizon.csv', index_col=0) - horizon_data.set_index('horizon_azimuth', inplace=True) - assert_frame_equal(pvgis_data, horizon_data) + horizon_data = horizon_data['horizon_elevation'] + assert_series_equal(pvgis_data, horizon_data) @pytest.mark.remote_data +@pytest.mark.flaky(reruns=RERUNS, reruns_delay=RERUNS_DELAY) def test_read_pvgis_horizon_invalid_coords(): with pytest.raises(requests.HTTPError, match='lat: Incorrect value'): _, _ = get_pvgis_horizon(100, 50) # unfeasible latitude From 22e09883039d735543b69da3b405f166fb53ae3f Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Thu, 25 May 2023 16:12:02 +0200 Subject: [PATCH 94/94] Remove duplicate north value (0, 360) --- pvlib/data/test_read_pvgis_horizon.csv | 1 - pvlib/iotools/pvgis.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/pvlib/data/test_read_pvgis_horizon.csv b/pvlib/data/test_read_pvgis_horizon.csv index d85dbb6e6f..0b3ed27bbc 100644 --- a/pvlib/data/test_read_pvgis_horizon.csv +++ b/pvlib/data/test_read_pvgis_horizon.csv @@ -47,4 +47,3 @@ horizon_azimuth,horizon_elevation 337.5,5 345,6.5 352.5,9.2 -360,9.9 diff --git a/pvlib/iotools/pvgis.py b/pvlib/iotools/pvgis.py index 3a19262c8c..16bfee7cee 100644 --- a/pvlib/iotools/pvgis.py +++ b/pvlib/iotools/pvgis.py @@ -717,4 +717,5 @@ def get_pvgis_horizon(latitude, longitude, url=URL, **kwargs): data['horizon_azimuth'] += 180 data.set_index('horizon_azimuth', inplace=True) data = data['horizon_elevation'] # convert to pd.Series + data = data[data.index < 360] # remove duplicate north point (0 and 360) return data, metadata