diff --git a/docs/sphinx/source/reference/index.rst b/docs/sphinx/source/reference/index.rst index 0f1902325a..b96ca4114f 100644 --- a/docs/sphinx/source/reference/index.rst +++ b/docs/sphinx/source/reference/index.rst @@ -20,3 +20,4 @@ API reference modelchain bifacial scaling + location diff --git a/docs/sphinx/source/reference/location.rst b/docs/sphinx/source/reference/location.rst new file mode 100644 index 0000000000..28ef46812a --- /dev/null +++ b/docs/sphinx/source/reference/location.rst @@ -0,0 +1,11 @@ +.. currentmodule:: pvlib + +Location +======== + +Methods for information about locations. + +.. autosummary:: + :toctree: generated/ + + location.lookup_altitude diff --git a/docs/sphinx/source/whatsnew/v0.9.3.rst b/docs/sphinx/source/whatsnew/v0.9.3.rst index d916ffb558..2f3c643d05 100644 --- a/docs/sphinx/source/whatsnew/v0.9.3.rst +++ b/docs/sphinx/source/whatsnew/v0.9.3.rst @@ -9,7 +9,9 @@ Deprecations Enhancements ~~~~~~~~~~~~ - +* Low resolution altitude lookup map + :py:func:`~pvlib.location.lookup_altitude` + (:issue:`1516`, :pull:`1518`) Bug fixes ~~~~~~~~~ @@ -33,4 +35,5 @@ Requirements Contributors ~~~~~~~~~~~~ -João Guilherme (:ghuser:`joaoguilhermeS`) +* João Guilherme (:ghuser:`joaoguilhermeS`) +* Nicolas Martinez (:ghuser:`nicomt`) diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py index 2f5cd68c62..ffc7ac9b55 100644 --- a/pvlib/clearsky.py +++ b/pvlib/clearsky.py @@ -14,6 +14,7 @@ import h5py from pvlib import atmosphere, tools +from pvlib.tools import _degrees_to_index def ineichen(apparent_zenith, airmass_absolute, linke_turbidity, @@ -286,67 +287,6 @@ def _calendar_month_middles(year): return middles -def _degrees_to_index(degrees, coordinate): - """Transform input degrees to an output index integer. The Linke - turbidity lookup tables have three dimensions, latitude, longitude, and - month. Specify a degree value and either 'latitude' or 'longitude' to get - the appropriate index number for the first two of these index numbers. - - Parameters - ---------- - degrees : float or int - Degrees of either latitude or longitude. - coordinate : string - Specify whether degrees arg is latitude or longitude. Must be set to - either 'latitude' or 'longitude' or an error will be raised. - - Returns - ------- - index : np.int16 - The latitude or longitude index number to use when looking up values - in the Linke turbidity lookup table. - """ - # Assign inputmin, inputmax, and outputmax based on degree type. - if coordinate == 'latitude': - inputmin = 90 - inputmax = -90 - outputmax = 2160 - elif coordinate == 'longitude': - inputmin = -180 - inputmax = 180 - outputmax = 4320 - else: - raise IndexError("coordinate must be 'latitude' or 'longitude'.") - - inputrange = inputmax - inputmin - scale = outputmax/inputrange # number of indices per degree - center = inputmin + 1 / scale / 2 # shift to center of index - outputmax -= 1 # shift index to zero indexing - index = (degrees - center) * scale - err = IndexError('Input, %g, is out of range (%g, %g).' % - (degrees, inputmin, inputmax)) - - # If the index is still out of bounds after rounding, raise an error. - # 0.500001 is used in comparisons instead of 0.5 to allow for a small - # margin of error which can occur when dealing with floating point numbers. - if index > outputmax: - if index - outputmax <= 0.500001: - index = outputmax - else: - raise err - elif index < 0: - if -index <= 0.500001: - index = 0 - else: - raise err - # If the index wasn't set to outputmax or 0, round it and cast it as an - # integer so it can be used in integer-based indexing. - else: - index = int(np.around(index)) - - return index - - def haurwitz(apparent_zenith): ''' Determine clear sky GHI using the Haurwitz model. diff --git a/pvlib/data/Altitude.h5 b/pvlib/data/Altitude.h5 new file mode 100644 index 0000000000..2921233212 Binary files /dev/null and b/pvlib/data/Altitude.h5 differ diff --git a/pvlib/location.py b/pvlib/location.py index ddd2d1b96a..0fb3b42a04 100644 --- a/pvlib/location.py +++ b/pvlib/location.py @@ -4,13 +4,16 @@ # Will Holmgren, University of Arizona, 2014-2016. +import os import datetime import warnings import pandas as pd import pytz +import h5py from pvlib import solarposition, clearsky, atmosphere, irradiance +from pvlib.tools import _degrees_to_index class Location: """ @@ -356,3 +359,88 @@ def get_sun_rise_set_transit(self, times, method='pyephem', **kwargs): 'one of pyephem, spa, geometric' .format(method)) return result + + +def lookup_altitude(latitude, longitude): + """ + Look up location altitude from low-resolution altitude map + supplied with pvlib. The data for this map comes from multiple open data + sources with varying resolutions aggregated by Mapzen. + + More details can be found here + https://github.com/tilezen/joerd/blob/master/docs/data-sources.md + + Altitudes from this map are a coarse approximation and can have + significant errors (100+ meters) introduced by downsampling and + source data resolution. + + Parameters + ---------- + latitude : float. + Positive is north of the equator. + Use decimal degrees notation. + + longitude : float. + Positive is east of the prime meridian. + Use decimal degrees notation. + + Returns + ------- + altitude : float + The altitude of the location in meters. + + Notes + ----------- + Attributions: + + * ArcticDEM terrain data DEM(s) were created from DigitalGlobe, Inc., + imagery and funded under National Science Foundation awards 1043681, + 1559691, and 1542736; + * Australia terrain data © Commonwealth of Australia + (Geoscience Australia) 2017; + * Austria terrain data © offene Daten Österreichs - Digitales + Geländemodell (DGM) Österreich; + * Canada terrain data contains information licensed under the Open + Government Licence - Canada; + * Europe terrain data produced using Copernicus data and information + funded by the European Union - EU-DEM layers; + * Global ETOPO1 terrain data U.S. National Oceanic and Atmospheric + Administration + * Mexico terrain data source: INEGI, Continental relief, 2016; + * New Zealand terrain data Copyright 2011 Crown copyright (c) Land + Information New Zealand and the New Zealand Government + (All rights reserved); + * Norway terrain data © Kartverket; + * United Kingdom terrain data © Environment Agency copyright and/or + database right 2015. All rights reserved; + * United States 3DEP (formerly NED) and global GMTED2010 and SRTM + terrain data courtesy of the U.S. Geological Survey. + + References + ---------- + .. [1] `Mapzen, Linux foundation project for open data maps + `_ + .. [2] `Joerd, tool for downloading and processing DEMs, Used by Mapzen + `_ + .. [3] `AWS, Open Data Registry Terrain Tiles + `_ + + """ + + pvlib_path = os.path.dirname(os.path.abspath(__file__)) + filepath = os.path.join(pvlib_path, 'data', 'Altitude.h5') + + latitude_index = _degrees_to_index(latitude, coordinate='latitude') + longitude_index = _degrees_to_index(longitude, coordinate='longitude') + + with h5py.File(filepath, 'r') as alt_h5_file: + alt = alt_h5_file['Altitude'][latitude_index, longitude_index] + + # 255 is a special value that means nodata. Fallback to 0 if nodata. + if alt == 255: + return 0 + # Altitude is encoded in 28 meter steps from -450 meters to 6561 meters + # There are 0-254 possible altitudes, with 255 reserved for nodata. + alt *= 28 + alt -= 450 + return float(alt) diff --git a/pvlib/tests/test_clearsky.py b/pvlib/tests/test_clearsky.py index d603cbcdfe..6e620f3c79 100644 --- a/pvlib/tests/test_clearsky.py +++ b/pvlib/tests/test_clearsky.py @@ -511,13 +511,6 @@ def monthly_lt_nointerp(lat, lon, time=months): monthly_lt_nointerp(38.2, -181) # exceeds min longitude -def test_degrees_to_index_1(): - """Test that _degrees_to_index raises an error when something other than - 'latitude' or 'longitude' is passed.""" - with pytest.raises(IndexError): # invalid value for coordinate argument - clearsky._degrees_to_index(degrees=22.0, coordinate='width') - - @pytest.fixture def detect_clearsky_data(): data_file = DATA_DIR / 'detect_clearsky_data.csv' diff --git a/pvlib/tests/test_location.py b/pvlib/tests/test_location.py index a4b69e5cb2..b3c1576bdf 100644 --- a/pvlib/tests/test_location.py +++ b/pvlib/tests/test_location.py @@ -12,7 +12,7 @@ from pytz.exceptions import UnknownTimeZoneError import pvlib -from pvlib.location import Location +from pvlib.location import Location, lookup_altitude from pvlib.solarposition import declination_spencer71 from pvlib.solarposition import equation_of_time_spencer71 from .conftest import requires_ephem @@ -326,3 +326,23 @@ def test_get_sun_rise_set_transit_valueerror(golden): def test_extra_kwargs(): with pytest.raises(TypeError, match='arbitrary_kwarg'): Location(32.2, -111, arbitrary_kwarg='value') + + +def test_lookup_altitude(): + max_alt_error = 125 + # location name, latitude, longitude, altitude + test_locations = [ + ('Tucson, USA', 32.2540, -110.9742, 724), + ('Lusaka, Zambia', -15.3875, 28.3228, 1253), + ('Tokio, Japan', 35.6762, 139.6503, 40), + ('Canberra, Australia', -35.2802, 149.1310, 566), + ('Bogota, Colombia', 4.7110, -74.0721, 2555), + ('Dead Sea, West Bank', 31.525849, 35.449214, -415), + ('New Delhi, India', 28.6139, 77.2090, 214), + ('Null Island, Atlantic Ocean', 0, 0, 0), + ] + + for name, lat, lon, expected_alt in test_locations: + alt_found = lookup_altitude(lat, lon) + assert abs(alt_found - expected_alt) < max_alt_error, \ + f'Max error exceded for {name} - e: {expected_alt} f: {alt_found}' diff --git a/pvlib/tests/test_tools.py b/pvlib/tests/test_tools.py index e6b8e9091a..167dca8cec 100644 --- a/pvlib/tests/test_tools.py +++ b/pvlib/tests/test_tools.py @@ -72,3 +72,10 @@ def test__golden_sect_DataFrame_nans(): v, x = tools._golden_sect_DataFrame(params, lower, upper, _obj_test_golden_sect) assert np.allclose(x, expected, atol=1e-8, equal_nan=True) + + +def test_degrees_to_index_1(): + """Test that _degrees_to_index raises an error when something other than + 'latitude' or 'longitude' is passed.""" + with pytest.raises(IndexError): # invalid value for coordinate argument + tools._degrees_to_index(degrees=22.0, coordinate='width') diff --git a/pvlib/tools.py b/pvlib/tools.py index dd5730ee6f..991568f9e0 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -412,3 +412,61 @@ def _get_sample_intervals(times, win_length): 'periods, leap days, etc.' ) raise NotImplementedError(message) + + +def _degrees_to_index(degrees, coordinate): + """Transform input degrees to an output index integer. + Specify a degree value and either 'latitude' or 'longitude' to get + the appropriate index number for these two index numbers. + Parameters + ---------- + degrees : float or int + Degrees of either latitude or longitude. + coordinate : string + Specify whether degrees arg is latitude or longitude. Must be set to + either 'latitude' or 'longitude' or an error will be raised. + Returns + ------- + index : np.int16 + The latitude or longitude index number to use when looking up values + in the Linke turbidity lookup table. + """ + # Assign inputmin, inputmax, and outputmax based on degree type. + if coordinate == 'latitude': + inputmin = 90 + inputmax = -90 + outputmax = 2160 + elif coordinate == 'longitude': + inputmin = -180 + inputmax = 180 + outputmax = 4320 + else: + raise IndexError("coordinate must be 'latitude' or 'longitude'.") + + inputrange = inputmax - inputmin + scale = outputmax/inputrange # number of indices per degree + center = inputmin + 1 / scale / 2 # shift to center of index + outputmax -= 1 # shift index to zero indexing + index = (degrees - center) * scale + err = IndexError('Input, %g, is out of range (%g, %g).' % + (degrees, inputmin, inputmax)) + + # If the index is still out of bounds after rounding, raise an error. + # 0.500001 is used in comparisons instead of 0.5 to allow for a small + # margin of error which can occur when dealing with floating point numbers. + if index > outputmax: + if index - outputmax <= 0.500001: + index = outputmax + else: + raise err + elif index < 0: + if -index <= 0.500001: + index = 0 + else: + raise err + # If the index wasn't set to outputmax or 0, round it and cast it as an + # integer so it can be used in integer-based indexing. + else: + index = int(np.around(index)) + + return index