diff --git a/src/2d/shallow/src2.f90 b/src/2d/shallow/src2.f90 index dc603c36e..376daf7e0 100644 --- a/src/2d/shallow/src2.f90 +++ b/src/2d/shallow/src2.f90 @@ -6,7 +6,7 @@ subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) use geoclaw_module, only: manning_break, num_manning use geoclaw_module, only: spherical_distance, coordinate_system use geoclaw_module, only: RAD2DEG, pi, dry_tolerance, DEG2RAD - use geoclaw_module, only: rho_air + use geoclaw_module, only: rho_air, rho use geoclaw_module, only: earth_radius, sphere_source use storm_module, only: wind_forcing, pressure_forcing, wind_drag @@ -40,10 +40,6 @@ subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) ! friction source term real(kind=8), parameter :: depth_tolerance = 1.0d-30 - ! Physics - ! Nominal density of water - real(kind=8), parameter :: rho = 1025.d0 - ! ---------------------------------------------------------------- ! Spherical geometry source term(s) ! @@ -154,7 +150,7 @@ subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) endif wind_speed = sqrt(aux(wind_index,i,j)**2 & + aux(wind_index+1,i,j)**2) - tau = wind_drag(wind_speed, phi) * rho_air * wind_speed / rho + tau = wind_drag(wind_speed, phi) * rho_air * wind_speed / rho(1) q(2,i,j) = q(2,i,j) + dt * tau * aux(wind_index,i,j) q(3,i,j) = q(3,i,j) + dt * tau * aux(wind_index+1,i,j) endif @@ -195,8 +191,8 @@ subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) ! ! Modify momentum in each layer ! if (h > dry_tolerance) then - ! q(2, i, j) = q(2, i, j) - dt * h * P_gradient(1) / rho - ! q(3, i, j) = q(3, i, j) - dt * h * P_gradient(2) / rho + ! q(2, i, j) = q(2, i, j) - dt * h * P_gradient(1) / rho(1) + ! q(3, i, j) = q(3, i, j) - dt * h * P_gradient(2) / rho(1) ! end if ! enddo ! enddo diff --git a/src/2d/shallow/surge/model_storm_module.f90 b/src/2d/shallow/surge/model_storm_module.f90 index cce7445a9..b49e98120 100644 --- a/src/2d/shallow/surge/model_storm_module.f90 +++ b/src/2d/shallow/surge/model_storm_module.f90 @@ -51,6 +51,19 @@ module model_storm_module end type model_storm_type + ! How to deterimine which way a storm should be made to spin + ! The default defers simply to assuming y is a latitude + ! Here either an integer or bool can be returned but as implemented 1 + ! refers to the Northern Hemisphere and therefore causes the storm to spin + ! in a counter-clockwise direction. + abstract interface + logical pure function rotation_def(x, y) + implicit none + real(kind=8), intent(in) :: x, y + end function rotation_def + end interface + procedure(rotation_def), pointer :: rotation + ! Interal tracking variables for storm integer, private :: last_storm_index @@ -68,12 +81,6 @@ module model_storm_module ! track to be close to but not equal the start time of the simulation real(kind=8), parameter :: TRACKING_TOLERANCE = 1d-10 - ! Global constants #TODO: Some of these are in geoclaw already - real(kind=8) :: pi=3.1415927 - real(kind=8) :: omega=7.2921e-5 - real(kind=8) :: chi=1.0 - real(kind=8) :: alpha=1.0 - contains @@ -153,8 +160,8 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit) 0.5d0 * (x(1) + y(1)), y(2)) storm%velocity(2, i) = sign(ds / dt, y(2) - x(2)) else - storm%velocity(1, i) = abs(x(2) - x(1)) / dt - storm%velocity(2, i) = abs(y(2) - y(1)) / dt + storm%velocity(1, i) = (y(1) - x(1)) / dt + storm%velocity(2, i) = (y(2) - x(2)) / dt end if end do @@ -338,7 +345,7 @@ subroutine get_storm_data(t, storm, location, velocity, max_wind_radius, & max_wind_speed, central_pressure, & radius, central_pressure_change) - use geoclaw_module, only: deg2rad, latlon2xy, xy2latlon + use geoclaw_module, only: deg2rad, latlon2xy, xy2latlon, coordinate_system implicit none @@ -382,13 +389,20 @@ subroutine get_storm_data(t, storm, location, velocity, max_wind_radius, & ! Convert coordinates temporarily to meters so that we can use ! the pre-calculated m/s velocities from before - x = latlon2xy(storm%track(2:3,i),storm%track(2:3,i)) - x = x + (t - storm%track(1,i)) * storm%velocity(:,i) - - fn = [xy2latlon(x,storm%track(2:3,i)), & - storm%velocity(:,i), storm%max_wind_radius(i), & - storm%max_wind_speed(i), storm%central_pressure(i), & - storm%radius(i), storm%central_pressure_change(i)] + if (coordinate_system == 2) then + x = latlon2xy(storm%track(2:3,i),storm%track(2:3,i)) + x = x + (t - storm%track(1,i)) * storm%velocity(:,i) + fn = [xy2latlon(x,storm%track(2:3,i)), & + storm%velocity(:,i), storm%max_wind_radius(i), & + storm%max_wind_speed(i), storm%central_pressure(i), & + storm%radius(i), storm%central_pressure_change(i)] + else + x = x + (t - storm%track(1,i)) * storm%velocity(:,i) + fn = [x, & + storm%velocity(:,i), storm%max_wind_radius(i), & + storm%max_wind_speed(i), storm%central_pressure(i), & + storm%radius(i), storm%central_pressure_change(i)] + end if else ! Inbetween two forecast time points (the function storm_index ! ensures that we are not before the first data point, i.e. i > 1) @@ -430,11 +444,8 @@ end subroutine get_storm_data pure subroutine adjust_max_wind(tv, mws, mod_mws, convert_height) real (kind=8), intent(inout) :: tv(2) - real (kind=8), intent(in) :: mws - logical, intent(in) :: convert_height - real (kind=8), intent(out) :: mod_mws real (kind=8) :: trans_speed, trans_mod @@ -521,21 +532,21 @@ end subroutine post_process_wind_estimate ! ========================================================================== pure subroutine calculate_polar_coordinate(x, y, sloc, r, theta) - use geoclaw_module, only: deg2rad, coordinate_system - use geoclaw_module, only: spherical_distance + use geoclaw_module, only: deg2rad, coordinate_system + use geoclaw_module, only: spherical_distance - real(kind=8), intent(in) :: x, y, sloc(2) - real(kind=8), intent(out) :: r, theta + real(kind=8), intent(in) :: x, y, sloc(2) + real(kind=8), intent(out) :: r, theta - if (coordinate_system == 2) then - ! lat-long coordinates, uses Haversine formula - r = spherical_distance(x, y, sloc(1), sloc(2)) - theta = atan2((y - sloc(2)) * DEG2RAD,(x - sloc(1)) * DEG2RAD) - else - ! Strictly cartesian - r = sqrt( (x - sloc(1))**2 + (y - sloc(2))**2) - theta = atan2(y - sloc(2), x - sloc(1)) - end if + if (coordinate_system == 2) then + ! lat-long coordinates, uses Haversine formula + r = spherical_distance(x, y, sloc(1), sloc(2)) + theta = atan2((y - sloc(2)),(x - sloc(1))) + else + ! Strictly cartesian + r = sqrt( (x - sloc(1))**2 + (y - sloc(2))**2) + theta = atan2((y - sloc(2)), (x - sloc(1))) + end if end subroutine calculate_polar_coordinate @@ -653,7 +664,7 @@ subroutine set_holland_1980_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -725,7 +736,7 @@ subroutine set_holland_2008_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -817,7 +828,7 @@ subroutine set_holland_2010_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -974,6 +985,7 @@ real(kind=8) function evaluate_inner_derivative(f,r_m,v_m,r_a,v_a) result(dMa) ! Variables real(kind=8) :: denominator, M_a + real(kind=8), parameter :: alpha=1.0 ! Calculate necessary components of the derivative to make ! calculations easier @@ -992,6 +1004,7 @@ real(kind=8) function evaluate_v_a(f,r_m,v_m,r_a) result(v_a) ! Input real(kind=8), intent(in) :: f, r_m, v_m, r_a + real(kind=8), parameter :: alpha=1.0 v_a = ((2.0*(r_a/r_m)**2)/(2.0-alpha+alpha*(r_a/r_m)**2))**(1.0/(2.0-alpha)) v_a = v_a*(0.5*f*r_m**2 + r_m*v_m) @@ -1061,6 +1074,7 @@ subroutine integrate_m_out(f,r_0,r_a,res,m_out) ! Parameters and Other variables real(kind=8) :: V_m, dr, r_guess, r_p integer :: i + real(kind=8), parameter :: chi=1.0 ! Intialize m_out(res) = 0.5*f*r_0**2 @@ -1104,6 +1118,8 @@ subroutine solve_hurricane_wind_parameters(f, r_m, v_m, res, r_0, r_a) real(kind=8) :: dr, r real(kind=8) :: inner_res, outer_res real(kind=8), dimension(1:res) :: v_out, v_in + real(kind=8), parameter :: chi=1.0 + real(kind=8), parameter :: alpha=1.0 ! Initialize guesses for merge and r_0 r_a = 2.0*r_m @@ -1255,8 +1271,8 @@ subroutine set_SLOSH_fields(maux, mbc, mx, my, xlower, ylower, & trans_speed_x = tv(1) * mwr * r / (mwr**2.d0 + r**2.d0) trans_speed_y = tv(2) * mwr * r / (mwr**2.d0 + r**2.d0) - aux(wind_index,i,j) = wind * merge(-1, 1, y >= 0) * sin(theta) + trans_speed_x - aux(wind_index+1,i,j) = wind * merge(1, -1, y >= 0) * cos(theta) + trans_speed_y + aux(wind_index,i,j) = wind * merge(-1, 1, rotation(x, y)) * sin(theta) + trans_speed_x + aux(wind_index+1,i,j) = wind * merge(1, -1, rotation(x, y)) * cos(theta) + trans_speed_y enddo enddo @@ -1331,7 +1347,7 @@ subroutine set_rankine_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -1419,7 +1435,7 @@ subroutine set_modified_rankine_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -1492,7 +1508,7 @@ subroutine set_deMaria_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo @@ -1584,7 +1600,7 @@ subroutine set_willoughby_fields(maux, mbc, mx, my, xlower, ylower, & call post_process_wind_estimate(maux, mbc, mx, my, i, j, wind, aux, & wind_index, pressure_index, r, radius, tv, mod_mws, theta, & - convert_height, y >= 0) + convert_height, rotation(x, y)) enddo enddo diff --git a/src/2d/shallow/surge/storm_module.f90 b/src/2d/shallow/surge/storm_module.f90 index 7af8b46ae..632cbff90 100644 --- a/src/2d/shallow/surge/storm_module.f90 +++ b/src/2d/shallow/surge/storm_module.f90 @@ -11,7 +11,7 @@ module storm_module - use model_storm_module, only: model_storm_type + use model_storm_module, only: model_storm_type, rotation use data_storm_module, only: data_storm_type implicit none @@ -126,7 +126,7 @@ subroutine set_storm(data_file) ! Locals integer, parameter :: unit = 13 - integer :: i, drag_law + integer :: i, drag_law, rotation_override character(len=200) :: storm_file_path, line if (.not.module_setup) then @@ -153,6 +153,17 @@ subroutine set_storm(data_file) stop "*** ERROR *** Invalid wind drag law." end select read(unit,*) pressure_forcing + read(unit,*) rotation_override + select case(rotation_override) + case(0) + rotation => hemisphere_rotation + case(1) + rotation => N_rotation + case(2) + rotation => S_rotation + case default + stop " *** ERROR *** Roation override invalid." + end select read(unit,*) ! Set some parameters @@ -480,4 +491,26 @@ subroutine output_storm_location(t) end subroutine output_storm_location + ! ========================================================================== + ! Default to assuming y is a latitude and if y >= 0 we are want to spin + ! counter-clockwise + ! ========================================================================== + logical pure function hemisphere_rotation(x, y) result(rotation) + implicit none + real(kind=8), intent(in) :: x, y + rotation = (y >= 0.d0) + end function hemisphere_rotation + ! This version just returns the user defined direction + logical pure function N_rotation(x, y) result(rotation) + implicit none + real(kind=8), intent(in) :: x, y + rotation = .true. + end function N_rotation + ! This version just returns the user defined direction + logical pure function S_rotation(x, y) result(rotation) + implicit none + real(kind=8), intent(in) :: x, y + rotation = .false. + end function S_rotation + end module storm_module diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index e3d514ffc..83d9fa769 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -545,17 +545,19 @@ class SurgeData(clawpack.clawutil.data.ClawData): 'SLOSH': 4, 'rankine': 5, 'modified-rankine': 6, - 'DeMaria': 7 + 'DeMaria': 7, + 'willoughby': 9, } - storm_spec_not_implemented = ['CLE'] + storm_spec_not_implemented = ['CLE', 'willoughby'] def __init__(self): - super(SurgeData,self).__init__() + super(SurgeData, self).__init__() # Source term controls - self.add_attribute('wind_forcing',False) - self.add_attribute('drag_law',1) - self.add_attribute('pressure_forcing',False) + self.add_attribute('wind_forcing', False) + self.add_attribute('drag_law', 1) + self.add_attribute('pressure_forcing', False) + self.add_attribute('rotation_override', 0) # Algorithm parameters - Indexing is python based self.add_attribute("wind_index", 4) @@ -563,8 +565,8 @@ def __init__(self): self.add_attribute("display_landfall_time", False) # AMR parameters - self.add_attribute('wind_refine',[20.0,40.0,60.0]) - self.add_attribute('R_refine',[60.0e3,40e3,20e3]) + self.add_attribute('wind_refine', [20.0,40.0,60.0]) + self.add_attribute('R_refine', [60.0e3,40e3,20e3]) # Storm parameters self.add_attribute('storm_type', None) # Backwards compatibility @@ -572,7 +574,7 @@ def __init__(self): self.add_attribute("storm_file", None) # File(s) containing data - def write(self,out_file='surge.data',data_source="setrun.py"): + def write(self, out_file='surge.data', data_source="setrun.py"): """Write out the data file to the path given""" # print "Creating data file %s" % out_file @@ -582,6 +584,19 @@ def write(self,out_file='surge.data',data_source="setrun.py"): self.data_write('drag_law', description='(Type of drag law to use)') self.data_write('pressure_forcing', description="(Pressure source term used)") + if isinstance(self.rotation_override, str): + if self.rotation_override.lower() == "normal": + self.rotation_override = 0 + elif "n" in self.rotation_override.lower(): + self.rotation_override = 1 + elif "s" in self.rotation_override.lower(): + self.rotation_override = 2 + else: + raise ValueError("Unknown rotation_override specification.") + else: + self.rotation_override = int(self.rotation_override) + self.data_write('rotation_override', + description="(Override storm rotation)") self.data_write() self.data_write("wind_index", value=self.wind_index + 1, diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py index 372ca1073..c2ada663f 100644 --- a/src/python/geoclaw/surge/plot.py +++ b/src/python/geoclaw/surge/plot.py @@ -41,6 +41,7 @@ surge_data = geodata.SurgeData() + class track_data(object): """Read in storm track data from run output""" @@ -83,14 +84,16 @@ def gauge_locations(current_data, gaugenos='all'): add_labels=True, xoffset=0.02, yoffset=0.02) + def gauge_dry_regions(cd, dry_tolerance=1e-16): """Masked array of zeros where gauge is dry.""" - return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) > dry_tolerance, - np.zeros(cd.gaugesoln.q[0,:].shape)) + return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) > dry_tolerance, + np.zeros(cd.gaugesoln.q[0, :].shape)) + def gauge_surface(cd, dry_tolerance=1e-16): """Sea surface at gauge masked when dry.""" - return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) < dry_tolerance, + return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) < dry_tolerance, cd.gaugesoln.q[3, :]) @@ -497,19 +500,21 @@ def plot_track(t, x, y, wind_radius, wind_speed, Pc, name=None): # Returns axes # Storm with category plotting function # ======================================================================== -def add_track(Storm, axes, plot_package=None, category_color=None, legend_loc='best', + + +def add_track(Storm, axes, plot_package=None, category_color=None, legend_loc='best', intensity=False, categorization="NHC", limits=None, track_color='red'): if category_color is None: - category_color = {5: 'red', - 4: 'orange', - 3: 'yellow', - 2: 'blue', # edit color - 1: 'violet', - 0: 'black', - -1: 'gray'} + category_color = {5: 'red', + 4: 'orange', + 3: 'yellow', + 2: 'blue', # edit color + 1: 'violet', + 0: 'black', + -1: 'gray'} category = Storm.category(categorization=categorization) - + # make it if intensity = true # basic plotting @@ -525,52 +530,55 @@ def add_track(Storm, axes, plot_package=None, category_color=None, legend_loc='b axes.set_xlabel("Longitude") axes.set_ylabel("Latitude") - - categories_legend = [] - if intensity and categorization is "NHC": + if intensity and categorization == "NHC": categories_legend = [] - # plotitem = plotaxes.new_plotitem(name='category', plot_type='1d_plot') + # plotitem = plotaxes.new_plotitem(name='category', plot_type='1d_plot') if (-1 in category): - negativeone = mlines.Line2D([], [], color=category_color[-1], marker='s', ls='', label="Tropical Depression") + negativeone = mlines.Line2D( + [], [], color=category_color[-1], marker='s', ls='', label="Tropical Depression") categories_legend.append(negativeone) if (0 in category): - zero = mlines.Line2D([], [], color=category_color[0], marker='s', ls='', label="Tropical Storn") + zero = mlines.Line2D( + [], [], color=category_color[0], marker='s', ls='', label="Tropical Storn") categories_legend.append(zero) if (1 in category): - one = mlines.Line2D([], [], color=category_color[1], marker='s', ls='', label="Category 1") + one = mlines.Line2D([], [], color=category_color[1], + marker='s', ls='', label="Category 1") categories_legend.append(one) if (2 in category): - two = mlines.Line2D([], [], color=category_color[2], marker='s', ls='', label="Category 2") + two = mlines.Line2D([], [], color=category_color[2], + marker='s', ls='', label="Category 2") categories_legend.append(two) if (3 in category): - three = mlines.Line2D([], [], color=category_color[3], marker='s', ls='', label="Category 3") + three = mlines.Line2D( + [], [], color=category_color[3], marker='s', ls='', label="Category 3") categories_legend.append(three) if (4 in category): - four = mlines.Line2D([], [], color=category_color[4], marker='s', ls='', label="Category 4") + four = mlines.Line2D( + [], [], color=category_color[4], marker='s', ls='', label="Category 4") categories_legend.append(four) if (5 in category): - five = mlines.Line2D([], [], color=category_color[5], marker='s', ls='', label="Category 5") + five = mlines.Line2D( + [], [], color=category_color[5], marker='s', ls='', label="Category 5") categories_legend.append(five) plt.legend(handles=categories_legend, loc=legend_loc) - + # if bounds is not None: # plotitem.pcolor_cmin = bounds[0] # plotitem.pcolor_cmax = bounds[1] return axes - - # if plot_type == 'pcolor' or plot_type == 'imshow': # plotitem = plotaxes.new_plotitem(name='wind', plot_type='2d_pcolor') # plotitem.plot_var = wind_speed diff --git a/src/python/geoclaw/surge/quality.py b/src/python/geoclaw/surge/quality.py index 6584efbdb..e62db2c98 100644 --- a/src/python/geoclaw/surge/quality.py +++ b/src/python/geoclaw/surge/quality.py @@ -9,13 +9,14 @@ class InstabilityError(Exception): pass + def quality_check(model_output_dir, - regions_to_check = [[-81,22,-40,55], # atlantic - [-100,15,-82,32], # gulf - [-88.5,13.25,-70,19.75]], # carribean - frames_to_check = 4, - mean_tol_abs = .01, - min_depth = 300): + regions_to_check=[[-81, 22, -40, 55], # atlantic + [-100, 15, -82, 32], # gulf + [-88.5, 13.25, -70, 19.75]], # carribean + frames_to_check=4, + mean_tol_abs=.01, + min_depth=300): """Run a simple check to flag runs that were potentially numerically unstable. This function looks through the final *frames_to_check* frames of a model output and checks all cells in AMR level 1 that are below *min_depth* and @@ -23,7 +24,7 @@ def quality_check(model_output_dir, to encompass individual basins). If the absolute value of the mean surface height for these cells, which should not really have much surge, is above/mean_tol_abs, it raises an error. - + :Input: - *model_output_dir* (str) Path to the output directory for a given model - *regions_to_check* (list of lists of floats) A list of 4-element lists @@ -33,17 +34,17 @@ def quality_check(model_output_dir, height within any region in *regions_to_check* is above this value (meters). - *min_depth* (float) Only cells that are below *min_depth* (meters) are checked, in order to avoid including cells that potentially should have large surge values. - + :Raises: - InstabilityError if the model is flagged as potentially unstable in one of the regions. """ - + # find which frame is last - output_files = glob(join(model_output_dir,'fort.b*')) + output_files = glob(join(model_output_dir, 'fort.b*')) last_frame = int(output_files[-1].split('/')[-1][-4:]) # get sl_init - with open(join(model_output_dir,'geoclaw.data'),'r') as f: + with open(join(model_output_dir, 'geoclaw.data'), 'r') as f: for l in f: if '=: sea_level' in l: sl_init = float(l.strip().split()[0]) @@ -51,44 +52,45 @@ def quality_check(model_output_dir, # bring in solution soln = Solution() for frame in range(last_frame-frames_to_check+1, last_frame+1): - soln.read(frame, path = model_output_dir, file_format='binary', read_aux=False) + soln.read(frame, path=model_output_dir, + file_format='binary', read_aux=False) for r in regions_to_check: all_vals = np.array([]) for s in soln.states: # only looking at lowest AMR level - if s.patch.level >1: + if s.patch.level > 1: continue x = s.grid.c_centers[0] y = s.grid.c_centers[1] - eta = s.q[3,:,:] - topo = eta - s.q[0,:,:] - - # only count - eta = np.where(topo<(-min_depth),eta,np.nan) - in_reg = eta[np.ix_((x[:,0]>=r[0]) & (x[:,0]<=r[2]), - (y[0,:]>=r[1]) & (y[0,:]<=r[3]))] - all_vals = np.concatenate([all_vals,in_reg.flatten()]) + eta = s.q[3, :, :] + topo = eta - s.q[0, :, :] + + # only count + eta = np.where(topo < (-min_depth), eta, np.nan) + in_reg = eta[np.ix_((x[:, 0] >= r[0]) & (x[:, 0] <= r[2]), + (y[0, :] >= r[1]) & (y[0, :] <= r[3]))] + all_vals = np.concatenate([all_vals, in_reg.flatten()]) # adjust for sl_init all_vals = all_vals - sl_init - if all_vals.shape[0]>0: + if all_vals.shape[0] > 0: if abs(np.nanmean(all_vals)) > mean_tol_abs: raise InstabilityError("Model possibly unstable due to large magnitude deep " "ocean surge at end of run ({:.1f} cm in region {})".format( - np.nanmean(all_vals)*100, r)) + np.nanmean(all_vals)*100, r)) - -def get_max_boundary_fluxes(model_output_dir, max_refinement_depth_to_check = None, - frames_to_check = 1): + +def get_max_boundary_fluxes(model_output_dir, max_refinement_depth_to_check=None, + frames_to_check=1): """A common cause of instability is a persistant normal flux at the boundary. This code checks the last frame(s) of a model output and returns the maximum magnitude normal fluxes and currents at each boundary, as well as the maximum magnitude fluxes and currents observed within the entire domain. - + :Input: - *model_output_dir* (str) Path to the output directory for a given model - *max_refinement_depth_to_check* (int or None, optional) How many refinement levels to @@ -100,20 +102,20 @@ def get_max_boundary_fluxes(model_output_dir, max_refinement_depth_to_check = No """ # find which frame is last - output_files = glob(join(model_output_dir,'fort.q*')) + output_files = glob(join(model_output_dir, 'fort.q*')) frames = [int(i.split('/')[-1][-4:]) for i in output_files] last_frame = max(frames) - + # get domain and file format - with open(join(model_output_dir,'claw.data'),'r') as f: + with open(join(model_output_dir, 'claw.data'), 'r') as f: for l in f: if '=: lower' in l: - xl,yl = [float(i) for i in l.strip().split()[:2]] + xl, yl = [float(i) for i in l.strip().split()[:2]] elif '=: upper' in l: - xu,yu = [float(i) for i in l.strip().split()[:2]] + xu, yu = [float(i) for i in l.strip().split()[:2]] elif '=: output_format' in l: of = int(l.strip().split()[0]) - 1 - file_format = ['ascii','netcdf','binary'][of] + file_format = ['ascii', 'netcdf', 'binary'][of] maxhxl = -np.inf maxhyl = -np.inf @@ -128,10 +130,11 @@ def get_max_boundary_fluxes(model_output_dir, max_refinement_depth_to_check = No maxcurrx_overall = -np.inf maxcurry_overall = -np.inf - for f in range(last_frame+1-frames_to_check,last_frame+1): + for f in range(last_frame+1-frames_to_check, last_frame+1): soln = Solution() - soln.read(f, path = model_output_dir, file_format=file_format, read_aux=False) + soln.read(f, path=model_output_dir, + file_format=file_format, read_aux=False) for s in soln.states: @@ -139,50 +142,56 @@ def get_max_boundary_fluxes(model_output_dir, max_refinement_depth_to_check = No if max_refinement_depth_to_check is not None: if s.patch.level > max_refinement_depth_to_check: continue - + # get rounding error tolerance delta = s.grid.dimensions[0].delta edge_tol = delta * .001 - + x = s.grid.c_centers[0] xedges = s.grid.c_nodes[0] y = s.grid.c_centers[1] yedges = s.grid.c_nodes[1] - eta = s.q[3,:,:] - h = s.q[0,:,:] - hx = s.q[1,:,:] + eta = s.q[3, :, :] + h = s.q[0, :, :] + hx = s.q[1, :, :] curr_x = hx / h - hy = s.q[2,:,:] + hy = s.q[2, :, :] curr_y = hy / h - topo = eta - s.q[0,:,:] - maxhx_overall = np.nanmax([np.abs(hx).max(),maxhx_overall]) - maxhy_overall = np.nanmax([np.abs(hy).max(),maxhy_overall]) - maxcurrx_overall = np.nanmax([np.nanmax(np.abs(curr_x)),maxcurrx_overall]) - maxcurry_overall = np.nanmax([np.nanmax(np.abs(curr_y)),maxcurry_overall]) - - if abs(xedges[0,0] - xl) < edge_tol: - maxhxl = np.nanmax([maxhxl,np.nanmax(np.abs(hx[0,:]))]) - maxcurrxl = np.nanmax([maxcurrxl,np.nanmax(np.abs(curr_x[0,:]))]) - if abs(xedges[-1,0] - xu) < edge_tol: - maxhxu = np.nanmax([maxhxu,np.nanmax(np.abs(hx[-1,:]))]) - maxcurrxu = np.nanmax([maxcurrxu,np.nanmax(np.abs(curr_x[-1,:]))]) - if abs(yedges[0,0] - yl) < edge_tol: - maxhyl = np.nanmax([maxhyl,np.nanmax(np.abs(hy[:,0]))]) - maxcurryl = np.nanmax([maxcurryl,np.nanmax(np.abs(curr_y[:,0]))]) - if abs(yedges[0,-1] - yu) < edge_tol: - maxhyu = np.nanmax([maxhyu,np.nanmax(np.abs(hy[:,-1]))]) - maxcurryu = np.nanmax([maxcurryu,np.nanmax(np.abs(curr_y[:,-1]))]) - - return ({'max_normal_fluxes':{'W': maxhxl, - 'E': maxhxu, - 'N': maxhyu, - 'S': maxhyl}, - 'max_normal_currents':{'W': maxcurrxl, - 'E': maxcurrxu, - 'N': maxcurryu, - 'S': maxcurryl}, + topo = eta - s.q[0, :, :] + maxhx_overall = np.nanmax([np.abs(hx).max(), maxhx_overall]) + maxhy_overall = np.nanmax([np.abs(hy).max(), maxhy_overall]) + maxcurrx_overall = np.nanmax( + [np.nanmax(np.abs(curr_x)), maxcurrx_overall]) + maxcurry_overall = np.nanmax( + [np.nanmax(np.abs(curr_y)), maxcurry_overall]) + + if abs(xedges[0, 0] - xl) < edge_tol: + maxhxl = np.nanmax([maxhxl, np.nanmax(np.abs(hx[0, :]))]) + maxcurrxl = np.nanmax( + [maxcurrxl, np.nanmax(np.abs(curr_x[0, :]))]) + if abs(xedges[-1, 0] - xu) < edge_tol: + maxhxu = np.nanmax([maxhxu, np.nanmax(np.abs(hx[-1, :]))]) + maxcurrxu = np.nanmax( + [maxcurrxu, np.nanmax(np.abs(curr_x[-1, :]))]) + if abs(yedges[0, 0] - yl) < edge_tol: + maxhyl = np.nanmax([maxhyl, np.nanmax(np.abs(hy[:, 0]))]) + maxcurryl = np.nanmax( + [maxcurryl, np.nanmax(np.abs(curr_y[:, 0]))]) + if abs(yedges[0, -1] - yu) < edge_tol: + maxhyu = np.nanmax([maxhyu, np.nanmax(np.abs(hy[:, -1]))]) + maxcurryu = np.nanmax( + [maxcurryu, np.nanmax(np.abs(curr_y[:, -1]))]) + + return ({'max_normal_fluxes': {'W': maxhxl, + 'E': maxhxu, + 'N': maxhyu, + 'S': maxhyl}, + 'max_normal_currents': {'W': maxcurrxl, + 'E': maxcurrxu, + 'N': maxcurryu, + 'S': maxcurryl}, 'domain_maxs': {'hu': maxhx_overall, 'hv': maxhy_overall, 'u': maxcurrx_overall, - 'v': maxcurry_overall}}) \ No newline at end of file + 'v': maxcurry_overall}}) diff --git a/src/python/geoclaw/surge/storm.py b/src/python/geoclaw/surge/storm.py index 6e48eacc8..dbcfad3ff 100644 --- a/src/python/geoclaw/surge/storm.py +++ b/src/python/geoclaw/surge/storm.py @@ -37,10 +37,10 @@ import argparse import datetime -import numpy +import numpy as np +import pandas as pd import clawpack.geoclaw.units as units -import clawpack.clawutil.data # ============================================================================= # Common acronyms across formats @@ -61,12 +61,12 @@ TCVitals_Basins = {"L": "North Atlantic", "E": "North East Pacific", "C": "North Central Pacific", - "W": "North West Pacific", - "B": "Bay of Bengal (North Indian Ocean)", - "A": "Arabian Sea (North Indian Ocean)", - "Q": "South Atlantic", - "P": "South Pacific", - "S": "South Indian Ocean"} + "W": "North West Pacific", + "B": "Bay of Bengal (North Indian Ocean)", + "A": "Arabian Sea (North Indian Ocean)", + "Q": "South Atlantic", + "P": "South Pacific", + "S": "South Indian Ocean"} # Tropical Cyclone Designations # see https://www.nrlmry.navy.mil/atcf_web/docs/database/new/abrdeck.html @@ -113,6 +113,7 @@ missing_necessary_data_warning_str = """No storm points in the input file had both a max wind speed and a central pressure observation.""" + class NoDataError(ValueError): """Exception to raise when no valid data in input file""" pass @@ -136,25 +137,26 @@ class Storm(object): *TODO:* Add description of unit handling :Attributes: - - *t* (list(datetime.datetiem)) Contains the time at which each entry of - the other arrays are at. These are expected to be *datetime* objects. - Note that when written some formats require a *time_offset* to be set. - - *eye_location* (ndarray(:, :)) location of the eye of the storm. - Default units are in signed decimcal longitude and latitude. + - *t* (list(float) or list(datetime.datetiem)) Contains the time at which + each entry of the other arrays are at. These are expected to + be *datetime* objects. Note that when written some formats require + a *time_offset* to be set. + - *eye_location* (ndarray(:, :)) location of the eye of the storm. Default + units are in signed decimal longitude and latitude. - *max_wind_speed* (ndarray(:)) Maximum wind speed. Default units are meters/second. - *max_wind_radius* (ndarray(:)) Radius at which the maximum wind speed occurs. Default units are meters. - - *central_pressure* (ndarray(:)) Central pressure of storm. Default - units are Pascals. + - *central_pressure* (ndarray(:)) Central pressure of storm. Default units + are Pascals. - *storm_radius* (ndarray(:)) Radius of storm, often defined as the last closed iso-bar of pressure. Default units are meters. - *time_offset* (datetime.datetime) A date time that as an offset for the - simulation time. This will default to the beginning of the first of the - year that the first time point is found in. + simulation time. This will default to the beginning of the first of + the year that the first time point is found in. - *wind_speeds* (ndarray(:, :)) Wind speeds defined in every record, such - as 34kt, 50kt, 64kt, etc and their radii. Default units are meters/second - and meters. + as 34kt, 50kt, 64kt, etc and their radii. Default units are + meters/second and meters. :Initialization: 1. Read in existing file at *path*. @@ -194,7 +196,7 @@ def __init__(self, path=None, file_format="ATCF", **kwargs): self.wind_speeds = None # Storm descriptions - not all formats provide these - self.name = None + self.name = None # Possibly a list of a storm's names self.basin = None # Basin containing storm self.ID = None # ID code - depends on format self.classification = None # Classification of storm (e.g. HU) @@ -207,10 +209,12 @@ def __init__(self, path=None, file_format="ATCF", **kwargs): # Basic object support def __str__(self): r"""""" - output = "Name: %s" % self.name - output = "\n".join((output, "Dates: %s - %s" % (self.t[0].isoformat(), - self.t[-1].isoformat()) - )) + output = f"Name: {self.name}\n" + if isinstance(self.t[0], datetime.datetime): + output += f"Dates: {self.t[0].isoformat()}" + output += f" - {self.t[-1].isoformat()}" + else: + output += f"Dates: {self.t[0]} - {self.t[-1]}" return output def __repr__(self): @@ -270,15 +274,20 @@ def read_geoclaw(self, path, verbose=False): - *verbose* (bool) Output more info regarding reading. """ + # Attempt to get name from file if is follows the convention name.storm + if ".storm" in os.path.splitext(path): + self.name = os.path.splitext(os.path.basename(path))[0] + + # Read header with open(path, 'r') as data_file: num_casts = int(data_file.readline()) self.time_offset = datetime.datetime.strptime( - data_file.readline()[:19], - '%Y-%m-%dT%H:%M:%S') - - data = numpy.loadtxt(path, skiprows=3) + data_file.readline()[:19], + '%Y-%m-%dT%H:%M:%S') + # Read rest of data + data = np.loadtxt(path, skiprows=3) num_forecasts = data.shape[0] - self.eye_location = numpy.empty((2, num_forecasts)) + self.eye_location = np.empty((2, num_forecasts)) assert(num_casts == num_forecasts) self.t = [self.time_offset + datetime.timedelta(seconds=data[i, 0]) for i in range(num_forecasts)] @@ -299,70 +308,85 @@ def read_atcf(self, path, verbose=False): - *path* (string) Path to the file to be read. - *verbose* (bool) Output more info regarding reading. """ - try: - import pandas as pd - except ImportError as e: - print("read_atcf currently requires pandas to work.") - raise e # See here for the ATCF format documentation: # https://www.nrlmry.navy.mil/atcf_web/docs/database/new/abdeck.txt + + # Slightly more robust converter for ATCF data fields that can be + # missing + def num_converter(x): + if isinstance(x, str): + if len(x.strip()) == 0: + # Only whitespace + return np.nan + else: + # Assume this is still a number + return float(x) + elif x is None: + return np.nan + return float(x) + df = pd.read_csv(path, engine="python", sep=",+", names=[ - "BASIN", "CY", "YYYYMMDDHH", "TECHNUM", "TECH", "TAU", - "LAT", "LON", "VMAX", "MSLP", "TY", - "RAD", "WINDCODE", "RAD1", "RAD2", "RAD3", "RAD4", - "POUTER", "ROUTER", "RMW", "GUSTS", "EYE", "SUBREGION", - "MAXSEAS", "INITIALS", "DIR", "SPEED", "STORMNAME", "DEPTH", - "SEAS", "SEASCODE", "SEAS1", "SEAS2", "SEAS3", "SEAS4", - "USERDEFINE1", "userdata1", - "USERDEFINE2", "userdata2", - "USERDEFINE3", "userdata3", - "USERDEFINE4", "userdata4", - "USERDEFINE5", "userdata5", - ], + "BASIN", "CY", "YYYYMMDDHH", "TECHNUM", "TECH", "TAU", + "LAT", "LON", "VMAX", "MSLP", "TY", + "RAD", "WINDCODE", "RAD1", "RAD2", "RAD3", "RAD4", + "POUTER", "ROUTER", "RMW", "GUSTS", "EYE", "SUBREGION", + "MAXSEAS", "INITIALS", "DIR", "SPEED", "STORMNAME", "DEPTH", + "SEAS", "SEASCODE", "SEAS1", "SEAS2", "SEAS3", "SEAS4", + "USERDEFINE1", "userdata1", + "USERDEFINE2", "userdata2", + "USERDEFINE3", "userdata3", + "USERDEFINE4", "userdata4", + "USERDEFINE5", "userdata5", + ], converters={ "YYYYMMDDHH": lambda d: datetime.datetime( int(d[1:5]), int(d[5:7]), int(d[7:9]), int(d[9:11])), "TAU": lambda d: datetime.timedelta(hours=int(d)), "LAT": lambda d: (-.1 if d[-1] == "S" else .1) * int(d.strip("NS ")), "LON": lambda d: (-.1 if d[-1] == "W" else .1) * int(d.strip("WE ")), - }, + "RAD": num_converter, + "RAD": num_converter, + "RAD1": num_converter, + "RAD2": num_converter, + "RAD3": num_converter, + "RAD4": num_converter, + "ROUTER": num_converter, + "RMW": num_converter, + "STORMNAME": lambda d: (d.strip() if isinstance(d, str) else d) + }, dtype={ "BASIN": str, "CY": int, "VMAX": float, "MSLP": float, - "TY": str, - "RAD": float, - "RAD1": float, - "RAD2": float, - "RAD3": float, - "RAD4": float, - "ROUTER": float, - "RMW": float, - }) + "TY": str + }) # Grab data regarding basin and cyclone number from first row self.basin = ATCF_basins[df["BASIN"][0]] self.ID = df["CY"][0] + # Keep around the name as an array + self.name = df["STORMNAME"].to_numpy() + # Take forecast period TAU into consideration df['DATE'] = df["YYYYMMDDHH"] + df["TAU"] df = df[["DATE", "TAU", "TY", "LAT", "LON", "VMAX", "MSLP", - "ROUTER", "RMW", "RAD", "RAD1", "RAD2", "RAD3", "RAD4",]] + "ROUTER", "RMW", "RAD", "RAD1", "RAD2", "RAD3", "RAD4", ]] df = df.sort_values(by=["DATE", "TAU"]).reset_index(drop=True) - # For each DATE, choose best (smallest TAU) available data for c in ["LAT", "LON", "VMAX", "MSLP", "ROUTER", "RMW", "RAD", "RAD1", "RAD2", "RAD3", "RAD4"]: - df[c] = df[c].where(df[c] != 0, numpy.nan) # value 0 means NaN + df[c] = df[c].where(df[c] != 0, np.nan) # value 0 means NaN df[c] = df.groupby("DATE")[c].bfill() df = df.groupby("DATE").first() # Wind profile (occasionally missing for older ATCF storms) # Wind speeds and their radii - df["RAD_MEAN"] = df[["RAD1", "RAD2", "RAD3", "RAD4"]].mean(axis=1, skipna=True) + df["RAD_MEAN"] = df[["RAD1", "RAD2", "RAD3", "RAD4"]].mean( + axis=1, skipna=True) df = df.drop(["TAU", "RAD1", "RAD2", "RAD3", "RAD4"], axis=1) df = df.dropna(how="any", subset=["LAT", "LON"]) @@ -380,25 +404,17 @@ def read_atcf(self, path, verbose=False): # max_wind_radius - convert from nm to m - 1.8520000031807990 * 1000.0 # central_pressure - convert from mbar to Pa - 100.0 # Radius of last isobar contour - convert from nm to m - 1.852000003180799d0 * 1000.0 - self.max_wind_speed = units.convert(df["VMAX"].to_numpy(), 'knots', 'm/s') - self.central_pressure = units.convert(df["MSLP"].to_numpy(), 'mbar', 'Pa') + self.max_wind_speed = units.convert( + df["VMAX"].to_numpy(), 'knots', 'm/s') + self.central_pressure = units.convert( + df["MSLP"].to_numpy(), 'mbar', 'Pa') self.max_wind_radius = units.convert(df["RMW"].to_numpy(), 'nmi', 'm') self.storm_radius = units.convert(df["ROUTER"].to_numpy(), 'nmi', 'm') - self.wind_speeds = df[["RAD","RAD_MEAN"]].to_numpy() - self.wind_speeds[:,0] = units.convert(self.wind_speeds[:,0], 'knots', 'm/s') - self.wind_speeds[:,1] = units.convert(self.wind_speeds[:,1], 'nmi', 'm') - - # Set NaNs to -1 to mark them as missing - for ar in [self.max_wind_speed, self.central_pressure, - self.max_wind_radius, self.storm_radius, self.wind_speeds]: - ar[numpy.isnan(ar)] = -1. - - if self.max_wind_speed.min() == -1: - warnings.warn('Some timesteps have missing max wind speed. These will not be written' - ' out to geoclaw format.') - if self.central_pressure.min() == -1: - warnings.warn('Some timesteps have missing central pressure. These will not be written' - ' out to geoclaw format.') + self.wind_speeds = df[["RAD", "RAD_MEAN"]].to_numpy() + self.wind_speeds[:, 0] = units.convert( + self.wind_speeds[:, 0], 'knots', 'm/s') + self.wind_speeds[:, 1] = units.convert( + self.wind_speeds[:, 1], 'nmi', 'm') def read_hurdat(self, path, verbose=False): r"""Read in HURDAT formatted storm file @@ -437,13 +453,13 @@ def read_hurdat(self, path, verbose=False): # Parse data block self.t = [] - self.event = numpy.empty(num_lines, dtype=str) - self.classification = numpy.empty(num_lines, dtype=str) - self.eye_location = numpy.empty((num_lines, 2)) - self.max_wind_speed = numpy.empty(num_lines) - self.central_pressure = numpy.empty(num_lines) - self.max_wind_radius = numpy.empty(num_lines) - self.storm_radius = numpy.empty(num_lines) + self.event = np.empty(num_lines, dtype=str) + self.classification = np.empty(num_lines, dtype=str) + self.eye_location = np.empty((num_lines, 2)) + self.max_wind_speed = np.empty(num_lines) + self.central_pressure = np.empty(num_lines) + self.max_wind_radius = np.empty(num_lines) + self.storm_radius = np.empty(num_lines) for (i, line) in enumerate(data_block): if len(line) == 0: @@ -480,28 +496,30 @@ def read_hurdat(self, path, verbose=False): # Intensity information - radii are not included directly in this # format and instead radii of winds above a threshold are included - self.max_wind_speed[i] = units.convert(float(data[6]), 'knots', 'm/s') - self.central_pressure[i] = units.convert(float(data[7]), 'mbar', 'Pa') + self.max_wind_speed[i] = units.convert( + float(data[6]), 'knots', 'm/s') + self.central_pressure[i] = units.convert( + float(data[7]), 'mbar', 'Pa') warnings.warn(missing_data_warning_str) self.max_wind_radius[i] = -1 self.storm_radius[i] = -1 def read_ibtracs(self, path, sid=None, storm_name=None, year=None, start_date=None, - agency_pref = ['wmo', - 'usa', - 'tokyo', - 'newdelhi', - 'reunion', - 'bom', - 'nadi', - 'wellington', - 'cma', - 'hko', - 'ds824', - 'td9636', - 'td9635', - 'neumann', - 'mlc']): + agency_pref=['wmo', + 'usa', + 'tokyo', + 'newdelhi', + 'reunion', + 'bom', + 'nadi', + 'wellington', + 'cma', + 'hko', + 'ds824', + 'td9636', + 'td9635', + 'neumann', + 'mlc']): r"""Read in IBTrACS formatted storm file This reads in the netcdf-formatted IBTrACS v4 data. You must either pass @@ -543,7 +561,8 @@ def read_ibtracs(self, path, sid=None, storm_name=None, year=None, start_date=No # only allow one method for specifying storms if (sid is not None) and ((storm_name is not None) or (year is not None)): - raise ValueError('Cannot specify both *sid* and *storm_name* or *year*.') + raise ValueError( + 'Cannot specify both *sid* and *storm_name* or *year*.') with xr.open_dataset(path) as ds: @@ -554,7 +573,7 @@ def read_ibtracs(self, path, sid=None, storm_name=None, year=None, start_date=No else: storm_name = storm_name.upper() # in case storm is unnamed - if storm_name.upper() in ['UNNAMED','NO-NAME']: + if storm_name.upper() in ['UNNAMED', 'NO-NAME']: storm_name = 'NOT_NAMED' storm_match = (ds.name == storm_name.encode()) year_match = (ds.time.dt.year == year).any(dim='date_time') @@ -567,10 +586,11 @@ def read_ibtracs(self, path, sid=None, storm_name=None, year=None, start_date=No raise ValueError('Storm/year not found in provided file') else: # see if a date was provided for multiple unnamed storms - assert start_date is not None, ValueError('Multiple storms identified and no start_date specified.') + assert start_date is not None, ValueError( + 'Multiple storms identified and no start_date specified.') start_times = ds.time.isel(date_time=0) - start_date = numpy.datetime64(start_date) + start_date = np.datetime64(start_date) # find storm with start date closest to provided storm_ix = abs(start_times - start_date).argmin() @@ -583,36 +603,35 @@ def read_ibtracs(self, path, sid=None, storm_name=None, year=None, start_date=No raise ValueError('No valid wind speeds found for this storm.') ds = ds.sel(date_time=valid_t) - # list of the agencies that correspond to 'usa_*' variables usa_agencies = [b'atcf', b'hurdat_atl', b'hurdat_epa', b'jtwc_ep', - b'nhc_working_bt', b'tcvightals', b'tcvitals'] + b'nhc_working_bt', b'tcvightals', b'tcvitals'] - - ## Create mapping from wmo_ or usa_agency - ## to the appropriate variable - agency_map = {b'':agency_pref.index('wmo')} + # Create mapping from wmo_ or usa_agency + # to the appropriate variable + agency_map = {b'': agency_pref.index('wmo')} # account for multiple usa agencies for a in usa_agencies: agency_map[a] = agency_pref.index('usa') # map all other agencies to themselves - for i in [a for a in agency_pref if a not in ['wmo','usa']]: + for i in [a for a in agency_pref if a not in ['wmo', 'usa']]: agency_map[i.encode('utf-8')] = agency_pref.index(i) # fill in usa as provider if usa_agency is # non-null when wmo_agency is null - provider = ds.wmo_agency.where(ds.wmo_agency!=b'',ds.usa_agency) + provider = ds.wmo_agency.where(ds.wmo_agency != b'', ds.usa_agency) # get index into from agency that is wmo_provider def map_val_to_ix(a): - func = lambda x: agency_map[x] - return xr.apply_ufunc(func,a, vectorize=True) - pref_agency_ix=map_val_to_ix(provider) + def func(x): return agency_map[x] + return xr.apply_ufunc(func, a, vectorize=True) + pref_agency_ix = map_val_to_ix(provider) - ## GET MAX WIND SPEED and PRES + # GET MAX WIND SPEED and PRES pref_vals = {} - for v in ['wind','pres']: - all_vals = ds[['{}_{}'.format(i,v) for i in agency_pref]].to_array(dim='agency') + for v in ['wind', 'pres']: + all_vals = ds[['{}_{}'.format(i, v) for i in agency_pref]].to_array( + dim='agency') # get wmo value val_pref = ds['wmo_'+v] @@ -630,29 +649,26 @@ def map_val_to_ix(a): # add to dict pref_vals[v] = val_pref - - ## THESE CANNOT BE MISSING SO DROP - ## IF EITHER MISSING + # THESE CANNOT BE MISSING SO DROP + # IF EITHER MISSING valid = pref_vals['wind'].notnull() & pref_vals['pres'].notnull() if not valid.any(): raise NoDataError(missing_necessary_data_warning_str) ds = ds.sel(date_time=valid) - for i in ['wind','pres']: + for i in ['wind', 'pres']: pref_vals[i] = pref_vals[i].sel(date_time=valid) - - ## GET RMW and ROCI - ## (these can be missing) - for r in ['rmw','roci']: - order = ['{}_{}'.format(i,r) for i in agency_pref if - '{}_{}'.format(i,r) in ds.data_vars.keys()] + # GET RMW and ROCI + # (these can be missing) + for r in ['rmw', 'roci']: + order = ['{}_{}'.format(i, r) for i in agency_pref if + '{}_{}'.format(i, r) in ds.data_vars.keys()] vals = ds[order].to_array(dim='agency') best_ix = vals.notnull().argmax(dim='agency') val_pref = vals.isel(agency=best_ix) pref_vals[r] = val_pref - - ## CONVERT TO GEOCLAW FORMAT + # CONVERT TO GEOCLAW FORMAT # assign basin to be the basin where track originates # in case track moves across basins @@ -664,36 +680,40 @@ def map_val_to_ix(a): self.t = [] for d in ds.time: t = d.dt - self.t.append(datetime.datetime(t.year,t.month,t.day,t.hour,t.minute,t.second)) + self.t.append(datetime.datetime(t.year, t.month, + t.day, t.hour, t.minute, t.second)) - ## events + # events self.event = ds.usa_record.values.astype(str) - ## time offset - if (self.event=='L').any(): + # time offset + if (self.event == 'L').any(): # if landfall, use last landfall - self.time_offset = numpy.array(self.t)[self.event=='L'][-1] + self.time_offset = np.array(self.t)[self.event == 'L'][-1] else: - #if no landfall, use last time of storm + # if no landfall, use last time of storm self.time_offset = self.t[-1] # Classification, note that this is not the category of the storm self.classification = ds.usa_status.values - self.eye_location = numpy.array([ds.lon,ds.lat]).T + self.eye_location = np.array([ds.lon, ds.lat]).T # Intensity information - for now, including only common, basic intensity # info. # TODO: add more detailed info for storms that have it - self.max_wind_speed = units.convert(pref_vals['wind'],'knots','m/s').where(pref_vals['wind'].notnull(),-1).values - self.central_pressure = units.convert(pref_vals['pres'],'mbar','Pa').where(pref_vals['pres'].notnull(),-1).values - self.max_wind_radius = units.convert(pref_vals['rmw'],'nmi','m').where(pref_vals['rmw'].notnull(),-1).values - self.storm_radius = units.convert(pref_vals['roci'],'nmi','m').where(pref_vals['roci'].notnull(),-1).values + self.max_wind_speed = units.convert( + pref_vals['wind'], 'knots', 'm/s').where(pref_vals['wind'].notnull(), -1).values + self.central_pressure = units.convert(pref_vals['pres'], 'mbar', 'Pa').where( + pref_vals['pres'].notnull(), -1).values + self.max_wind_radius = units.convert(pref_vals['rmw'], 'nmi', 'm').where( + pref_vals['rmw'].notnull(), -1).values + self.storm_radius = units.convert(pref_vals['roci'], 'nmi', 'm').where( + pref_vals['roci'].notnull(), -1).values # warn if you have missing vals for RMW or ROCI if (self.max_wind_radius.max()) == -1 or (self.storm_radius.max() == -1): warnings.warn(missing_data_warning_str) - def read_jma(self, path, verbose=False): r"""Read in JMA formatted storm file @@ -726,13 +746,13 @@ def read_jma(self, path, verbose=False): # Parse data block self.t = [] - self.event = numpy.empty(num_lines, dtype=str) - self.classification = numpy.empty(num_lines, dtype=str) - self.eye_location = numpy.empty((num_lines, 2)) - self.max_wind_speed = numpy.empty(num_lines) - self.central_pressure = numpy.empty(num_lines) - self.max_wind_radius = numpy.empty(num_lines) - self.storm_radius = numpy.empty(num_lines) + self.event = np.empty(num_lines, dtype=str) + self.classification = np.empty(num_lines, dtype=str) + self.eye_location = np.empty((num_lines, 2)) + self.max_wind_speed = np.empty(num_lines) + self.central_pressure = np.empty(num_lines) + self.max_wind_radius = np.empty(num_lines) + self.storm_radius = np.empty(num_lines) for (i, line) in enumerate(data_block): if len(line) == 0: break @@ -754,13 +774,14 @@ def read_jma(self, path, verbose=False): # Intensity information - current the radii are not directly given # Available data includes max/min of radius of winds of 50 and # 30 kts instead - self.central_pressure[i] = units.convert(float(data[5]), 'hPa', 'Pa') - self.max_wind_speed[i] = units.convert(float(data[6]), 'knots', 'm/s') + self.central_pressure[i] = units.convert( + float(data[5]), 'hPa', 'Pa') + self.max_wind_speed[i] = units.convert( + float(data[6]), 'knots', 'm/s') warnings.warn(missing_data_warning_str) self.max_wind_radius[i] = -1 self.storm_radius[i] = -1 - def read_imd(self, path, verbose=False): r"""Extract relevant hurricane data from IMD file and update storm fields with proper values. @@ -774,7 +795,6 @@ def read_imd(self, path, verbose=False): "implemented yet but is planned for a ", "future release.")) - def read_tcvitals(self, path, verbose=False): r"""Extract relevant hurricane data from TCVITALS file and update storm fields with proper values. @@ -801,12 +821,12 @@ def read_tcvitals(self, path, verbose=False): # Central_pressure - convert from mbar to Pa - 100.0 # Radius of last isobar contour - convert from km to m - 1000.0 self.t = [] - self.classification = numpy.empty(num_lines, dtype=str) - self.eye_location = numpy.empty((num_lines, 2)) - self.max_wind_speed = numpy.empty(num_lines) - self.central_pressure = numpy.empty(num_lines) - self.max_wind_radius = numpy.empty(num_lines) - self.storm_radius = numpy.empty(num_lines) + self.classification = np.empty(num_lines, dtype=str) + self.eye_location = np.empty((num_lines, 2)) + self.max_wind_speed = np.empty(num_lines) + self.central_pressure = np.empty(num_lines) + self.max_wind_radius = np.empty(num_lines) + self.storm_radius = np.empty(num_lines) for (i, data) in enumerate(data_block): # End at an empty lines - skips lines at the bottom of a file @@ -836,13 +856,14 @@ def read_tcvitals(self, path, verbose=False): # Intensity Information self.max_wind_speed[i] = float(data[12]) - self.central_pressure[i] = units.convert(float(data[9]), 'mbar', 'Pa') + self.central_pressure[i] = units.convert( + float(data[9]), 'mbar', 'Pa') self.max_wind_radius[i] = units.convert(float(data[13]), 'km', 'm') self.storm_radius[i] = units.convert(float(data[11]), 'km', 'm') - # ========================================================================= # Write Routines + def write(self, path, file_format="geoclaw", **kwargs): r"""Write out the storm data to *path* in format *file_format* @@ -865,105 +886,120 @@ def write(self, path, file_format="geoclaw", **kwargs): getattr(self, 'write_%s' % file_format.lower())(path, **kwargs) - def write_geoclaw(self, path, verbose=False, max_wind_radius_fill=None, - storm_radius_fill=None): + def write_geoclaw(self, path, force=False, skip=True, verbose=False, + fill_dict={}, **kwargs): r"""Write out a GeoClaw formatted storm file GeoClaw storm files are read in by the GeoClaw Fortran code. :Input: - *path* (string) Path to the file to be written. + - *skip* (bool) Skip a time if NaNs are found and are not replaced. + Default is `True`. + - *force* (bool) Force output of storm even if there is missing data. + Default is `False`. - *verbose* (bool) Print out additional information when writing. - - *max_wind_radius_fill* (func) Function that can be used to fill in - missing data for `max_wind_radius` values. This defaults to simply - setting the value to -1. The function signature should be - `max_wind_radius(t, storm)` where t is the time of the forecast and - `storm` is the storm object. Note that if this or `storm_radius` - field remains -1 that this data line will be assumed to be redundant - and not be written out. - - *storm_radius_fill* (func) Function that can be used to fill in - missing data for `storm_radius` values. This defaults to simply - setting the value to -1. The function signature should be - `storm_radius(t, storm)` where t is the time of the forecast and - `storm` is the storm object. Note that if this or `max_wind_radius` - field remains -1 that this data line will be assumed to be redundant - and not be written + Default is `False`. + - *fill_dict* (dict) Dictionary of functions to use to fill in missing + data represented by NaNs. The keys are the field to be filled and + the function signature should be `my_func(t, storm)` where t is the + time of the forecast and `storm` is the storm object. If the + field remains a NaN or a function is not provided these lines will + be assumed redundant and will be ommitted. Note that the older + keyword arguments are put in this dictionary. Currently the one + default function is for `storm_radius`, which sets the value to + 500 km. """ - if max_wind_radius_fill is None: - max_wind_radius_fill = lambda t, storm: -1 - if storm_radius_fill is None: - storm_radius_fill = lambda t, storm: -1 - - # Create list for output - # Leave this first line blank as we need to count the actual valid lines - # that will be left in the file below + # If a filling function is not provided we will provide some defaults + fill_dict.update({"storm_radius": lambda t, storm: 500e3}) + # Handle older interface that had specific fill functions + if "max_wind_radius_fill" in kwargs.keys(): + fill_dict.update( + {"max_wind_radius": kwargs['max_wind_radius_fill']}) + if "storm_radius_fill" in kwargs.keys(): + fill_dict.update({"storm_radius": kwargs['storm_radius_fill']}) + + # Loop through each line of data and if the line is valid, perform the + # necessary work to write it out. Otherwise either raise an exception + # or skip it num_casts = 0 - data_string = [""] - if self.time_offset is None: - data_string.append("None") - self.time_offset = self.t[0] - else: - data_string.append("%s\n\n" % self.time_offset.isoformat()) + data = [] for n in range(len(self.t)): - # Remove duplicate times - if n > 0: - if self.t[n] == self.t[n - 1]: - continue - - format_string = ("{:19,.8e} " * 7)[:-1] + "\n" - data = [] - if not isinstance(self.time_offset, float): - data.append((self.t[n] - self.time_offset).total_seconds()) - else: - data.append(self.t[n] - self.time_offset) - data.append(self.eye_location[n, 0]) - data.append(self.eye_location[n, 1]) - - if self.max_wind_speed[n] == -1: + if self.t[n] == self.t[n - 1]: + # Skip this time continue - data.append(self.max_wind_speed[n]) - - # Allow custom function to set max wind radius if not - # available - if self.max_wind_radius[n] == -1: - new_wind_radius = max_wind_radius_fill(self.t[n], self) - if new_wind_radius == -1: - continue - else: - data.append(new_wind_radius) - else: - data.append(self.max_wind_radius[n]) - if self.central_pressure[n] == -1: + # Check each value we need for this time to make sure it is valid + valid = True + for name in ["max_wind_speed", "central_pressure", + "max_wind_radius", "storm_radius"]: + if np.isnan(getattr(self, name)[n]): + if name in fill_dict.keys(): + # Fill value with function provided + getattr(self, name)[n] = fill_dict[name]( + self.t[n], self) + elif skip: + # Skip this line + valid = False + if verbose: + # Just warn that a NaN was found but continue + msg = ("*** WARNING: The value {} at {} is a " + + "NaN. Skipping this line.") + warnings.warn(msg.format(name, self.t[n])) + elif not force: + # If we are not asked to force to write raise an + # exception given the NaN + msg = ("The value {} at {} is a NaN and the storm " + + "will not be written in GeoClaw format. If " + + "you want to fill in the value provide a " + + "function or set `force=True`.") + raise ValueError(msg.format(name, self.t[n])) + if not valid: continue - data.append(self.central_pressure[n]) - # Allow custom function to set storm radius if not available - if self.storm_radius[n] == -1: - new_storm_radius = storm_radius_fill(self.t[n], self) - if new_storm_radius == -1: - continue - else: - data.append(new_storm_radius) - else: - data.append(self.storm_radius[n]) - - data_string.append(format_string.format(*data)) + # Succeeded, add this time to the output num_casts += 1 + data.append(np.empty(7)) + # If we do not have a time offset use the first valid row as the + # offset time + if self.time_offset is None: + self.time_offset = self.t[n] - # Write to actual file now that we know exactly how many lines it will - # contain + # Time + if not isinstance(self.time_offset, float): + data[-1][0] = (self.t[n] - self.time_offset).total_seconds() + else: + data[-1][0] = self.t[n] - self.time_offset + # Eye-location + data[-1][1:3] = self.eye_location[n, :] + # Max wind speed + data[-1][3] = self.max_wind_speed[n] + # Max wind radius + data[-1][4] = self.max_wind_radius[n] + # Central pressure + data[-1][5] = self.central_pressure[n] + # Outer storm radius + data[-1][6] = self.storm_radius[n] + + # Write out file + format_string = ("{:19,.8e} " * 7)[:-1] + "\n" try: - # Update number of forecasts here - data_string[0] = "%s\n" % num_casts with open(path, "w") as data_file: - for data_line in data_string: - data_file.write(data_line) + # Write header + data_file.write(f"{num_casts}\n") + if isinstance(self.time_offset, datetime.datetime): + data_file.write(f"{self.time_offset.isoformat()}\n\n") + else: + data_file.write(f"{str(self.time_offset)}\n\n") + + # Write data lines + for line in data: + data_file.write(format_string.format(*line)) except Exception as e: - # Remove possibly partially generated file if not successful + # If an exception occurs clean up a partially generated file if os.path.exists(path): os.remove(path) raise e @@ -981,24 +1017,24 @@ def write_atcf(self, path, verbose=False): with open(path, 'w') as data_file: for n in range(len(self.t)): data_file.write("".join((", " * 2, - "%s" % seconds2date(self.t[n]), - ", " * 4, - "%s" % (int(self.eye_location[n, 0] * - 10.0)), - ", ", - "%s" % (int(self.eye_location[n, 1] * - 10.0)), - ", ", - "%s" % self.max_wind_speed[n], - ", ", - "%s" % self.central_pressure[n], - ", ", - ", " * 8, - "%s" % self.storm_radius[n], - ", ", - "%s" % self.max_wind_radius[n], - ", " * 10, - "\n"))) + "%s" % seconds2date(self.t[n]), + ", " * 4, + "%s" % (int(self.eye_location[n, 0] * + 10.0)), + ", ", + "%s" % (int(self.eye_location[n, 1] * + 10.0)), + ", ", + "%s" % self.max_wind_speed[n], + ", ", + "%s" % self.central_pressure[n], + ", ", + ", " * 8, + "%s" % self.storm_radius[n], + ", ", + "%s" % self.max_wind_radius[n], + ", " * 10, + "\n"))) except Exception as e: # Remove possiblly partially generated file if not successful if os.path.exists(path): @@ -1026,34 +1062,34 @@ def write_hurdat(self, path, verbose=False): # Convert latitude to proper Hurdat format e.g 12.0N if latitude > 0: - latitude = str(numpy.abs(latitude)) + 'N' + latitude = str(np.abs(latitude)) + 'N' else: - latitude = str(numpy.abs(latitude)) + 'S' + latitude = str(np.abs(latitude)) + 'S' # Convert longitude to proper Hurdat format e.g 12.0W if longitude > 0: - longitude = str(numpy.abs(longitude)) + 'E' + longitude = str(np.abs(longitude)) + 'E' else: - longitude = str(numpy.abs(longitude)) + 'W' + longitude = str(np.abs(longitude)) + 'W' data_file.write("".join(("%s" % self.seconds2date( - self.t[n])[0:-2], - "%s00" % self.seconds2date( - self.t[n])[-2:], - ", " * 3, - "%s" % (latitude), - ", ", - "%s" % (longitude), - ", ", - "%s" % self.max_wind_speed[n], - ", ", - "%s" % self.central_pressure[n], - ", ", - "%s" % self.storm_radius[n], - ", ", - "%s" % self.max_wind_radius[n], - ", " * 10, - "\n"))) + self.t[n])[0:-2], + "%s00" % self.seconds2date( + self.t[n])[-2:], + ", " * 3, + "%s" % (latitude), + ", ", + "%s" % (longitude), + ", ", + "%s" % self.max_wind_speed[n], + ", ", + "%s" % self.central_pressure[n], + ", ", + "%s" % self.storm_radius[n], + ", ", + "%s" % self.max_wind_radius[n], + ", " * 10, + "\n"))) except Exception as e: # Remove possiblly partially generated file if not successful if os.path.exists(path): @@ -1073,23 +1109,23 @@ def write_jma(self, path, verbose=False): with open(path, 'w') as data_file: for n in range(self.t.shape[0]): data_file.write("".join(("%s" % self.seconds2date(self.t[n]), - " " * 4, - "%s" % (int(self.eye_location[n, 0] * - 10.0)), - ", ", - "%s" % (int(self.eye_location[n, 1] * - 10.0)), - ", ", - "%s" % self.max_wind_speed[n], - ", ", - "%s" % self.central_pressure[n], - ", ", - ", " * 8, - "%s" % self.storm_radius[n], - ", ", - "%s" % self.max_wind_radius[n], - ", " * 10, - "\n"))) + " " * 4, + "%s" % (int(self.eye_location[n, 0] * + 10.0)), + ", ", + "%s" % (int(self.eye_location[n, 1] * + 10.0)), + ", ", + "%s" % self.max_wind_speed[n], + ", ", + "%s" % self.central_pressure[n], + ", ", + ", " * 8, + "%s" % self.storm_radius[n], + ", ", + "%s" % self.max_wind_radius[n], + ", " * 10, + "\n"))) except Exception as e: # Remove possiblly partially generated file if not successful if os.path.exists(path): @@ -1120,7 +1156,7 @@ def write_tcvitals(self, path, verbose=False): # ========================================================================= # Other Useful Routines - + def category(self, categorization="NHC", cat_names=False): r"""Categorizes storm based on relevant storm data @@ -1145,7 +1181,7 @@ def category(self, categorization="NHC", cat_names=False): if categorization.upper() == "BEAUFORT": # Beaufort scale below uses knots speeds = units.convert(self.max_wind_speed, "m/s", "knots") - category = (numpy.zeros(speeds.shape) + + category = (np.zeros(speeds.shape) + (speeds >= 1) * (speeds < 4) * 1 + (speeds >= 4) * (speeds < 7) * 2 + (speeds >= 7) * (speeds < 11) * 3 + @@ -1158,16 +1194,16 @@ def category(self, categorization="NHC", cat_names=False): (speeds >= 48) * (speeds < 56) * 10 + (speeds >= 56) * (speeds < 64) * 11 + (speeds >= 64) * 12) - cat_map = { 0: "Calm", - 1: "Light air", - 2: "Light breeze", - 3: "Gentle breeze", - 4: "Moderate breeze", - 5: "Fresh breeze", - 6: "Strong breeze", - 7: "High wind", - 8: "Gale", - 9: "Strong gale", + cat_map = {0: "Calm", + 1: "Light air", + 2: "Light breeze", + 3: "Gentle breeze", + 4: "Moderate breeze", + 5: "Fresh breeze", + 6: "Strong breeze", + 7: "High wind", + 8: "Gale", + 9: "Strong gale", 10: "Whole gale", 11: "Violent storm", 12: "Hurricane"} @@ -1177,7 +1213,7 @@ def category(self, categorization="NHC", cat_names=False): # Definitely not in the correct format now # TODO: Add TD and TS designations speeds = units.convert(self.max_wind_speed, "m/s", "knots") - category = (numpy.zeros(speeds.shape) + + category = (np.zeros(speeds.shape) + (speeds < 30) * -1 + (speeds >= 64) * (speeds < 83) * 1 + (speeds >= 83) * (speeds < 96) * 2 + @@ -1185,12 +1221,12 @@ def category(self, categorization="NHC", cat_names=False): (speeds >= 113) * (speeds < 135) * 4 + (speeds >= 135) * 5) cat_map = {-1: "Tropical Depression", - 0: "Tropical Storm", - 1: "Category 1 Hurricane", - 2: "Category 2 Hurricane", - 3: "Category 3 Hurricane", - 4: "Category 4 Hurricane", - 5: "Category 5 Hurricane"} + 0: "Tropical Storm", + 1: "Category 1 Hurricane", + 2: "Category 2 Hurricane", + 3: "Category 3 Hurricane", + 4: "Category 4 Hurricane", + 5: "Category 5 Hurricane"} elif categorization.upper() == "JTWC": raise NotImplementedError("JTWC categorization not implemented.") @@ -1314,17 +1350,17 @@ def fill_rad_w_other_source(t, storm_targ, storm_fill, var, interp_kwargs={}): print("fill_rad_w_other_source currently requires xarray to work.") raise e - fill_da = xr.DataArray(getattr(storm_fill,var), - coords = {'t': getattr(storm_fill,'t')}, - dims = ('t',)) + fill_da = xr.DataArray(getattr(storm_fill, var), + coords={'t': getattr(storm_fill, 't')}, + dims=('t',)) # convert -1 to nan - fill_da = fill_da.where(fill_da>0,numpy.nan) + fill_da = fill_da.where(fill_da > 0, np.nan) # if not all missing, try using storm_fill to fill if fill_da.notnull().any(): - #remove duplicates + # remove duplicates fill_da = fill_da.groupby('t').first() # remove NaNs @@ -1335,19 +1371,19 @@ def fill_rad_w_other_source(t, storm_targ, storm_fill, var, interp_kwargs={}): # try replacing with storm_fill # (assuming atcf has more data points than ibtracs) - if not numpy.isnan(fill_interp): + if not np.isnan(fill_interp): return fill_interp # next, try just interpolating other ibtracs values - targ_da = xr.DataArray(getattr(storm_targ,var), - coords = {'t': getattr(storm_targ,'t')}, - dims = ('t',)) - targ_da = targ_da.where(targ_da>0,numpy.nan) + targ_da = xr.DataArray(getattr(storm_targ, var), + coords={'t': getattr(storm_targ, 't')}, + dims=('t',)) + targ_da = targ_da.where(targ_da > 0, np.nan) if targ_da.notnull().any(): targ_da = targ_da.groupby('t').first() targ_da = targ_da.dropna('t') targ_interp = targ_da.interp(t=[t], kwargs=interp_kwargs).item() - if not numpy.isnan(targ_interp): + if not np.isnan(targ_interp): return targ_interp # if nothing worked, return the missing value (-1) @@ -1392,19 +1428,23 @@ def make_multi_structure(path): fileWrite.writelines(line) else: fileWrite.close() - stormDict[curTime].update({curTrack: Storm(path=os.path.join(os.path.expandvars(os.getcwd()), "Clipped_ATCFs", curTime, curTrack), file_format="ATCF")}) + stormDict[curTime].update({curTrack: Storm(path=os.path.join(os.path.expandvars( + os.getcwd()), "Clipped_ATCFs", curTime, curTrack), file_format="ATCF")}) curTrack = lineArr[4] - fileWrite = open("Clipped_ATCFs/" + curTime + "/" + curTrack, 'w') + fileWrite = open("Clipped_ATCFs/" + + curTime + "/" + curTrack, 'w') fileWrite.writelines(line) else: if curTime != "test": fileWrite.close() - stormDict[curTime].update({curTrack: Storm(path=os.path.join(os.path.expandvars(os.getcwd()), "Clipped_ATCFs", curTime, curTrack), file_format="ATCF")}) + stormDict[curTime].update({curTrack: Storm(path=os.path.join(os.path.expandvars( + os.getcwd()), "Clipped_ATCFs", curTime, curTrack), file_format="ATCF")}) curTime = lineArr[2] curTrack = lineArr[4] stormDict[curTime] = {} os.mkdir("Clipped_ATCFs/" + curTime) - fileWrite = open("Clipped_ATCFs/" + curTime + "/" + curTrack, 'w') + fileWrite = open("Clipped_ATCFs/" + + curTime + "/" + curTrack, 'w') fileWrite.writelines(line) return stormDict