Skip to content

Commit

Permalink
Merge pull request #155 from NREL/bnb/log_interp
Browse files Browse the repository at this point in the history 
Bnb/log interp
  • Loading branch information
@bnb32
bnb32 authored Aug 18, 2023
2 parents 0146035 + 3f056c6 commit 117816a
Show file tree
Hide file tree
Showing 20 changed files with 4,364 additions and 1,213 deletions.
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
matplotlib>=3.1
NREL-rex>=0.2.82
NREL-phygnn>=0.0.23
NREL-rev>=0.6.6
NREL-rev<0.8.0
NREL-farms>=1.0.4
pytest>=5.2
pillow
Expand Down
144 changes: 109 additions & 35 deletions sup3r/bias/bias_transforms.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# -*- coding: utf-8 -*-
"""Bias correction transformation functions."""
import logging
import os
from warnings import warn

import numpy as np
Expand All @@ -10,6 +11,62 @@
logger = logging.getLogger(__name__)


def get_spatial_bc_factors(lat_lon, feature_name, bias_fp, threshold=0.1):
"""Get bc factors (scalar/adder) for the given feature for the given
domain (specified by lat_lon).
Parameters
----------
lat_lon : ndarray
Array of latitudes and longitudes for the domain to bias correct
(n_lats, n_lons, 2)
feature_name : str
Name of feature that is being corrected. Datasets with names
"{feature_name}_scalar" and "{feature_name}_adder" will be retrieved
from bias_fp.
bias_fp : str
Filepath to bias correction file from the bias calc module. Must have
datasets "{feature_name}_scalar" and "{feature_name}_adder" that are
the full low-resolution shape of the forward pass input that will be
sliced using lr_padded_slice for the current chunk.
threshold : float
Nearest neighbor euclidean distance threshold. If the coordinates are
more than this value away from the bias correction lat/lon, an error is
raised.
"""
dset_scalar = f'{feature_name}_scalar'
dset_adder = f'{feature_name}_adder'
with Resource(bias_fp) as res:
lat = np.expand_dims(res['latitude'], axis=-1)
lon = np.expand_dims(res['longitude'], axis=-1)
lat_lon_bc = np.dstack((lat, lon))
diff = lat_lon_bc - lat_lon[:1, :1]
diff = np.hypot(diff[..., 0], diff[..., 1])
idy, idx = np.where(diff == diff.min())
slice_y = slice(idy[0], idy[0] + lat_lon.shape[0])
slice_x = slice(idx[0], idx[0] + lat_lon.shape[1])

if diff.min() > threshold:
msg = (
'The DataHandler top left coordinate of {} '
'appears to be {} away from the nearest '
'bias correction coordinate of {} from {}. '
'Cannot apply bias correction.'.format(
lat_lon,
diff.min(),
lat_lon_bc[idy, idx],
os.path.basename(bias_fp),
)
)
logger.error(msg)
raise RuntimeError(msg)

assert dset_scalar in res.dsets and dset_adder in res.dsets
scalar = res[dset_scalar, slice_y, slice_x]
adder = res[dset_adder, slice_y, slice_x]
return scalar, adder


