Merge pull request NHERI-SimCenter#398 from yisangriB/master

sy - polishing hurricane generator

Author: fmckenna

modules/performRegionalEventSimulation/regionalWindField/ComputeIntensityMeasure.py

@@ -156,9 +156,9 @@ def simulate_storm_cpp(  # noqa: C901, D103
             ][1]
         # updating landfall properties
         scenario_info['Storm']['LandingAngle'] = scenario_data[0]['CycloneParam'][2]
-        scenario_info['Storm']['Pressure'] = scenario_data[0]['CycloneParam'][3]
-        scenario_info['Storm']['Speed'] = scenario_data[0]['CycloneParam'][4]
-        scenario_info['Storm']['Radius'] = scenario_data[0]['CycloneParam'][5]
+        scenario_info['Storm']['Pressure'] = scenario_data[0]['Pressure']
+        scenario_info['Storm']['Speed'] = scenario_data[0]['Speed']
+        scenario_info['Storm']['Radius'] = scenario_data[0]['Radius']
 
         config = {'Scenario': scenario_info, 'Event': event_info}
         abs_path_config = os.path.abspath(os.path.join(input_dir, 'SimuConfig.json'))  # noqa: PTH100, PTH118
@@ -191,27 +191,34 @@ def simulate_storm_cpp(  # noqa: C901, D103
             )
             df = pd.DataFrame.from_dict(  # noqa: PD901
                 {
-                    'Lat': scenario_data[0]['TrackSimu'],
+                    'Lat': scenario_data[0]['TrackSimu']['Latitude'],
+                    'Lon': scenario_data[0]['TrackSimu']['Longitude'],
                 }
             )
             df.to_csv(abs_path_latw, sep=',', header=False, index=False)
+
         if scenario_info['Generator'] == 'SimulationHist':
             df = pd.DataFrame.from_dict(  # noqa: PD901
                 {
-                    'Lat': scenario_data[0]['TrackSimu'],
+                    'Lat': scenario_data[0]['TrackSimu']['Latitude'],
+                    'Lon': scenario_data[0]['TrackSimu']['Longitude'],
                 }
             )
+
             df.to_csv(abs_path_latw, sep=',', header=False, index=False)
         # terrain file
         if 'Terrain' in scenario_info.keys():  # noqa: SIM118
+            file_name = os.path.basename(scenario_info['Terrain'])
+
             abs_path_terrain = os.path.abspath(  # noqa: PTH100
-                os.path.join(input_dir, scenario_info['Terrain'])  # noqa: PTH118
+                os.path.join(input_dir, file_name)  # noqa: PTH118
             )
         else:
             # default terrain z0 = 0.01 everywhere for the defined domain
             abs_path_terrain = os.path.abspath(  # noqa: PTH100
                 os.path.join(input_dir, 'DefaultTerrain.geojson')  # noqa: PTH118
             )
