Skip to content

Commit

Permalink
Add user controlled perturber features and add a relevant unittest
Browse files Browse the repository at this point in the history
  • Loading branch information
@SorooshMani-NOAA
SorooshMani-NOAA committed Oct 4, 2024
1 parent 587a57c commit d834b8a
Show file tree
Hide file tree
Showing 22 changed files with 378 additions and 69 deletions.
118 changes: 49 additions & 69 deletions ensembleperturbation/perturbation/atcf.py
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
from abc import ABC
import concurrent.futures
from collections import namedtuple
from copy import copy
from difflib import get_close_matches
from datetime import datetime, timedelta
from enum import Enum
from enum import Enum, Flag, auto
from glob import glob
import json
from math import exp, inf, sqrt
import os
from os import PathLike
from pathlib import Path
from random import gauss, uniform
from typing import Dict, List, Mapping, Union
from typing import Dict, List, Mapping, Union, Optional
import warnings

import chaospy
Expand Down Expand Up @@ -712,20 +713,22 @@ def storm_errors(self, data_frame: DataFrame, loc_index: int = 0) -> DataFrame:
return self.historical_forecast_errors[storm_classification]


# IsotachRadius parent class for isotach radii dependent
# on the changes to Rmax, Vmax and track
class IsotachRadius(VortexPerturbedVariable):
IsotachQuadrant = namedtuple('IsotachQuadrant', ['angle', 'column'])

name = 'isotach_radius'
perturbation_type = None

def __init__(self):
super().__init__(
lower_bound=5,
upper_bound=400,
historical_forecast_errors={},
unit=units.nautical_mile,
)
class IsotachRadiusQuadrant(Enum):
SWQ = IsotachQuadrant(225, 'isotach_radius_for_SWQ')
SEQ = IsotachQuadrant(135, 'isotach_radius_for_SEQ')
NWQ = IsotachQuadrant(315, 'isotach_radius_for_NWQ')
NEQ = IsotachQuadrant(45, 'isotach_radius_for_NEQ')


# IsotachRadius adjuster dependent on the changes to Rmax, Vmax and track
class IsotachRadius:

def __init__(self, rad_quad_info: IsotachRadiusQuadrant):
self.quadrant_angle = rad_quad_info.value.angle
self.column = rad_quad_info.value.column

