Refactor Projection Argument and Time Metadata; Fix Kajiura Filter and Coordinate Arrays

donorModels/donorInterface.py

@@ -13,14 +13,14 @@
 """
 
 
-def get_donorModel(choose_donormodel, donor_output, spatial_resolution, CRS_reference, bathy_file, include_horizontal=False):
+def get_donorModel(choose_donormodel, donor_output, spatial_resolution, projection, bathy_file, include_horizontal=False):
   """
   This function provides the data for the respective donor model.
 
   :param choose_donormodel: donor model 
   :param donor_output: name of the output file or path from the donor model
   :param spatial_resolution: spatial resolution for which the data will be given as a structured mesh (resolution is the same in both horizontal directions)
-  :param CRS_reference: CRS reference coordinates (list of longitude and latitude of lower left corner of the domain)  
+  *param projection: projection parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates
   :param bathy_file: file for the bathymetry (has to be larger than the mesh provided by the donor mesh)
   :param include_horizontal: handle whether to include horizontal deformations (False by default)
   """
@@ -33,10 +33,10 @@ def get_donorModel(choose_donormodel, donor_output, spatial_resolution, CRS_refe
     donor_deformation, donor_x, donor_y, donor_time = get_bingclaw(donor_output, resolution)
     return donor_deformation, donor_x, donor_y, donor_time, donor_bathy, resolution
   if (choose_donormodel == 'seissol'):
-    donor_deformation, donor_x, donor_y, donor_time = get_seissol(donor_output, resolution, CRS_reference, include_horizontal)
+    donor_deformation, donor_x, donor_y, donor_time = get_seissol(donor_output, resolution, projection, include_horizontal)
     return donor_deformation, donor_x, donor_y, donor_time, donor_bathy, resolution
   if (choose_donormodel == 'shaltop'):
-    donor_deformation, donor_x, donor_y, donor_time = get_shaltop(donor_output, resolution, CRS_reference)
+    donor_deformation, donor_x, donor_y, donor_time = get_shaltop(donor_output, resolution, projection)
     return donor_deformation, donor_x, donor_y, donor_time, donor_bathy, resolution
   else:
     raise NotImplementedError("The provided donor model is unknown. Possible donor models include: bingclaw, seissol or shaltop")

donorModels/donorResolution.py

@@ -17,7 +17,6 @@ def get_bathyResolutionlnMeters(bathy_resolution):
   Calculate bathymetry resolution in m (from transformation).
 
   :param bathy_resolution: resolution of bathymetry file in lat/lon
-  :param CRS_reference_coordinates: CRS reference coordinates (list of longitude and latitude of lower left corner of the domain)
 
   """
   # Define CRS  
@@ -51,7 +50,10 @@ def donor_chooseResolution(spatial_resolution, bathy_file):
   except:
     bathy_x = bathy_data.variables['lon']
     bathy_y = bathy_data.variables['lat']
-    bathy = bathy_data.variables['elevation']
+    try:
+      bathy = bathy_data.variables['elevation']
+    except:
+      bathy = bathy_data.variables['z']
 
   # Check spatial resolution and set it equal to the bathymetry resolution if 0
   if (spatial_resolution == 0.0):

donorModels/seissol2exchangegrid.py