def global_linear_bc(input, scalar, adder, out_range=None):
"""Bias correct data using a simple global *scalar +adder method.
Expand Down Expand Up @@ -37,8 +94,15 @@ def global_linear_bc(input, scalar, adder, out_range=None):
return out


def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
out_range=None, smoothing=0):
def local_linear_bc(
input,
lat_lon,
feature_name,
bias_fp,
lr_padded_slice,
out_range=None,
smoothing=0,
):
"""Bias correct data using a simple annual (or multi-year) *scalar +adder
method on a site-by-site basis.
Expand All @@ -47,6 +111,9 @@ def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
input : np.ndarray
Sup3r input data to be bias corrected, assumed to be 3D with shape
(spatial, spatial, temporal) for a single feature.
lat_lon : ndarray
Array of latitudes and longitudes for the domain to bias correct
(n_lats, n_lons, 2)
feature_name : str
Name of feature that is being corrected. Datasets with names
"{feature_name}_scalar" and "{feature_name}_adder" will be retrieved
Expand Down Expand Up @@ -77,12 +144,7 @@ def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
out = input * scalar + adder
"""

scalar = f'{feature_name}_scalar'
adder = f'{feature_name}_adder'
with Resource(bias_fp) as res:
scalar = res[scalar]
adder = res[adder]

scalar, adder = get_spatial_bc_factors(lat_lon, feature_name, bias_fp)
# 3D bias correction factors have seasonal/monthly correction in last axis
if len(scalar.shape) == 3 and len(adder.shape) == 3:
scalar = scalar.mean(axis=-1)
Expand All @@ -94,8 +156,10 @@ def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
adder = adder[spatial_slice]

if np.isnan(scalar).any() or np.isnan(adder).any():
msg = ('Bias correction scalar/adder values had NaNs for "{}" from: {}'
.format(feature_name, bias_fp))
msg = (
'Bias correction scalar/adder values had NaNs for '
f'"{feature_name}" from: {bias_fp}'
)
logger.warning(msg)
warn(msg)

Expand All @@ -107,12 +171,12 @@ def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,

if smoothing > 0:
for idt in range(scalar.shape[-1]):
scalar[..., idt] = gaussian_filter(scalar[..., idt],
smoothing,
mode='nearest')
adder[..., idt] = gaussian_filter(adder[..., idt],
smoothing,
mode='nearest')
scalar[..., idt] = gaussian_filter(
scalar[..., idt], smoothing, mode='nearest'
)
adder[..., idt] = gaussian_filter(
adder[..., idt], smoothing, mode='nearest'
)

out = input * scalar + adder
if out_range is not None:
Expand All @@ -122,9 +186,17 @@ def local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
return out


def monthly_local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
time_index, temporal_avg=True, out_range=None,
smoothing=0):
def monthly_local_linear_bc(
input,
lat_lon,
feature_name,
bias_fp,
lr_padded_slice,
time_index,
temporal_avg=True,
out_range=None,
smoothing=0,
):
"""Bias correct data using a simple monthly *scalar +adder method on a
site-by-site basis.
Expand All @@ -133,6 +205,9 @@ def monthly_local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
input : np.ndarray
Sup3r input data to be bias corrected, assumed to be 3D with shape
(spatial, spatial, temporal) for a single feature.
lat_lon : ndarray
Array of latitudes and longitudes for the domain to bias correct
(n_lats, n_lons, 2)
feature_name : str
Name of feature that is being corrected. Datasets with names
"{feature_name}_scalar" and "{feature_name}_adder" will be retrieved
Expand Down Expand Up @@ -172,12 +247,7 @@ def monthly_local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
out : np.ndarray
out = input * scalar + adder
"""

scalar = f'{feature_name}_scalar'
adder = f'{feature_name}_adder'
with Resource(bias_fp) as res:
scalar = res[scalar]
adder = res[adder]
scalar, adder = get_spatial_bc_factors(lat_lon, feature_name, bias_fp)

assert len(scalar.shape) == 3, 'Monthly bias correct needs 3D scalars'
assert len(adder.shape) == 3, 'Monthly bias correct needs 3D adders'
Expand All @@ -199,25 +269,29 @@ def monthly_local_linear_bc(input, feature_name, bias_fp, lr_padded_slice,
scalar = np.repeat(scalar, input.shape[-1], axis=-1)
adder = np.repeat(adder, input.shape[-1], axis=-1)
if len(time_index.month.unique()) > 2:
msg = ('Bias correction method "monthly_local_linear_bc" was used '
'with temporal averaging over a time index with >2 months.')
msg = (
'Bias correction method "monthly_local_linear_bc" was used '
'with temporal averaging over a time index with >2 months.'
)
warn(msg)
logger.warning(msg)

if np.isnan(scalar).any() or np.isnan(adder).any():
msg = ('Bias correction scalar/adder values had NaNs for "{}" from: {}'
.format(feature_name, bias_fp))
msg = (
'Bias correction scalar/adder values had NaNs for '
f'"{feature_name}" from: {bias_fp}'
)
logger.warning(msg)
warn(msg)

if smoothing > 0:
for idt in range(scalar.shape[-1]):
scalar[..., idt] = gaussian_filter(scalar[..., idt],
smoothing,
mode='nearest')
adder[..., idt] = gaussian_filter(adder[..., idt],
smoothing,
mode='nearest')
scalar[..., idt] = gaussian_filter(
scalar[..., idt], smoothing, mode='nearest'
)
adder[..., idt] = gaussian_filter(
adder[..., idt], smoothing, mode='nearest'
)

out = input * scalar + adder
if out_range is not None:
Expand Down
33 changes: 23 additions & 10 deletions sup3r/pipeline/forward_pass.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -849,6 +849,14 @@ def get_time_index(self, file_paths, max_workers=None, **kwargs):
----------
file_paths : list
List of file paths for source data
max_workers : int | None
Number of workers to use to extract the time index from the given
files. This is used when a large number of single timestep netcdf
files is provided.
**kwargs : dict
Dictionary of kwargs passed to the resource handler opening the
given file_paths. For netcdf files this is xarray.open_mfdataset().
For h5 files this is usually rex.Resource().
Returns
-------
Expand Down Expand Up @@ -1061,7 +1069,8 @@ def __init__(self, strategy, chunk_index=0, node_index=0):
self.data_handler.load_cached_data()
self.input_data = self.data_handler.data

self.input_data = self.bias_correct_source_data(self.input_data)
self.input_data = self.bias_correct_source_data(
self.input_data, self.strategy.lr_lat_lon)

exo_s_en = self.exo_kwargs.get('s_enhancements', None)
out = self.pad_source_data(self.input_data, self.pad_width,
Expand Down Expand Up @@ -1400,7 +1409,7 @@ def pad_source_data(input_data, pad_width, exo_data,

return out, exo_data

def bias_correct_source_data(self, data):
def bias_correct_source_data(self, data, lat_lon):
"""Bias correct data using a method defined by the bias_correct_method
input to ForwardPassStrategy
Expand All @@ -1409,6 +1418,10 @@ def bias_correct_source_data(self, data):
data : np.ndarray
Any source data to be bias corrected, with the feature channel in
the last axis.
lat_lon : np.ndarray
Latitude longitude array for the given data. Used to get the
correct bc factors for the appropriate domain.
(n_lats, n_lons, 2)
Returns
-------
Expand All @@ -1433,7 +1446,8 @@ def bias_correct_source_data(self, data):
'using function: {} with kwargs: {}'
.format(feature, idf, method, feature_kwargs))