+
             dict_dt = {
                 'type': 'FeatureCollection',
                 'features': [
@@ -292,6 +299,7 @@ def simulate_storm_cpp(  # noqa: C901, D103
                 if os.path.exists(output_subdir):  # noqa: PTH110
                     shutil.rmtree(output_subdir)
                 os.makedirs(output_subdir)  # noqa: PTH103
+
                 args = [
                     windsimu_bin,
                     '--config',
@@ -310,19 +318,26 @@ def simulate_storm_cpp(  # noqa: C901, D103
                     output_subdir,
                     '--output',
                     output_subdir,
-                ]
+                ]                
+
 
                 pert_list.append(abs_path_pert)
                 args_list.append(args)
                 odir_list.append(output_subdir)
             # running
             print('ComputeIntensityMeaure: running analysis.')  # noqa: T201
+            print(args_list)  # noqa: T201
             procs_list = [
                 subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)  # noqa: S603
                 for cmd in args_list
             ]
             for proc in procs_list:
-                proc.communicate()
+                stdout, stderr = proc.communicate()  # Wait for the process to finish and get its output
+                if proc.returncode == 0:
+                    print(f"Output of process (PID {proc.pid}): {stdout.decode()}")
+                else:
+                    print(f"Error in process (PID {proc.pid}): {stderr.decode()}")
+                    
             # loading output
             print('ComputeIntensityMeaure: postprocessing simulation data.')  # noqa: T201
             for j in range(num_per_site):
@@ -352,9 +367,9 @@ def simulate_storm_cpp(  # noqa: C901, D103
                 )
                 station_res['PWS']['windspeed'] = df.values.tolist()  # noqa: PD011
                 res.append(station_res)
-                shutil.rmtree(odir_list[j])
+                #shutil.rmtree(odir_list[j])
         # house-keeping
-        os.remove(abs_path_config)  # noqa: PTH107
+        #os.remove(abs_path_config)  # noqa: PTH107
     else:
         print(  # noqa: T201
             'ComputeIntensityMeasure: currently only supporting LinearAnalytical model'
@@ -419,13 +434,13 @@ def convert_wind_speed(event_info, simu_res):  # noqa: D103
             zg_t = 274.32
         # conversion
         pws_raw = interp_wind_by_height(pws_raw, measure_height, reference_height)
-        print(np.max(pws_raw))  # noqa: T201
+        # print(np.max(pws_raw))  # noqa: T201
         # computing gradient-height wind speed
         pws_tmp = pws_raw * (zg / reference_height) ** (1.0 / alpha)
         # converting exposure
         pws_tmp = pws_tmp * (reference_height / zg_t) ** (1.0 / alpha_t)
         pws = pws_tmp * gust_factor_ESDU(gust_duration_simu, gust_duration)
-        print(np.max(pws))  # noqa: T201
+        # print(np.max(pws))  # noqa: T201
         # appending to pws_mr
         pws_mr.append(pws)
 

modules/performRegionalEventSimulation/regionalWindField/CreateScenario.py

@@ -47,6 +47,7 @@
 def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noqa: C901, D103
     # Number of scenarios
     source_num = scenario_info.get('Number', 1)
+
     # Directly defining earthquake ruptures
     if scenario_info['Generator'] == 'Simulation':
         # Collecting site locations
@@ -81,9 +82,16 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
                 header=None,
                 index_col=None,
             )
-            track_simu = df.iloc[:, 0].values.tolist()  # noqa: PD011
+            #track_simu = df.iloc[:, 0].values.tolist()  # noqa: PD011
+            track_simu = {
+                'Latitude': df.iloc[:, 0].values.tolist(),  # noqa: PD011
+                'Longitude': df.iloc[:, 1].values.tolist(),  # noqa: PD011
+            }
         else:
-            track_simu = track['Latitude']
+            #track_simu = track['Latitude']
+            track_simu = track
+
+
         # Reading Terrain info (if provided)
         terrain_file = scenario_info.get('Terrain', None)
         if terrain_file:
@@ -97,11 +105,21 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
             param.append(scenario_info['Storm']['Landfall']['Latitude'])
             param.append(scenario_info['Storm']['Landfall']['Longitude'])
             param.append(scenario_info['Storm']['Landfall']['LandingAngle'])
-            param.append(scenario_info['Storm']['Landfall']['Pressure'])
-            param.append(scenario_info['Storm']['Landfall']['Speed'])
-            param.append(scenario_info['Storm']['Landfall']['Radius'])
+            #param.append(scenario_info['Storm']['Landfall']['Pressure'])
+            #param.append(scenario_info['Storm']['Landfall']['Speed'])
+            #param.append(scenario_info['Storm']['Landfall']['Radius'])
+            landfall_prs = 1013.0 - float(scenario_info['Storm']['Landfall']['Pressure']) 
+            landfall_spd = float(scenario_info['Storm']['Landfall']['Speed'])* 1.852  # from kts to km/h
+            landfall_rad = float(scenario_info['Storm']['Landfall']['Radius'])* 1.852 # from nmile to km
         except:  # noqa: E722
             print('CreateScenario: please provide all needed landfall properties.')  # noqa: T201
+
+
+
+        track_simu,dummy,dummy,dummy = interpolate_track(track_simu)
+        prs_list, spd_list, rad_list = landfall2track(landfall_prs, landfall_spd, landfall_rad, track_simu)
+
+
         # Monte-Carlo
         # del_par = [0, 0, 0] # default
         # Parsing mesh configurations
@@ -117,6 +135,9 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
                     i: {
                         'Type': 'Wind',
                         'CycloneParam': param,
+                        'Pressure': prs_list,
+                        'Speed': spd_list,
+                        'Radius': rad_list,
                         'StormTrack': track,
                         'StormMesh': mesh_info,
                         'Terrain': terrain_data,
@@ -212,6 +233,7 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
             tmploc = dist2land.index(
                 0
             )  # the first landing point in case the storm sway back and forth
+        
         # simulation track
         track_simu_file = scenario_info['Storm'].get('TrackSimu', None)
         if track_simu_file:
@@ -221,20 +243,30 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
                     header=None,
                     index_col=None,
                 )
-                track_simu = df.iloc[:, 0].values.tolist()  # noqa: PD011
+                #track_simu = df.iloc[:, 0].values.tolist()  # noqa: PD011
+
+                track_simu = {
+                    'Latitude': df.iloc[:, 0].values.tolist(),  # noqa: PD011
+                    'Longitude': df.iloc[:, 1].values.tolist(),  # noqa: PD011
+                }
+
             except:  # noqa: E722
                 print(  # noqa: T201
                     'CreateScenario: warning - TrackSimu file not found, using the full track.'
                 )
-                track_simu = track_lat
+                #track_simu = track_lat
+                track_simu = track
         else:
             print(  # noqa: T201
                 'CreateScenario: warning - no truncation defined, using the full track.'
             )
             # tmp = track_lat
             # track_simu = tmp[max(0, tmploc - 5): len(dist2land) - 1]
             # print(track_simu)
-            track_simu = track_lat
+            #track_simu = track_lat
+            track_simu = track
+
+
         # Reading data
         try:
             landfall_lat = float(df_chs[('USA_LAT', 'degrees_north')].iloc[tmploc])  # noqa: RUF031, RUF100
@@ -249,46 +281,71 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
             print('CreateScenario: error - no landing angle is found.')  # noqa: T201
         if landfall_ang > 180.0:  # noqa: PLR2004
             landfall_ang = landfall_ang - 360.0
-        landfall_prs = (
-            1013.0
-            - np.min(
-                [
-                    float(x)
-                    for x in df_chs[('USA_PRES', 'mb')]  # noqa: PD011, RUF031, RUF100
-                    .iloc[tmploc - 5 :]
-                    .values.tolist()
-                    if x != ' '
-                ]
-            )
-        )
-        landfall_spd = (
-            float(df_chs[('STORM_SPEED', 'kts')].iloc[tmploc]) * 0.51444  # noqa: RUF031, RUF100
-        )  # convert knots/s to km/s
+        # prs_list = (
+        #     1013.0
+        #     - np.min(
+        #         [
+        #             float(x)
+        #             for x in df_chs[('USA_PRES', 'mb')]  # noqa: PD011, RUF031, RUF100
+        #             .iloc[:]
+        #             .values.tolist()
+        #             if x != ' '
+        #         ]
+        #     )
+        # )
+
+        air_pressure = df_chs[('USA_PRES', 'mb')].iloc[:].values.tolist() 
+        prs_list = []
+        for x in air_pressure:
+            if x == ' ':
+                prs_list += [0] # no record == insignificant dp
+            else:
+                prs_list += [1013.0  - float(x)]
+
+
+        # prs_list = [1013.0 - float(x) for x in df_chs[('USA_PRES', 'mb')].iloc[:].values.tolist() ]
+
+        # spd_list = (
+        #     # float(df_chs[('STORM_SPEED', 'kts')].iloc[tmploc]) * 0.51444  # noqa: RUF031, RUF100 
+        #     float(df_chs[('STORM_SPEED', 'kts')].iloc[:]) * 0.00051444  # noqa: RUF031, RUF100  
+        # )  # convert knots to km/h // TODO: check if it was intended to be km/s or m/s. If it should be km/s this is correct.
+
+        spd_list = [float(x)* 1.852   for x in df_chs[('STORM_SPEED', 'kts')].iloc[:].values.tolist() ] 
+
         try:
-            landfall_rad = (
-                float(df_chs[('USA_RMW', 'nmile')].iloc[tmploc]) * 1.60934  # noqa: RUF031, RUF100
-            )  # convert nmile to km
+            # rad_list = (
+            #     float(df_chs[('USA_RMW', 'nmile')].iloc[:]) * 1.852  # noqa: RUF031, RUF100   
+            # )  # convert nmile to km // sy - from nmile to km is 1.852 not 1.609
+            rad_list = [float(x)* 1.852 for x in df_chs[('USA_RMW', 'nmile')].iloc[:].values.tolist() ]
+
         except:  # noqa: E722
             # No available radius of maximum wind is found
             print('CreateScenario: warning - switching to REUNION_RMW.')  # noqa: T201
             try:
                 # If the default option (USA_RMW) is not available, switching to REUNION_RMW
-                landfall_rad = (
-                    float(df_chs[('REUNION_RMW', 'nmile')].iloc[tmploc]) * 1.60934  # noqa: RUF031, RUF100
-                )  # convert nmile to km
+                # rad_list = (
+                #     float(df_chs[('REUNION_RMW', 'nmile')].iloc[:]) * 1.852# noqa: RUF031, RUF100  
+                # )  # convert nmile to km // sy - from nmile to km is 1.852 not 1.609
+                rad_list = [float(x)* 1.852 for x in df_chs[('REUNION_RMW', 'nmile')].iloc[:].values.tolist() ]
+
             except:  # noqa: E722
                 # No available radius of maximum wind is found
                 print(  # noqa: T201
                     'CreateScenario: warning - no available radius of maximum wind is found, using a default 50 km.'
                 )
-                landfall_rad = 50
+                rad_list = 50
+
+
+        track_simu, prs_list, spd_list, rad_list = interpolate_track(track_simu, prs_list, spd_list, rad_list )
+
+
         param = []
         param.append(landfall_lat)
         param.append(landfall_lon)
         param.append(landfall_ang)
-        param.append(landfall_prs)
-        param.append(landfall_spd)
-        param.append(landfall_rad)
+        #param.append(prs_list)
+        #param.append(spd_list)
+        #param.append(rad_list)
         # Monte-Carlo
         # del_par = [0, 0, 0] # default
         # Parsing mesh configurations
@@ -304,6 +361,9 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
                     i: {
                         'Type': 'Wind',
                         'CycloneParam': param,
+                        'Pressure': prs_list,
+                        'Speed': spd_list,
+                        'Radius': rad_list,
                         'StormTrack': track,
                         'StormMesh': mesh_info,
                         'Terrain': terrain_data,
@@ -318,3 +378,62 @@ def create_wind_scenarios(scenario_info, event_info, stations, data_dir):  # noq
 
     else:
         print('CreateScenario: currently only supporting Simulation generator.')  # noqa: T201, RET503
+
+
+
+def interpolate_track(track_simu, prs_list=None, spd_list=None, rad_list=None):  # noqa: C901, D103
+
+    # if maximum of (dlat, dlong) is larger than the threshold value, we will add more intermediate points
+
+    track_lat = track_simu['Latitude']
+    track_long = track_simu['Longitude']
+    thres = 0.1
+
+    track_lat_new = [track_lat[0]]
+    track_long_new = [track_long[0]]
+
+    if prs_list is not None:
+        prs_new =  [prs_list[0]]
+        spd_new =  [spd_list[0]]
+        rad_new =  [rad_list[0]]
+    else:
+        prs_new, spd_new, rad_new = None, None, None
+        
+    for nn in range(len(track_lat)-1):
+        dlat = (track_lat[nn+1]-track_lat[nn])
+        dlon = (track_long[nn+1]-track_long[nn])
+
+        if prs_list is not None:
+            dp = (prs_list[nn+1]-prs_list[nn])
+            ds = (spd_list[nn+1]-spd_list[nn])
+            dr = (rad_list[nn+1]-rad_list[nn])
+
+
+        npoints = int(max(dlat, dlon)//thres) + 1 # // is same as Ceil: This should be 1 if interval is sufficiently small
+        for np in range(npoints):
+            track_lat_new += [track_lat[nn] + dlat/npoints*(np+1)]
+            track_long_new += [track_long[nn] + dlon/npoints*(np+1)]
+
+            if prs_list is not None:
+                prs_new += [prs_list[nn] + dp/npoints*(np+1)]
+                spd_new += [spd_list[nn] + ds/npoints*(np+1)]
+                rad_new += [rad_list[nn] + dr/npoints*(np+1)]
+
+    track_simu_interpolated = {
+                'Latitude': track_lat_new,
+                'Longitude': track_long_new
+            }       
+
+    return track_simu_interpolated, prs_new, spd_new, rad_new
+
+
+def landfall2track(landfall_prs, landfall_spd, landfall_rad, track_simu):
+    
+    # TODO filling model should replace this function
+    #print(landfall_prs)
+    #print(track_simu)
+    prs_list = [landfall_prs] * len(track_simu['Latitude'])
+    psd_list = [landfall_spd] * len(track_simu['Latitude'])
+    rad_list = [landfall_rad] * len(track_simu['Latitude'])
+
+    return prs_list, psd_list, rad_list, 

modules/performRegionalEventSimulation/regionalWindField/HurricaneSimulation.py

@@ -110,6 +110,7 @@
     if scenario_info['Type'] == 'Wind':
         if 'Simulation' in scenario_info['Generator']:
             if scenario_info['ModelType'] == 'LinearAnalyticalPy':
+                # TODO: sy - the below python code need to be updated because now scenarios['track_simu'] contains both latitude and longitude. Only c++ code updated for now
                 # simulating storm
                 storm_simu = simulate_storm(  # noqa: F405
                     scenarios, event_info, 'LinearAnalytical'