@@ -97,19 +97,15 @@ def project_coordinates(x, y, inputCRS):
   :param inputCRS: input CRS
   """
   transformer = Transformer.from_crs(inputCRS, basicCRS, always_xy=True)
-
-  # move x- and y-coordinates so that the lower left corner lies at [0,0]
-  x_tmp = x + np.abs(x[0])
-  y_tmp = y + np.abs(y[0])
-
+
   # transform x-coordinates
-  y_lowerRow = np.repeat(y_tmp[0], len(x))
-  xnew, ynew = transformer.transform(x_tmp, y_lowerRow)
+  y_lowerRow = np.repeat(y[0], len(x))
+  xnew, ynew = transformer.transform(x, y_lowerRow)
   x_proj = xnew  
 
     # transform y-coordinates
-  x_leftColumn = np.repeat(x_tmp[0], len(y))
-  xnew, ynew = transformer.transform(x_leftColumn, y_tmp)
+  x_leftColumn = np.repeat(x[0], len(y))
+  xnew, ynew = transformer.transform(x_leftColumn, y)
   y_proj = ynew
 
   return x_proj, y_proj
@@ -218,10 +214,11 @@ def interpolate_seissol2structured(sx, dx, coord_min, coord_max, inputCRS, inclu
   nTime = sx.ReadNdt()  # number of time steps in the Seissol file
 
   # Choose which data interpolate (only vertical or all components) based on the input handle  
+  is_new_format = 'u3' in sx.ReadAvailableDataFields()
   if (include_horizontal):
-    data = ['u1', 'u2', 'u3']
+    data = ['u1', 'u2', 'u3'] if is_new_format else ['U', 'V', 'W']
   else:
-    data = ['u3']    
+    data = ['u3'] if is_new_format else ['W']   
 
   # Create probe filter and get projected coordinates
   probeFilter, projDataShape, x_proj, y_proj = setUp_grid_interpolation(coord_min, coord_max, dx, inputCRS)
@@ -243,20 +240,16 @@ def interpolate_seissol2structured(sx, dx, coord_min, coord_max, inputCRS, inclu
 
 
 
-def get_seissol(filename, spatial_resolution, CRS_reference_coordinates, include_horizontal):
+def get_seissol(filename, spatial_resolution, projection, include_horizontal):
   """
   Actual part that will be used by the main donor model functionality to get the deformation data.
 
   :param filename: filename for the SeisSol data (has to be an XDMF file)
   :param spatial_resolution: spatial_resolution in meters
-  :param CRS_reference_coordinates: CRS reference coordinates (list of longitude and latitude of lower left corner of the domain)  
+  *param projection: projection parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates
   :param include_horizontal: boolean handle whether the output will only contain the vertical deformation (deformation) or also include the horizontal deformation
   """
 
-  # CRS of the 2d mesh
-  inputCRS = "+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=" + CRS_reference_coordinates[0] + \
-                      " +lat_0=" + CRS_reference_coordinates[1]  
-
   print("Getting output data from SeisSol.\n".center(column_size))
 
   sx = seissolxdmfExtended(filename) # get seissolxdmf from provided XDMF file
@@ -266,6 +259,6 @@ def get_seissol(filename, spatial_resolution, CRS_reference_coordinates, include
   coordinate_min = geom.min(0)
   coordinate_max = geom.max(0) 
 
-  donor_deformation, donor_x, donor_y, donor_time = interpolate_seissol2structured(sx, spatial_resolution, coordinate_min, coordinate_max, inputCRS, include_horizontal)
+  donor_deformation, donor_x, donor_y, donor_time = interpolate_seissol2structured(sx, spatial_resolution, coordinate_min, coordinate_max, projection, include_horizontal)
 
   return donor_deformation, donor_x, donor_y, donor_time

donorModels/shaltop2exchangegrid.py

@@ -281,19 +281,15 @@ def interpolate_shaltop_data(shaltop_deformation, spatial_resolution, shaltop_x,
 
 
 
-def get_shaltop(donor_output_path, spatial_resolution, CRS_reference_coordinates):
+def get_shaltop(donor_output_path, spatial_resolution, projection):
   """
   Main donor model functionality to get the deformation data.
 
   :param donor_output_path: path to the directory where the SHALTOP output is stored
   :param spatial_resolution: spatial_resolution in meters
