diff --git a/hopp/simulation/technologies/resource/__init__.py b/hopp/simulation/technologies/resource/__init__.py index 7a8ea8f0b..0d3601457 100644 --- a/hopp/simulation/technologies/resource/__init__.py +++ b/hopp/simulation/technologies/resource/__init__.py @@ -5,3 +5,5 @@ from hopp.simulation.technologies.resource.resource import Resource from hopp.simulation.technologies.resource.greet_data import GREETData from hopp.simulation.technologies.resource.cambium_data import CambiumData +from hopp.simulation.technologies.resource.nsrdb_data import HPCSolarData +from hopp.simulation.technologies.resource.wind_toolkit_data import HPCWindData diff --git a/hopp/simulation/technologies/resource/nsrdb_data.py b/hopp/simulation/technologies/resource/nsrdb_data.py new file mode 100644 index 000000000..5ab132cb9 --- /dev/null +++ b/hopp/simulation/technologies/resource/nsrdb_data.py @@ -0,0 +1,178 @@ +from rex import NSRDBX +from rex.sam_resource import SAMResource +import numpy as np +from hopp.simulation.technologies.resource.resource import Resource +from typing import Optional, Union +from pathlib import Path +import os +# import pandas as pd +NSRDB_DEP = "/kfs2/datasets/NSRDB/deprecated_v3/nsrdb_" +NSRDB_NEW = "/kfs2/datasets/NSRDB/current/nsrdb_" +# Pull Solar Resource Data directly from NSRDB on HPC +# To be called instead of SolarResource from hopp.simulation.technologies.resource +class HPCSolarData(Resource): + + def __init__( + self, + lat: float, + lon: float, + year: int, + nsrdb_source_path: Union[str,Path] = "", + filepath: str = "", + **kwargs): + """ + Input: + lat (float): site latitude + lon (float): site longitude + resource_year (int): year to get resource data for + filepath (str): filepath for nsrdb data on HPC + + Output: self.data + dictionary: + tz: float + elev: float + lat: float + lon: float + year: list of floats + month: list of floats + day: list of floats + hour: list of floats + minute: list of floats + dn: list of floats + df: list of floats + gh: list of floats + wspd: list of floats + tdry: list of floats + pres: list of floats + tdew: list of floats + """ + # NOTE: self.data must be compatible with PVWatts.SolarResource.solar_resource_data to https://nrel-pysam.readthedocs.io/en/main/modules/Pvwattsv8.html#PySAM.Pvwattsv8.Pvwattsv8.SolarResource + self.latitude = lat + self.longitude = lon + self.year = year + super().__init__(lat, lon, year) + + if filepath == "" and nsrdb_source_path=="": + self.nsrdb_file = NSRDB_NEW + "{}.h5".format(self.year) + elif filepath != "" and nsrdb_source_path == "": + self.nsrdb_file = filepath + elif filepath=="" and nsrdb_source_path !="": + self.nsrdb_file = os.path.join(str(nsrdb_source_path),"nsrdb_{}.h5".format(self.year)) + else: + self.nsrdb_file = NSRDB_NEW + "{}.h5".format(self.year) + # Pull data from HPC NSRDB dataset + # self.extract_resource() + self.download_resource() + + # Set solar resource data into SAM/PySAM digestible format + self.format_data() + + # Define final dictionary + + + # def extract_resource(self): + def download_resource(self): + # Define file to download from + # NOTE: HOPP is currently calling an old deprecated version of PSM v3.1 which corresponds to /api/nsrdb/v2/solar/psm3-download + + + # NOTE: Current version of PSM v3.2.2 which corresponds to /api/nsrdb/v2/solar/psm3-2-2-download + + + # Open file with rex NSRDBX object + with NSRDBX(self.nsrdb_file, hsds=False) as f: + # get gid of location closest to given lat/lon coordinates + site_gid = f.lat_lon_gid((self.latitude,self.longitude)) + + # extract timezone, elevation, latitude and longitude from meta dataset with gid + self.time_zone = f.meta['timezone'].iloc[site_gid] + self.elevation = f.meta['elevation'].iloc[site_gid] + self.nsrdb_latitude = f.meta['latitude'].iloc[site_gid] + self.nsrdb_longitude = f.meta['longitude'].iloc[site_gid] + + # extract remaining datapoints: year, month, day, hour, minute, dn, df, gh, wspd,tdry, pres, tdew + # NOTE: datasets have readings at 0 and 30 minutes each hour, HOPP/SAM workflow requires only 30 minute reading values -> filter 0 minute readings with [1::2] + # NOTE: datasets are not auto shifted by timezone offset -> wrap extraction in SAMResource.roll_timeseries(input_array, timezone, #steps in an hour=1) to roll timezones + # NOTE: solar_resource.py code references solar_zenith_angle and RH = relative_humidity but I couldn't find them actually being utilized. Captured them below just in case. + self.year_arr = f.time_index.year.values[1::2] + self.month_arr = f.time_index.month.values[1::2] + self.day_arr = f.time_index.day.values[1::2] + self.hour_arr = f.time_index.hour.values[1::2] + self.minute_arr = f.time_index.minute.values[1::2] + self.dni_arr = SAMResource.roll_timeseries((f['dni', :, site_gid][1::2]), self.time_zone, 1) + self.dhi_arr = SAMResource.roll_timeseries((f['dhi', :, site_gid][1::2]), self.time_zone, 1) + self.ghi_arr = SAMResource.roll_timeseries((f['ghi', :, site_gid][1::2]), self.time_zone, 1) + self.wspd_arr = SAMResource.roll_timeseries((f['wind_speed', :, site_gid][1::2]), self.time_zone, 1) + self.tdry_arr = SAMResource.roll_timeseries((f['air_temperature', :, site_gid][1::2]), self.time_zone, 1) + # self.relative_humidity_arr = SAMResource.roll_timeseries((f['relative_humidity', :, site_gid][1::2]), self.time_zone, 1) + # self.solar_zenith_arr = SAMResource.roll_timeseries((f['solar_zenith_angle', :, site_gid][1::2]), self.time_zone, 1) + self.pres_arr = SAMResource.roll_timeseries((f['surface_pressure', :, site_gid][1::2]), self.time_zone, 1) + self.tdew_arr = SAMResource.roll_timeseries((f['dew_point', :, site_gid][1::2]), self.time_zone, 1) + + self.site_gid = site_gid + + + def format_data(self): + # Remove data from feb29 on leap years + if (self.year % 4) == 0: + feb29 = np.arange(1416,1440) + self.year_arr = np.delete(self.year_arr, feb29) + self.month_arr = np.delete(self.month_arr, feb29) + self.day_arr = np.delete(self.day_arr, feb29) + self.hour_arr = np.delete(self.hour_arr, feb29) + self.minute_arr = np.delete(self.minute_arr, feb29) + self.dni_arr = np.delete(self.dni_arr, feb29) + self.dhi_arr = np.delete(self.dhi_arr, feb29) + self.ghi_arr = np.delete(self.ghi_arr, feb29) + self.wspd_arr = np.delete(self.wspd_arr, feb29) + self.tdry_arr = np.delete(self.tdry_arr, feb29) + # self.relative_humidity_arr = np.delete(self.relative_humidity_arr, feb29) + # self.solar_zenith_arr = np.delete(self.solar_zenith_arr, feb29) + self.pres_arr = np.delete(self.pres_arr, feb29) + self.tdew_arr = np.delete(self.tdew_arr, feb29) + + # round to desired precision and convert to desired data type + # NOTE: unsure if SAM/PySAM is sensitive to data types and decimal precision. If not sensitive, can remove .astype() and round() to increase computational efficiency + self.time_zone = float(self.time_zone) + self.elevation = round(float(self.elevation), 0) + self.nsrdb_latitude = round(float(self.nsrdb_latitude), 2) + self.nsrdb_longitude = round(float(self.nsrdb_longitude),2) + self.year_arr = list(self.year_arr.astype(float, copy=False)) + self.month_arr = list(self.month_arr.astype(float, copy=False)) + self.day_arr = list(self.day_arr.astype(float, copy=False)) + self.hour_arr = list(self.hour_arr.astype(float, copy=False)) + self.minute_arr = list(self.minute_arr.astype(float, copy=False)) + self.dni_arr = list(self.dni_arr.astype(float, copy=False)) + self.dhi_arr = list(self.dhi_arr.astype(float, copy=False)) + self.ghi_arr = list(self.ghi_arr.astype(float, copy=False)) + self.wspd_arr = list(self.wspd_arr.astype(float, copy=False)) + self.tdry_arr = list(self.tdry_arr.astype(float, copy=False)) + # self.relative_humidity_arr = list(np.round(self.relative_humidity_arr, decimals=1)) + # self.solar_zenith_angle_arr = list(np.round(self.solar_zenith_angle_arr, decimals=1)) + self.pres_arr = list(self.pres_arr.astype(float, copy=False)) + self.tdew_arr = list(self.tdew_arr.astype(float, copy=False)) + self.data = {} #unsure if this will cause problem + @Resource.data.setter + def data(self,data_dict): + dic = { + # 'site_gid': self.site_gid, + # 'nsrdb_lat':self.nsrdb_latitude, + # 'nsrdb_lon':self.nsrdb_longitude, + 'tz' : self.time_zone, + 'elev' : self.elevation, + 'lat' : self.nsrdb_latitude, + 'lon' : self.nsrdb_longitude, + 'year' : self.year_arr, + 'month' : self.month_arr, + 'day' : self.day_arr, + 'hour' : self.hour_arr, + 'minute' : self.minute_arr, + 'dn' : self.dni_arr, + 'df' : self.dhi_arr, + 'gh' : self.ghi_arr, + 'wspd' : self.wspd_arr, + 'tdry' : self.tdry_arr, + 'pres' : self.pres_arr, + 'tdew' : self.tdew_arr + } + self._data = dic \ No newline at end of file diff --git a/hopp/simulation/technologies/resource/wind_toolkit_data.py b/hopp/simulation/technologies/resource/wind_toolkit_data.py new file mode 100644 index 000000000..44743ef66 --- /dev/null +++ b/hopp/simulation/technologies/resource/wind_toolkit_data.py @@ -0,0 +1,165 @@ +from rex import WindX +from rex.sam_resource import SAMResource +import numpy as np +from typing import Optional, Union +from pathlib import Path +import os +from hopp.simulation.technologies.resource.resource import Resource +WTK_V10_BASE = "/kfs2/datasets/WIND/conus/v1.0.0/wtk_conus_" +WTK_V11_BASE = "/kfs2/datasets/WIND/conus/v1.1.0/wtk_conus_" +class HPCWindData(Resource): + def __init__( + self, + lat: float, + lon: float, + year: int, + hub_height_meters: float, + wtk_source_path: Union[str,Path] = "", + filepath: str = "", + **kwargs + ): + """ + Input: + lat (float): site latitude + lon (float): site longitude + resource_year (int): year to get resource data for + wtk_source_path (str): directory of wind resource data on HPC + filepath (str): filepath for wind toolkit h5 file on HPC + """ + + + + self.latitude = lat + self.longitude = lon + self.year = year + super().__init__(lat, lon, year) + + self.hub_height_meters = hub_height_meters + self.allowed_hub_heights_meters = [10, 40, 60, 80, 100, 120, 140, 160, 200] + self.data_hub_heights = self.calculate_heights_to_download() + + + if filepath == "" and wtk_source_path=="": + if self.year < 2014: + self.wtk_file = WTK_V10_BASE + "{}.h5".format(self.year) + # wtk_file = '/datasets/WIND/conus/v1.0.0/wtk_conus_{year}.h5'.format(year=self.year) + elif self.year == 2014: + self.wtk_file = WTK_V11_BASE + "{}.h5".format(self.year) + elif filepath != "" and wtk_source_path == "": + self.wtk_file = filepath + elif filepath=="" and wtk_source_path !="": + self.wtk_file = os.path.join(str(wtk_source_path),"wtk_conus_{}.h5".format(self.year)) + else: + if self.year < 2014: + self.wtk_file = WTK_V10_BASE + "{}.h5".format(self.year) + # wtk_file = '/datasets/WIND/conus/v1.0.0/wtk_conus_{year}.h5'.format(year=self.year) + elif self.year == 2014: + self.wtk_file = WTK_V11_BASE + "{}.h5".format(self.year) + + # self.extract_resource() + self.download_resource() + self.format_data() + + # self.data = {'heights': [float(h) for h in self.data_hub_heights for i in range(4)], + # 'fields': [1, 2, 3, 4] * len(self.data_hub_heights), + # 'data': self.combined_data + # } + + + + def calculate_heights_to_download(self): + """ + Given the system hub height, and the available hubheights from WindToolkit, + determine which heights to download to bracket the hub height + """ + hub_height_meters = self.hub_height_meters + + # evaluate hub height, determine what heights to download + heights = [hub_height_meters] + if hub_height_meters not in self.allowed_hub_heights_meters: + height_low = self.allowed_hub_heights_meters[0] + height_high = self.allowed_hub_heights_meters[-1] + for h in self.allowed_hub_heights_meters: + if h < hub_height_meters: + height_low = h + elif h > hub_height_meters: + height_high = h + break + heights[0] = height_low + heights.append(height_high) + + return heights + + # def extract_resource(self): + def download_resource(self): + # Define file to download from + # NOTE: Current setup of files on HPC WINDToolkit v1.0.0 = 2007-2013, v1.1.0 = 2014 + + + # Open file with rex WindX object + with WindX(self.wtk_file, hsds=False) as f: + # get gid of location closest to given lat/lon coordinates and timezone offset + site_gid = f.lat_lon_gid((self.latitude, self.longitude)) + time_zone = f.meta['timezone'].iloc[site_gid] + + # instantiate temp dictionary to hold each attributes dataset + self.wind_dict = {} + # loop through hub heights to download, capture datasets + # NOTE: datasets are not auto shifted by timezone offset -> wrap extraction in SAMResource.roll_timeseries(input_array, timezone, #steps in an hour=1) to roll timezones + # NOTE: pressure datasets unit = Pa, convert to atm via division by 101325 + for h in self.data_hub_heights: + self.wind_dict['temperature_{height}m_arr'.format(height=h)] = SAMResource.roll_timeseries((f['temperature_{height}m'.format(height=h), :, site_gid]), time_zone, 1) + self.wind_dict['pressure_{height}m_arr'.format(height=h)] = SAMResource.roll_timeseries((f['pressure_{height}m'.format(height=h), :, site_gid]/101325), time_zone, 1) + self.wind_dict['windspeed_{height}m_arr'.format(height=h)] = SAMResource.roll_timeseries((f['windspeed_{height}m'.format(height=h), :, site_gid]), time_zone, 1) + self.wind_dict['winddirection_{height}m_arr'.format(height=h)] = SAMResource.roll_timeseries((f['winddirection_{height}m'.format(height=h), :, site_gid]), time_zone, 1) + + self.site_gid = site_gid + def format_data(self): + # Remove data from feb29 on leap years + if (self.year % 4) == 0: + feb29 = np.arange(1416,1440) + for key, value in self.wind_dict.items(): + self.wind_dict[key] = np.delete(value, feb29) + + # round to desired precision and concatenate data into format needed for data dictionary + if len(self.data_hub_heights) == 2: + # NOTE: Unsure if SAM/PySAM is sensitive to data types ie: floats with long precision vs to 2 or 3 decimals. If not sensitive, can remove following 8 lines of code to increase computational efficiency + self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=1) + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=2) + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=3) + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=1) + self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[1])] = np.round((self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[1])]), decimals=1) + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[1])] = np.round((self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[1])]), decimals=2) + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[1])] = np.round((self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[1])]), decimals=3) + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[1])] = np.round((self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[1])]), decimals=1) + # combine all data into one 2D list + self.combined_data = [list(a) for a in zip(self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[1])], + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[1])], + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[1])], + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[1])])] + + elif len(self.data_hub_heights) == 1: + # NOTE: Unsure if SAM/PySAM is sensitive to data types ie: floats with long precision vs to 2 or 3 decimals. If not sensitive, can remove following 4 lines of code to increase computational efficiency + self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=1) + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=2) + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=3) + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])] = np.round((self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])]), decimals=1) + # combine all data into one 2D list + self.combined_data = [list(a) for a in zip(self.wind_dict['temperature_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['pressure_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['windspeed_{h}m_arr'.format(h=self.data_hub_heights[0])], + self.wind_dict['winddirection_{h}m_arr'.format(h=self.data_hub_heights[0])])] + self.data = self.combined_data + + @Resource.data.setter + def data(self, data_file): + dic = { + 'heights': [float(h) for h in self.data_hub_heights for i in range(4)], + 'fields': [1, 2, 3, 4] * len(self.data_hub_heights), + 'data': data_file #self.combined_data + } + self._data = dic \ No newline at end of file diff --git a/hopp/simulation/technologies/sites/site_info.py b/hopp/simulation/technologies/sites/site_info.py index b1c40ed4d..3d16ed2ab 100644 --- a/hopp/simulation/technologies/sites/site_info.py +++ b/hopp/simulation/technologies/sites/site_info.py @@ -17,7 +17,9 @@ SolarResource, WindResource, WaveResource, - ElectricityPrices + ElectricityPrices, + HPCWindData, + HPCSolarData, ) from hopp.tools.layout.plot_tools import plot_shape from hopp.utilities.log import hybrid_logger as logger @@ -29,6 +31,7 @@ from hopp.simulation.base import BaseClass from hopp.utilities.validators import contains +from hopp import ROOT_DIR def plot_site(verts, plt_style, labels): for i in range(len(verts)): if i == 0: @@ -65,6 +68,11 @@ class SiteInfo(BaseClass): wind_resource_file: Union[Path, str] = field(default="", converter=resource_file_converter) wave_resource_file: Union[Path, str] = field(default="", converter=resource_file_converter) grid_resource_file: Union[Path, str] = field(default="", converter=resource_file_converter) + + path_resource: Optional[Union[Path, str]] = field(default=ROOT_DIR / "simulation" / "resource_files") + wtk_source_path: Optional[Union[Path,str]] = field(default = "") + nsrdb_source_path: Optional[Union[Path,str]] = field(default = "") + hub_height: hopp_float_type = field(default=97., converter=hopp_float_type) capacity_hours: NDArray = field(default=[], converter=converter(bool)) desired_schedule: NDArrayFloat = field(default=[], converter=converter()) @@ -73,6 +81,7 @@ class SiteInfo(BaseClass): solar: bool = field(default=True) wind: bool = field(default=True) wave: bool = field(default=False) + renewable_resource_origin: str = field(default="API", validator=contains(["API", "HPC"])) wind_resource_origin: str = field(default="WTK", validator=contains(["WTK", "TAP"])) # Set in post init hook @@ -81,8 +90,8 @@ class SiteInfo(BaseClass): lon: hopp_float_type = field(init=False) year: int = field(init=False, default=2012) tz: Optional[int] = field(init=False, default=None) - solar_resource: Optional[SolarResource] = field(init=False, default=None) - wind_resource: Optional[WindResource] = field(init=False, default=None) + solar_resource: Optional[Union[SolarResource,HPCSolarData]] = field(default=None) + wind_resource: Optional[Union[WindResource,HPCWindData]] = field(default=None) wave_resoure: Optional[WaveResource] = field(init=False, default=None) elec_prices: Optional[ElectricityPrices] = field(init=False, default=None) n_periods_per_day: int = field(init=False) @@ -135,7 +144,11 @@ def __attrs_post_init__(self): self.tz = data['tz'] if self.solar: - self.solar_resource = SolarResource(data['lat'], data['lon'], data['year'], filepath=self.solar_resource_file) + if self.solar_resource is None: + if self.renewable_resource_origin=="API": + self.solar_resource = SolarResource(data['lat'], data['lon'], data['year'], path_resource=self.path_resource, filepath=self.solar_resource_file) + else: + self.solar_resource = HPCSolarData(data['lat'], data['lon'], data['year'],nsrdb_source_path = self.nsrdb_source_path, filepath=self.solar_resource_file) self.n_timesteps = len(self.solar_resource.data['gh']) // 8760 * 8760 if self.wave: self.wave_resource = WaveResource(data['lat'], data['lon'], data['year'], filepath = self.wave_resource_file) @@ -143,8 +156,13 @@ def __attrs_post_init__(self): if self.wind: # TODO: allow hub height to be used as an optimization variable - self.wind_resource = WindResource(data['lat'], data['lon'], data['year'], wind_turbine_hub_ht=self.hub_height, - filepath=self.wind_resource_file, source=self.wind_resource_origin) + if self.wind_resource is None: + if self.renewable_resource_origin=="API": + self.wind_resource = WindResource(data['lat'], data['lon'], data['year'], wind_turbine_hub_ht=self.hub_height, + path_resource=self.path_resource, filepath=self.wind_resource_file, source=self.wind_resource_origin) + else: + self.wind_resource = HPCWindData(data['lat'], data['lon'], data['year'], hub_height_meters=self.hub_height, + wtk_source_path=self.wtk_source_path, filepath=self.wind_resource_file) n_timesteps = len(self.wind_resource.data['data']) // 8760 * 8760 if self.n_timesteps is None: self.n_timesteps = n_timesteps diff --git a/pyproject.toml b/pyproject.toml index 76156a499..94a822091 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,7 +48,8 @@ dependencies = [ "attrs", "utm", "pyyaml-include", - "profast" + "profast", + "NREL-rex" ] keywords = [ "python3",