Enable independent HC and RC in wind downscaling#2324
Conversation
… and function names, etc
82e18ff to
41f9625
Compare
brhooper
left a comment
There was a problem hiding this comment.
Thanks @mo-AliceLake , this looks to have made a lot of progress towards allowing independent application of roughness and height corrections for wind downscaling.
I've added a few comments suggesting ways in which I think the code might be improved. I haven't gotten as far as reviewing the unit tests yet, I hope to review them tomorrow. I have noted that the acceptance tests need updating following changes to the names of things.
| Returns: | ||
| 3D float32 array of wind speed after applying RC and HC. |
There was a problem hiding this comment.
This statement needs updating.
| def _mask_missing_data(self, height_above_orog, wspeed_original): | ||
| """Return a boolean mask: True where either RC or HC may be applied. |
There was a problem hiding this comment.
Have the type hints been removed here on purpose?
| if wspeed_original.ndim == 3: | ||
| missing_w = (wspeed_original == RMDI).any(axis=2) | ||
| else: | ||
| missing_w = wspeed_original == RMDI |
There was a problem hiding this comment.
| if wspeed_original.ndim == 3: | |
| missing_w = (wspeed_original == RMDI).any(axis=2) | |
| else: | |
| missing_w = wspeed_original == RMDI | |
| missing_w = wspeed_original == RMDI | |
| if wspeed_original.ndim == 3: | |
| missing_w = missing_w.any(axis=2) |
I think this would work.
|
|
||
| return valid | ||
|
|
||
| def do_rc(self, height_above_orog, wspeed_original): |
There was a problem hiding this comment.
This function deserves a doc-string and type-hints.
| height_above_orog, wspeed_original, mask_rc | ||
| ) | ||
|
|
||
| def do_hc(self, height_above_orog, wspeed_original): |
There was a problem hiding this comment.
This function deserves a doc-string and type-hints.
| # Mask where missing data in height and wind fields | ||
| valid = self._mask_missing_data(height_above_orog, wspeed_original) | ||
| mask_rc = np.copy(self.rc_mask) | ||
| mask_rc[~valid] = False | ||
|
|
||
| return self.calc_roughness_correction( | ||
| height_above_orog, wspeed_original, mask_rc | ||
| ) |
There was a problem hiding this comment.
Slightly reduces data copying.
| # Mask where missing data in height and wind fields | |
| valid = self._mask_missing_data(height_above_orog, wspeed_original) | |
| mask_rc = np.copy(self.rc_mask) | |
| mask_rc[~valid] = False | |
| return self.calc_roughness_correction( | |
| height_above_orog, wspeed_original, mask_rc | |
| ) | |
| # Mask where missing data in height and wind fields | |
| valid = self._mask_missing_data(height_above_orog, wspeed_original) | |
| valid[~self.rc_mask] = False | |
| return self.calc_roughness_correction( | |
| height_above_orog, wspeed_original, valid | |
| ) |
| @@ -715,17 +741,30 @@ | |||
| uhref_orig, height_above_orog, mask_hc, onemfrac | |||
| ) | |||
There was a problem hiding this comment.
Same suggested change as in do_rc() to remove a data copy step.
| # Mask where missing data in height and wind fields | |
| valid = self._mask_missing_data(height_above_orog, wspeed_original) | |
| valid[~self.hc_mask] = False | |
| # Height correction | |
| # Requires wind speed at the reference height, so interpolate first | |
| z_ref = 1.0 / self.wavenumber | |
| if wspeed_original.ndim == 3 and wspeed_original.shape[2] > 1: | |
| uhref_orig = self._interpolate_wspeed_to_height( | |
| wspeed_original, | |
| height_above_orog, | |
| z_ref, | |
| valid, | |
| ) | |
| else: | |
| # Single level (e.g. 10m wind) so no interpolation possible | |
| uhref_orig = np.copy(wspeed_original) | |
| # HC only where u(h_ref) is positive | |
| valid[uhref_orig <= 0.0] = False | |
| # Setting this value to 1, is equivalent to setting the | |
| # Bessel function to 1. (Friedrich, 2016) | |
| # Example usage if the Bessel function was not set to 1 is: | |
| # onemfrac = 1.0 - BfuncFrac(nx,ny,nz,heightvec,z_0,waveno, Ustar, UI) | |
| onemfrac = 1.0 | |
| hc_add = self._calc_height_corr( | |
| uhref_orig, height_above_orog, valid, onemfrac | |
| ) |
| model_z0_cube: | ||
| 2D vegetative roughness length (m), representing drag from vegetation | ||
| and land cover. Controls the near-surface wind profile. | ||
| Historically static, but may now be time-varying (e.g. from StaGE). | ||
|
|
||
| height_levels_cube: | ||
| 1D or 3D height levels of the input wind field (m). | ||
|
|
||
| mode: | ||
| Which correction(s) to apply: "hc_and_rc" (default), "hc", or "rc". |
There was a problem hiding this comment.
I briefly wondered whether the mode input could be removed in favour of assuming providing model_z0_cube implied that a roughness correction should be applied and similarly that providing height_levels_cube implied a height correction should be applied.
I see now, however, that height_levels_cube is not actually necessary for the height correction, so this wouldn't work. Leaving this comment here as evidence that, sometimes, I do try to think during reviews. There's nothing to action.
| if np.isnan(zwp): | ||
| # Reorder wind cube so dimensions are consistently (y, x, [, t]) | ||
| if np.isnan(twp): | ||
| input_cube.transpose([ywp, xwp]) | ||
| else: | ||
| input_cube.transpose([ywp, xwp, twp]) | ||
| else: | ||
| input_cube.transpose([ywp, xwp, zwp, twp]) | ||
| # Reorder wind cube so dimensions are consistently (y, x, z [, t]) | ||
| if np.isnan(twp): | ||
| input_cube.transpose([ywp, xwp, zwp]) | ||
| else: | ||
| input_cube.transpose([ywp, xwp, zwp, twp]) |
There was a problem hiding this comment.
I think this can be made a bit shorter.
| if np.isnan(zwp): | |
| # Reorder wind cube so dimensions are consistently (y, x, [, t]) | |
| if np.isnan(twp): | |
| input_cube.transpose([ywp, xwp]) | |
| else: | |
| input_cube.transpose([ywp, xwp, twp]) | |
| else: | |
| input_cube.transpose([ywp, xwp, zwp, twp]) | |
| # Reorder wind cube so dimensions are consistently (y, x, z [, t]) | |
| if np.isnan(twp): | |
| input_cube.transpose([ywp, xwp, zwp]) | |
| else: | |
| input_cube.transpose([ywp, xwp, zwp, twp]) | |
| # Reorder wind cube so dimensions are consistently (y, x), (y, x, [, t]) or | |
| # (y, x, z [, t]) depending on what dimensions are present. | |
| coord_order = [ywp, xwp] | |
| for coord in [zwp, twp]: | |
| if not np.isnan(coord): | |
| coord_order.append(coord) | |
| input_cube.transpose(coord_order) |
| dims = [ywp, xwp] | ||
| if not np.isnan(zwp): | ||
| dims.append(zwp) | ||
| input_dims = dims.copy() | ||
| output_dims = dims.copy() | ||
| if not np.isnan(twp): | ||
| input_dims.append(twp) | ||
| output_dims.insert(0, twp) | ||
| input_cube.transpose(np.argsort(input_dims)) | ||
| output_cube.transpose(np.argsort(output_dims)) |
There was a problem hiding this comment.
This might be slightly more readable.
| dims = [ywp, xwp] | |
| if not np.isnan(zwp): | |
| dims.append(zwp) | |
| input_dims = dims.copy() | |
| output_dims = dims.copy() | |
| if not np.isnan(twp): | |
| input_dims.append(twp) | |
| output_dims.insert(0, twp) | |
| input_cube.transpose(np.argsort(input_dims)) | |
| output_cube.transpose(np.argsort(output_dims)) | |
| input_dims = [ywp, xwp, zwp, twp] | |
| output_dims = [twp, ywp, xwp, zwp] | |
| for coord in [zwp, twp]: | |
| if np.isnan(coord): | |
| input_dims.remove(coord) | |
| output_dims.remove(coord) | |
| input_cube.transpose(np.argsort(input_dims)) | |
| output_cube.transpose(np.argsort(output_dims)) |
brhooper
left a comment
There was a problem hiding this comment.
Thanks @mo-AliceLake, one additional comment that I should have added yesterday. I think the unit tests look fine.
| mode: | ||
| Which correction(s) to apply: "hc_and_rc" (default), "hc", or "rc". |
There was a problem hiding this comment.
One more thought: it might be worth explicitly defining what each of rc, hc, and hc_and_rc actually do.
Description
This PR restructures the WindTerrainAdjustment plugin (previously RoughnessCorrection) to provide clean, explicit control over applying Roughness Correction (RC), Height Correction (HC), or their combined behaviour. The goal is to make HC‑only and RC‑only configurations explicitly selectable, allowing clearer scientific experimentation and operational configurations (e.g. HC‑only over UK-domain, etc).
Previously, HC and RC were tightly coupled within a single routine. Running only one of the corrections was possible but relied on implicit side‑effects (e.g. suppressing RC by passing z0=None), which was neither obvious nor robust. This PR replaces those implicit behaviours with a clear interface.
Previously, HC and RC were tightly coupled inside a single routine (do_rc_hc_all). This made it difficult to:
Summary of Changes
Introduced explicit correction modes:
do_rc_hc_allroutine into three:do_hc()do_rc()do_hc_and_rc()(chained, preserving legacy HC+RC behaviour)modeargument to the plugin constructor, allowing users to select:"hc_and_rc"(default, identical to previous behaviour)"hc""rc"Clarified RC dependency on z0:
In the new interface, requesting RC (mode="rc" or "hc_and_rc") now requires a valid roughness‑length (z0) field. If it is missing, the plugin raises a clear ValueError.
This replaces the old implicit behaviour where suppressing RC was only achievable indirectly by passing z0=None. The new behaviour is explicit and predictable.
Unified masking and validation
Missing‑data, sea‑point, and invalid‑height masking has been consolidated and made consistent across RC, HC, and combined modes.
Testing