data[..., idf] = method(data[..., idf], **feature_kwargs)
data[..., idf] = method(data[..., idf], lat_lon,
**feature_kwargs)

return data

Expand Down Expand Up @@ -1716,8 +1730,8 @@ def _single_proc_run(cls, strategy, node_index, chunk_index):

@classmethod
def run(cls, strategy, node_index):
"""This routine runs forward passes on all spatiotemporal chunks for
the given node index.
"""Runs forward passes on all spatiotemporal chunks for the given node
index.
Parameters
----------
Expand All @@ -1738,8 +1752,8 @@ def run(cls, strategy, node_index):

@classmethod
def _run_serial(cls, strategy, node_index):
"""This routine runs forward passes, on all spatiotemporal chunks for
the given node index, in serial.
"""Runs forward passes, on all spatiotemporal chunks for the given node
index, in serial.
Parameters
----------
Expand Down Expand Up @@ -1772,9 +1786,8 @@ def _run_serial(cls, strategy, node_index):

@classmethod
def _run_parallel(cls, strategy, node_index):
"""This routine runs forward passes, on all spatiotemporal chunks for
the given node index, with data extraction and forward pass routines in
parallel.
"""Runs forward passes, on all spatiotemporal chunks for the given node
index, with data extraction and forward pass routines in parallel.
Parameters
----------
Expand Down
Loading

0 comments on commit 117816a

Please sign in to comment.