-  :param CRS_reference_coordinates: CRS reference coordinates (list of longitude and latitude of lower left corner of the domain)
+  *param projection: projection parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates
   """
 
-  # CRS of the 2d mesh
-  inputCRS = "+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=" + CRS_reference_coordinates[0] + \
-                      " +lat_0=" + CRS_reference_coordinates[1]  
-
   print("Getting output data from SHALTOP.\n".center(column_size))
 
   shaltop_deformation, shaltop_x, shaltop_y, donor_time = obtain_shaltop_data(donor_output_path)
@@ -304,7 +300,7 @@ def get_shaltop(donor_output_path, spatial_resolution, CRS_reference_coordinates
 
   donor_deformation, interpolated_x, interpolated_y = interpolate_shaltop_data(shaltop_deformation, spatial_resolution, shaltop_x, shaltop_y, len(donor_time))
 
-  donor_x, donor_y = project_coordinates(interpolated_x, interpolated_y, inputCRS)
+  donor_x, donor_y = project_coordinates(interpolated_x, interpolated_y, projection)
   stop = time.time()    
   print(f"The interpolation took {stop - start} s.\n".center(column_size))
 

exchangeGrid/exchangeGridCreation.py

@@ -54,8 +54,8 @@ def donor2bathyDomain(bathy_file, donor_x, donor_y, donor_deformation):
   dx = donor_x[1] - donor_x[0]
 
   # Calculate new number of points for both x- and y-coordinates
-  new_Nx = int((bathy_x[-1] - bathy_x[0]) / dx) 
-  new_Ny = int((bathy_y[-1] - bathy_y[0]) / dx) 
+  new_Nx = int((bathy_x[-1] - bathy_x[0]) / dx) + 1
+  new_Ny = int((bathy_y[-1] - bathy_y[0]) / dx) + 1
 
   # Calculate indices at which the "old"/donor data will start and end on the "new"/exchange grid
   donor_x_index_start = int((donor_x[0] - x_min)/dx)

exchangeGrid/filter.py

@@ -197,7 +197,7 @@ def use_kajiura_filter(deformation, bathymetry, spatial_resolution, filtering_de
   Nx = np.shape(deformation)[2]
 
   current_disp = np.zeros((Ny, Nx))
-  current_η = np.zeros((Ny, Nx))
+  current_η_diff = np.zeros((Ny, Nx))
   current_deformation_diff = np.zeros((Ny, Nx))
 
   print(f"Precomputing parts for the filtering.".center(column_size))
@@ -230,14 +230,14 @@ def use_kajiura_filter(deformation, bathymetry, spatial_resolution, filtering_de
     print(f"Filter is in timestep {t+1:{frmt}d} of {Ntime} with a filtering depth of {kajiura_depth} m.".center(column_size))
 
     # Start filterting
-    current_η[:] = 0.0
+    current_η_diff[:] = 0.0
     current_deformation_diff[:] = current_deformation - current_disp  
-    apply_kajiura_fft(current_bathymetry, current_deformation_diff, current_η, 
+    apply_kajiura_fft(current_bathymetry, current_deformation_diff, current_η_diff, 
                                -np.min(bathymetry), spatial_resolution, spatial_resolution, 
                                precalc_σ, precalc_R, kajiura_depth)
     current_disp = deformation[t]
 
-    filtered_deformation[t] = current_η
+    filtered_deformation[t] = filtered_deformation[max(0, t-1)] + current_η_diff
 
     stop = time.time()
     print(f"Timestep {t+1} took {stop - start} seconds.".center(column_size))

exchangeGrid/interpolateBathy.py

@@ -39,7 +39,10 @@ def get_interpolatedBathy(bathy_file, donor_x, donor_y, donor_deformation):
   except:
     bathy_x = bathy_data.variables['lon']
     bathy_y = bathy_data.variables['lat']
-    bathy = bathy_data.variables['elevation']
+    try:
+      bathy = bathy_data.variables['elevation']
+    except:
+      bathy = bathy_data.variables['z']
 
   # Create interpolation class 
   # Note that x and y-coordinates are flipped due to the storage in the netCDF file

interface_module.py

@@ -20,11 +20,11 @@
 This script can be run from the command line and needs some additional arguments.
 
 How to run the script: 
-  python interface_module.py --donor donor_model --CRS_reference lon, lat donor_output bathy_file (--resolution resolution --only_donor_domain --receiver receiver_model --filter filter --filtering_depth filtering_depth --casename casename --include_horizontal_deformation)
+  python interface_module.py --donor donor_model --projection proj4string donor_output bathy_file (--resolution resolution --only_donor_domain --receiver receiver_model --filter filter --filtering_depth filtering_depth --casename casename --include_horizontal_deformation)
 
 Arguments that need/can to be provided:
   * --donor donor_model           where donor_model = seissol, shaltop, bingclaw (all lower case!) 
-  * --CRS_reference                    CRS coordinates reference (lon, lat of lower left corner of domain)
+  * --projection proj4_string     parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates
   * donor_output                       name of the output file or path from the donor model
   * bathy_file                             name of the bathymetry file. Domain has to be larger compared to the domain from the donor model
   * --resolution resolution:          (optional) spatial resolution the donor output will be interpolated to (will be used for both x- and y-coordinates; has to be provided in meters)
@@ -50,7 +50,6 @@
 import os
 
 #TODO: include functionality for parameter file ?
-#TODO: Check if CRS_reference can be removed after SeisSol and SHALTOP have a CRS reference; In this case, the CRS_reference in each donor module needs to be changed
 
 # Some definitions for a nice print on the terminal
 column_size = os.get_terminal_size().columns
@@ -62,12 +61,10 @@
 parser.add_argument("-d", "--donor",
     help="donor model; seissol, shaltop, bingclaw (all lower case)",
     required=True,)
-parser.add_argument("-crs", "--CRS_reference", 
-    metavar=("[lon, lat]"),
-    nargs=2,
+parser.add_argument("-p", "--projection",
     type=str,
-    help="CRS reference coordinates (list of longitude and latitude of lower left corner of the domain); required for SHALTOP and SeisSol.",
-    default = ["0", "0"])
+    help="projection parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates",
+    default = "epsg:4326")
 parser.add_argument("donor_output", help="name of the output file(s) from the donor model")
 parser.add_argument("bathy_file", help="name of the bathymetry file")
 parser.add_argument(
@@ -99,7 +96,7 @@
 filtername = args.filter
 filtering_depth = float(args.filtering_depth)
 casename = args.casename
-CRS_reference = args.CRS_reference
+projection = args.projection
 
 print(asterisk_fill + "\n")
 print("WP6 Interface module\n".center(column_size))
@@ -115,7 +112,7 @@
 
 start = time.time()
 
-donor_deformation, donor_x, donor_y, donor_time, donor_bathy, eg_resolution = donorInterface.get_donorModel(args.donor, args.donor_output, spatial_resolution, CRS_reference, args.bathy_file, incl_horizontal)
+donor_deformation, donor_x, donor_y, donor_time, donor_bathy, eg_resolution = donorInterface.get_donorModel(args.donor, args.donor_output, spatial_resolution, projection, args.bathy_file, incl_horizontal)
 
 stop = time.time()
 print((f"Stage 1 completed. It took {stop - start} seconds.\n").center(column_size))
@@ -150,8 +147,8 @@
 
 start = time.time()
 
-if (spatial_resolution > 0.0):
-  writeBathy.write_interpolatedBathy(args.receiver, eg_bathymetry, eg_x, eg_y, casename,  Ntime=np.shape(eg_deformation)[0])
+# if (spatial_resolution > 0.0):
+writeBathy.write_interpolatedBathy(args.receiver, eg_bathymetry, eg_x, eg_y, casename,  Ntime=np.shape(eg_deformation)[0])
 writeUplift.write_deformation(eg_deformation, eg_x, eg_y, donor_time, args.receiver, args.donor, filtername, casename, spatial_resolution)
 
 stop = time.time()

receiverModel/deformation2hysea.py

@@ -47,7 +47,7 @@ def write2hysea(eg_deformation, eg_x, eg_y, eg_time, donor, filtername, casename
   time_dim = ds.createDimension('time', None)
 
   time = ds.createVariable('time', 'f4', ('time',))
-  time.units = 'time step'
+  time.units = 'seconds'
   latitude = ds.createVariable('y', 'f8', ('y',))
   latitude.units = 'degrees north (WGS84)'
   latitude.long_name = 'latitude'