def find_GAHM_parameters(self, B, Ro_inv) -> Quantity:
"""
Expand Down Expand Up @@ -831,7 +834,7 @@ def ts_factor(isotach_rad, Rmax):
dataframe['isotach_radius'].values
* MaximumSustainedWindSpeed.unit # LC12 factor and 20-deg CCW rotation
- ts_factor(
dataframe[self.name].values, dataframe[RadiusOfMaximumWinds.name].values
dataframe[self.column].values, dataframe[RadiusOfMaximumWinds.name].values
)
* translation_speed
* numpy.sin(numpy.deg2rad(self.quadrant_angle + 20))
Expand All @@ -849,7 +852,7 @@ def ts_factor(isotach_rad, Rmax):
# through both the isotach and the Rmax value
return Vmax, Vr, Rmax, Ro_inv, f, B

def perturb(
def adjust(
self, vortex_dataframe: DataFrame, original_dataframe, inplace: bool = False,
) -> DataFrame:
""" Compute new isotach radius based on Holland profile using perturbed Rmax/Vmax/speed values"""
Expand All @@ -862,7 +865,7 @@ def perturb(
dataframe=vortex_dataframe
)
# determine original isotach_rad and Rrat
isotach_rad = original_dataframe[self.name].values * RadiusOfMaximumWinds.unit
isotach_rad = original_dataframe[self.column].values * RadiusOfMaximumWinds.unit
isotach_rad[isotach_rad == 0] = numpy.nan
Rrat = Rmax_old / isotach_rad # [nmi/nmi] = []
Rrat[abs(Rrat - 1.0) < 1e-3] = 0.999 # ensure not exactly 1
Expand All @@ -888,51 +891,10 @@ def perturb(
)
isotach_rad_new = isotach_rad_new.magnitude
isotach_rad_new[numpy.isnan(isotach_rad_new)] = 0
vortex_dataframe[self.name] = isotach_rad_new
vortex_dataframe[self.column] = isotach_rad_new
return vortex_dataframe


# IsotachRadius subclass for each quadrant
class IsotachRadiusSWQ(IsotachRadius):
"""
``isotach_radius_for_SWQ``
"""

name = 'isotach_radius_for_SWQ'
perturbation_type = None
quadrant_angle = 225 # degrees


class IsotachRadiusSEQ(IsotachRadius):
"""
``isotach_radius_for_SEQ``
"""

name = 'isotach_radius_for_SEQ'
perturbation_type = None
quadrant_angle = 135 # degrees


class IsotachRadiusNWQ(IsotachRadius):
"""
``isotach_radius_for_NWQ``
"""

name = 'isotach_radius_for_NWQ'
perturbation_type = None
quadrant_angle = 315 # degrees


class IsotachRadiusNEQ(IsotachRadius):
"""
``isotach_radius_for_NEQ``
"""

name = 'isotach_radius_for_NEQ'
perturbation_type = None
quadrant_angle = 45 # degrees


class CrossTrack(VortexPerturbedVariable):
"""
``cross_track`` represents a perpendicular offset of the tropical cyclone center track, accomplished by moving each forecast time left or right perpedicular to the track line.
Expand Down Expand Up @@ -1283,6 +1245,12 @@ def perturb(
return vortex_dataframe


class PerturberFeatures(Flag):
NONE = 0
ISOTACH_ADJUSTMENT = auto()



class VortexPerturber:
"""
``VortexPerturber`` takes an ATCF track from an input storm and perturbs it based on several variables (of the class ``VortexPerturbedVariable``)
Expand Down Expand Up @@ -1332,6 +1300,7 @@ def __init__(
end_date: datetime = None,
file_deck: ATCF_FileDeck = None,
advisories: List[str] = None,
features: Optional[PerturberFeatures] = PerturberFeatures.ISOTACH_ADJUSTMENT,
):
"""
:param storm: NHC storm code, for instance `al062018`
Expand All @@ -1346,6 +1315,9 @@ def __init__(
self.__file_deck = None
self.__track = None
self.__previous_configuration = None
if features is None:
features = PerturberFeatures.NONE
self.__features = PerturberFeatures(features)

self.storm = storm
self.start_date = start_date
Expand All @@ -1363,6 +1335,7 @@ def from_file(
end_date: datetime = None,
file_deck: ATCF_FileDeck = None,
advisories: List[str] = None,
features: Optional[PerturberFeatures] = PerturberFeatures.ISOTACH_ADJUSTMENT,
) -> 'VortexPerturber':
"""
build storm perturber from an existing `fort.22` or ATCF file
Expand All @@ -1385,12 +1358,16 @@ def from_file(
track.forecast_time = (
track.data.track_start_time.min() if start_date is None else start_date
)
instance = VortexPerturber.from_track(track)
instance = cls.from_track(track, features)
instance.__filename = filename
return instance

@classmethod
def from_track(cls, track: VortexTrack) -> 'VortexPerturber':
def from_track(
cls,
track: VortexTrack,
features: Optional[PerturberFeatures] = PerturberFeatures.ISOTACH_ADJUSTMENT,
) -> 'VortexPerturber':
"""
build storm perturber from an existing `VortexTrack` object
Expand All @@ -1406,6 +1383,7 @@ def from_track(cls, track: VortexTrack) -> 'VortexPerturber':
end_date=track.end_date,
file_deck=track.file_deck,
advisories=track.advisories,
features=features,
)
instance.track = track
instance.__previous_configuration = {
Expand Down Expand Up @@ -1927,23 +1905,22 @@ def write_perturbed_track(
# Compute potential changes in the central pressure in accordance with Holland B relationship
dataframe[CentralPressure.name] = self.compute_pc_from_Vmax(dataframe)

if RadiusOfMaximumWinds in variables:
# If perturbing for GAHM, use isotach perturbation
if self.__features | PerturberFeatures.ISOTACH_ADJUSTMENT:
# Compute potential changes to r34/50,64 radii at all quadrants in accordance with the GAHM profile
quadrants = [
IsotachRadiusNEQ(),
IsotachRadiusSEQ(),
IsotachRadiusSWQ(),
IsotachRadiusNWQ(),
IsotachRadius(IsotachRadiusQuadrant.NEQ),
IsotachRadius(IsotachRadiusQuadrant.SEQ),
IsotachRadius(IsotachRadiusQuadrant.SWQ),
IsotachRadius(IsotachRadiusQuadrant.NWQ),
]
for quadrant in quadrants:
dataframe = quadrant.perturb(
dataframe = quadrant.adjust(
vortex_dataframe=dataframe,
original_dataframe=self.track.data,
inplace=True,
)
# remove any rows that now have all 0 isotach radii for the 50-kt or 64-kt isotach
quadrant_names = [q.name for q in quadrants]
quadrant_names = [q.column for q in quadrants]
drop_row = (dataframe[quadrant_names] == 0).all(axis=1) & (
dataframe['isotach_radius'].isin([50, 64])
)
Expand Down Expand Up @@ -2184,6 +2161,7 @@ def perturb_tracks(
advisories: List[str] = None,
overwrite: bool = False,
parallel: bool = True,
features: Optional[PerturberFeatures] = PerturberFeatures.ISOTACH_ADJUSTMENT,
):
"""
write a set of perturbed storm tracks
Expand Down Expand Up @@ -2226,6 +2204,7 @@ def perturb_tracks(
end_date=end_date,
file_deck=file_deck,
advisories=advisories,
features=features,
)
else:
raise FileNotFoundError
Expand All @@ -2239,6 +2218,7 @@ def perturb_tracks(
end_date=end_date,
file_deck=file_deck,
advisories=advisories,
features=features,
)

filenames = [directory / 'original.22']
Expand Down
23 changes: 23 additions & 0 deletions tests/data/reference/test_no_isotach_adj/original.22
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
AL, 09, 2021082800, 03, OFCL, 0, 225N, 835W, 70, 985, HU, 34, NEQ, 100, 60, 40, 70, 1010, 0, 20, 85, 0, , 0, JPC, 317, 7, , 1
AL, 09, 2021082800, 03, OFCL, 0, 225N, 835W, 70, 985, HU, 50, NEQ, 50, 30, 0, 30, 1010, 0, 20, 85, 0, , 0, JPC, 317, 7, , 1
AL, 09, 2021082800, 03, OFCL, 0, 225N, 835W, 70, 985, HU, 64, NEQ, 20, 0, 0, 0, 1010, 0, 20, 85, 0, , 0, JPC, 317, 7, , 1
AL, 09, 2021082800, 03, OFCL, 3, 230N, 840W, 70, 989, HU, 34, NEQ, 100, 60, 40, 70, 1010, 0, 18, 85, 0, , 0, JPC, 317, 7, , 2
AL, 09, 2021082800, 03, OFCL, 3, 230N, 840W, 70, 989, HU, 50, NEQ, 50, 30, 0, 30, 1010, 0, 18, 85, 0, , 0, JPC, 317, 7, , 2
AL, 09, 2021082800, 03, OFCL, 3, 230N, 840W, 70, 989, HU, 64, NEQ, 20, 0, 0, 0, 1010, 0, 18, 85, 0, , 0, JPC, 317, 7, , 2
AL, 09, 2021082800, 03, OFCL, 12, 244N, 856W, 85, 973, HU, 34, NEQ, 110, 80, 60, 90, 1010, 0, 16, 105, 0, , 0, JPC, 314, 7, , 3
AL, 09, 2021082800, 03, OFCL, 12, 244N, 856W, 85, 973, HU, 50, NEQ, 50, 40, 20, 40, 1010, 0, 16, 105, 0, , 0, JPC, 314, 7, , 3
AL, 09, 2021082800, 03, OFCL, 12, 244N, 856W, 85, 973, HU, 64, NEQ, 20, 15, 10, 20, 1010, 0, 16, 105, 0, , 0, JPC, 314, 7, , 3
AL, 09, 2021082800, 03, OFCL, 24, 261N, 877W, 105, 954, HU, 34, NEQ, 120, 100, 70, 100, 1010, 0, 13, 130, 0, , 0, JPC, 312, 7, , 4
AL, 09, 2021082800, 03, OFCL, 24, 261N, 877W, 105, 954, HU, 50, NEQ, 60, 50, 30, 40, 1010, 0, 13, 130, 0, , 0, JPC, 312, 7, , 4
AL, 09, 2021082800, 03, OFCL, 24, 261N, 877W, 105, 954, HU, 64, NEQ, 20, 20, 15, 20, 1010, 0, 13, 130, 0, , 0, JPC, 312, 7, , 4
AL, 09, 2021082800, 03, OFCL, 36, 278N, 896W, 120, 937, HU, 34, NEQ, 130, 110, 80, 110, 1010, 0, 12, 145, 0, , 0, JPC, 315, 6, , 5
AL, 09, 2021082800, 03, OFCL, 36, 278N, 896W, 120, 937, HU, 50, NEQ, 60, 60, 30, 50, 1010, 0, 12, 145, 0, , 0, JPC, 315, 6, , 5
AL, 09, 2021082800, 03, OFCL, 36, 278N, 896W, 120, 937, HU, 64, NEQ, 35, 30, 20, 30, 1010, 0, 12, 145, 0, , 0, JPC, 315, 6, , 5
AL, 09, 2021082800, 03, OFCL, 48, 292N, 908W, 110, 948, HU, 34, NEQ, 110, 110, 70, 70, 1010, 0, 15, 135, 0, , 0, JPC, 323, 5, , 6
AL, 09, 2021082800, 03, OFCL, 48, 292N, 908W, 110, 948, HU, 50, NEQ, 60, 60, 40, 40, 1010, 0, 15, 135, 0, , 0, JPC, 323, 5, , 6
AL, 09, 2021082800, 03, OFCL, 48, 292N, 908W, 110, 948, HU, 64, NEQ, 30, 30, 20, 20, 1010, 0, 15, 135, 0, , 0, JPC, 323, 5, , 6
AL, 09, 2021082800, 03, OFCL, 60, 307N, 910W, 60, 992, TS, 34, NEQ, 70, 100, 60, 50, 1010, 0, 21, 75, 0, , 0, JPC, 353, 4, , 7
AL, 09, 2021082800, 03, OFCL, 60, 307N, 910W, 60, 992, TS, 50, NEQ, 30, 30, 20, 20, 1010, 0, 21, 75, 0, , 0, JPC, 353, 4, , 7
AL, 09, 2021082800, 03, OFCL, 72, 323N, 906W, 35, 1004, TS, 34, NEQ, 0, 150, 0, 0, 1010, 0, 73, 45, 0, , 0, JPC, 12, 4, , 8
AL, 09, 2021082800, 03, OFCL, 96, 349N, 882W, 25, 1007, TD, 34, NEQ, 0, 0, 0, 0, 1010, 0, 100, 35, 0, , 0, JPC, 37, 4, , 9
AL, 09, 2021082800, 03, OFCL, 120, 368N, 839W, 20, 1008, LO, 34, NEQ, 0, 0, 0, 0, 1010, 0, 110, 30, 0, , 0, JPC, 60, 5, , 10
7 changes: 7 additions & 0 deletions tests/data/reference/test_no_isotach_adj/original.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
{
"cross_track": NaN,
"along_track": NaN,
"max_sustained_wind_speed": NaN,
"radius_of_maximum_winds": NaN,
"weight": 1.0
}
21 changes: 21 additions & 0 deletions tests/data/reference/test_no_isotach_adj/vortex_4_variable_korobov_1.22
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
AL, 09, 2021082800, 03, OFCL, 0, 225N, 834W, 66, 988, HU, 34, NEQ, 100, 62, 42, 70, 1010, 0, 22, 85, 0, , 0, JPC, 322, 7, , 1
AL, 09, 2021082800, 03, OFCL, 0, 225N, 834W, 66, 988, HU, 50, NEQ, 50, 32, 32, 31, 1010, 0, 22, 85, 0, , 0, JPC, 322, 7, , 1
AL, 09, 2021082800, 03, OFCL, 0, 225N, 834W, 66, 988, HU, 64, NEQ, 22, 0, 0, 0, 1010, 0, 22, 85, 0, , 0, JPC, 322, 7, , 1
AL, 09, 2021082800, 03, OFCL, 3, 230N, 839W, 65, 992, HU, 34, NEQ, 98, 61, 41, 68, 1010, 0, 21, 85, 0, , 0, JPC, 322, 7, , 2
AL, 09, 2021082800, 03, OFCL, 3, 230N, 839W, 65, 992, HU, 50, NEQ, 49, 31, 0, 28, 1010, 0, 21, 85, 0, , 0, JPC, 322, 7, , 2
AL, 09, 2021082800, 03, OFCL, 12, 245N, 853W, 76, 981, HU, 34, NEQ, 111, 83, 62, 89, 1010, 0, 21, 105, 0, , 0, JPC, 319, 7, , 3
AL, 09, 2021082800, 03, OFCL, 12, 245N, 853W, 76, 981, HU, 50, NEQ, 52, 43, 19, 41, 1010, 0, 21, 105, 0, , 0, JPC, 319, 7, , 3
AL, 09, 2021082800, 03, OFCL, 12, 245N, 853W, 76, 981, HU, 64, NEQ, 19, 20, 0, 0, 1010, 0, 21, 105, 0, , 0, JPC, 319, 7, , 3
AL, 09, 2021082800, 03, OFCL, 24, 263N, 872W, 91, 968, HU, 34, NEQ, 130, 111, 79, 108, 1010, 0, 19, 130, 0, , 0, JPC, 316, 6, , 4
AL, 09, 2021082800, 03, OFCL, 24, 263N, 872W, 91, 968, HU, 50, NEQ, 69, 58, 37, 47, 1010, 0, 19, 130, 0, , 0, JPC, 316, 6, , 4
AL, 09, 2021082800, 03, OFCL, 24, 263N, 872W, 91, 968, HU, 64, NEQ, 27, 27, 18, 26, 1010, 0, 19, 130, 0, , 0, JPC, 316, 6, , 4
AL, 09, 2021082800, 03, OFCL, 36, 279N, 889W, 104, 955, HU, 34, NEQ, 147, 127, 95, 126, 1010, 0, 19, 145, 0, , 0, JPC, 318, 6, , 5
AL, 09, 2021082800, 03, OFCL, 36, 279N, 889W, 104, 955, HU, 50, NEQ, 74, 73, 40, 62, 1010, 0, 19, 145, 0, , 0, JPC, 318, 6, , 5
AL, 09, 2021082800, 03, OFCL, 36, 279N, 889W, 104, 955, HU, 64, NEQ, 45, 40, 28, 39, 1010, 0, 19, 145, 0, , 0, JPC, 318, 6, , 5
AL, 09, 2021082800, 03, OFCL, 48, 291N, 898W, 92, 967, HU, 34, NEQ, 118, 118, 82, 82, 1010, 0, 23, 135, 0, , 0, JPC, 327, 4, , 6
AL, 09, 2021082800, 03, OFCL, 48, 291N, 898W, 92, 967, HU, 50, NEQ, 68, 67, 49, 49, 1010, 0, 23, 135, 0, , 0, JPC, 327, 4, , 6
AL, 09, 2021082800, 03, OFCL, 48, 291N, 898W, 92, 967, HU, 64, NEQ, 37, 37, 20, 20, 1010, 0, 23, 135, 0, , 0, JPC, 327, 4, , 6
AL, 09, 2021082800, 03, OFCL, 60, 301N, 899W, 40, 1002, TS, 34, NEQ, 49, 54, 33, 37, 1010, 0, 29, 75, 0, , 0, JPC, 353, 3, , 7
AL, 09, 2021082800, 03, OFCL, 72, 312N, 896W, 15, 1009, TS, 34, NEQ, 0, 0, 0, 0, 1010, 0, 82, 45, 0, , 0, JPC, 16, 3, , 8
AL, 09, 2021082800, 03, OFCL, 96, 329N, 876W, 15, 1008, TD, 34, NEQ, 0, 0, 0, 0, 1010, 0, 109, 35, 0, , 0, JPC, 43, 3, , 9
AL, 09, 2021082800, 03, OFCL, 120, 341N, 838W, 15, 1008, LO, 34, NEQ, 0, 0, 0, 0, 1010, 0, 120, 30, 0, , 0, JPC, 69, 4, , 10
7 changes: 7 additions & 0 deletions tests/data/reference/test_no_isotach_adj/vortex_4_variable_korobov_1.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
{
"cross_track": -1.2815515655446004,
"along_track": 0.5244005127080407,
"max_sustained_wind_speed": 1.2815515655446004,
"radius_of_maximum_winds": -0.5244005127080409,
"weight": 1.0
}
Loading

0 comments on commit d834b8a

Please sign in to comment.