diff --git a/README.md b/README.md index c22284c..f4e4d54 100644 --- a/README.md +++ b/README.md @@ -19,13 +19,13 @@ Latest changes made in 04/24 by M. Bänsch (UHAM) 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_output bathy_file --donor donor_model --projection projection_donor (--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) * 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 + * --donor donor_model where donor_model = seissol, shaltop, bingclaw (all lower case!) + * --projection projection_donor projection parameters (in Proj4 format) for converting from Cartesian to geographic (lon, lat) coordinates * --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) * --only_donor_domain (optional) handle to only use the domain given by the donor model (False by default) * --receiver receiver_model (optional) receiver model (as of now, only hysea is available) diff --git a/donorModels/donorInterface.py b/donorModels/donorInterface.py index 56808ce..87c61ba 100644 --- a/donorModels/donorInterface.py +++ b/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") diff --git a/donorModels/donorResolution.py b/donorModels/donorResolution.py index f16fb0d..c86106b 100644 --- a/donorModels/donorResolution.py +++ b/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): diff --git a/donorModels/seissol2exchangegrid.py b/donorModels/seissol2exchangegrid.py index 5f78da0..5aff5a5 100644 --- a/donorModels/seissol2exchangegrid.py +++ b/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 diff --git a/donorModels/shaltop2exchangegrid.py b/donorModels/shaltop2exchangegrid.py index 5113b3e..0e236d9 100644 --- a/donorModels/shaltop2exchangegrid.py +++ b/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)) diff --git a/exchangeGrid/exchangeGridCreation.py b/exchangeGrid/exchangeGridCreation.py index 45c31c7..ca30234 100644 --- a/exchangeGrid/exchangeGridCreation.py +++ b/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) diff --git a/exchangeGrid/filter.py b/exchangeGrid/filter.py index fa37b9c..d0d28e9 100644 --- a/exchangeGrid/filter.py +++ b/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)) diff --git a/exchangeGrid/interpolateBathy.py b/exchangeGrid/interpolateBathy.py index f25c31a..e6fd65a 100644 --- a/exchangeGrid/interpolateBathy.py +++ b/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 diff --git a/interface_module.py b/interface_module.py index dfd898b..f553af7 100644 --- a/interface_module.py +++ b/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() diff --git a/receiverModel/deformation2hysea.py b/receiverModel/deformation2hysea.py index 3656528..524d59e 100644 --- a/receiverModel/deformation2hysea.py +++ b/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'