From 8d084b199fb7b0c70e632f255b874e987660c553 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthieu=20S=C3=A9cher?= <matthieu.secher@edf.fr> Date: Wed, 10 May 2023 16:22:26 +0200 Subject: [PATCH] [scripts] Add features to MascaretGeoFile class Add reading of topo-bath information which is given by a string in the geometry file for each points of each profile ("B" stand for bathymetry and "T" for topography). This information is used to mark the limits between minor and major bed. Add reading Courlis geometry files (geoC and georefC) Add an option to the class to activate or desactivate loading of the file + some syntax corrections --- .../data_manip/formats/mascaretgeo_file.py | 196 +++++++++++++----- 1 file changed, 143 insertions(+), 53 deletions(-) diff --git a/scripts/python3/data_manip/formats/mascaretgeo_file.py b/scripts/python3/data_manip/formats/mascaretgeo_file.py index 8d9582f106..249b1228e1 100644 --- a/scripts/python3/data_manip/formats/mascaretgeo_file.py +++ b/scripts/python3/data_manip/formats/mascaretgeo_file.py @@ -5,7 +5,8 @@ from data_manip.formats.mascaret_file import Reach, Section from utils.exceptions import MascaretException -class MascaretGeoFile: + +class MascaretGeoFile(): """ Parse Mascaret geometry file (geo/geoC/georef/georefC) Handles multiple reaches @@ -23,9 +24,10 @@ class MascaretGeoFile: """ OUTPUT_FLOAT_FMT = '%.6f' - def __init__(self, file_name, fformat=None, mode='read'): + def __init__(self, file_name, load_file=True, fformat=None, mode='read'): """ @param file_name (str) file name + @param load_file (bool) option to activate loading of the file @param fformat (str) file format ('opt' or 'rub') @param mode (str) define the mode for the class, 'read' by default to read a file, @@ -44,18 +46,19 @@ class MascaretGeoFile: self.fformat = os.path.splitext(file_name)[1][1:] else: self.fformat = fformat.lower().strip() - if self.fformat not in ('geo', 'georef'): + if self.fformat not in ('geo', 'geoC', 'georef', 'georefC'): raise NotImplementedError( 'Format `%s` not supported,\ - only geo and georef formats are supported as input' % - self.fformat) + only geo, geoC, georef and georefC formats are supported\ + as input' % self.fformat) self.has_ref = 'ref' in self.fformat # Layers for sediments (Courlis) self.has_layers = self.fformat.endswith('C') - # Load file content - self.load() + if load_file: + # Load file content + self.load() def load(self): """ @@ -67,8 +70,9 @@ class MascaretGeoFile: section_id = 0 section_name = '' section_pk = -1.0 - dist, x_list, y_list, z_list = [], [], [], [] - xa, ya = None, None + dist, x_list, y_list, z_list, topo_bath_list, layers_elev_list = \ + [], [], [], [], [], [] + x_axis, y_axis = None, None for line in filein: if line.startswith('#'): @@ -78,19 +82,25 @@ class MascaretGeoFile: if dist: # Add previous Section section = Section(section_id, section_pk, section_name) - if xa is not None and ya is not None: - section.set_axis(xa, ya) + if x_axis is not None and y_axis is not None: + section.set_axis(x_axis, y_axis) if self.has_ref: - section.set_points_from_xyz(x_list, y_list, z_list) + section.set_points_from_xyz(x_list, y_list, z_list, + topo_bath_list) else: - section.set_points_from_trans(dist, z_list) + section.set_points_from_trans(dist, z_list, + topo_bath_list) + section.distances = dist + if self.has_layers: + section.add_layers_from_elevations( + layers_elev_list) reach.add_section(section) if self.has_ref: _, reach_name, section_name, pk_str, x1, y1, x2, y2,\ - _, xa, ya = line.split() - xa = float(xa) - ya = float(ya) + _, x_axis, y_axis = line.split() + x_axis = float(x_axis) + y_axis = float(y_axis) else: _, reach_name, section_name, pk_str = line.split() @@ -108,28 +118,89 @@ class MascaretGeoFile: # Reset variables to store section section_pk = float(pk_str) - dist, x_list, y_list, z_list = [], [], [], [] + dist, x_list, y_list, z_list, topo_bath_list, \ + layers_elev_list = [], [], [], [], [], [] section_id += 1 else: - if self.has_ref: - dist_str, z_str, _, x1, y1 = line.split() - x_list.append(float(x1)) - y_list.append(float(y1)) + list_line = line.split() + + if self.has_ref and self.has_layers is False: + dist_str = list_line[0] + z_str = list_line[1] + if "T" in list_line or "B" in list_line: + shift = 0 + topo_bath = list_line[2] + else: + shift = -1 + topo_bath = "B" + + x_list.append(float(list_line[shift + 3])) + y_list.append(float(list_line[shift + 4])) + + elif self.has_ref and self.has_layers: + dist_str = list_line[0] + z_str = list_line[1] + if "T" in list_line or "B" in list_line: + shift = 0 + topo_bath = list_line[-3] + else: + shift = 1 + topo_bath = "B" + + layers_elev = list_line[2:shift - 3] + + x_list.append(float(list_line[-2])) + y_list.append(float(list_line[-1])) + + layers_elev = [float(elev) for elev in layers_elev] + + if self.nlayers == 0: + self.nlayers = len(layers_elev) + + elif self.has_ref is False and self.has_layers: + dist_str = list_line[0] + z_str = list_line[1] + if "T" in list_line or "B" in list_line: + layers_elev = list_line[2:-1] + topo_bath = list_line[-1] + else: + layers_elev = list_line[2:] + topo_bath = "B" + + layers_elev = [float(elev) for elev in layers_elev] + + if self.nlayers == 0: + self.nlayers = len(layers_elev) + else: - dist_str, z_str, _ = line.split() + dist_str = list_line[0] + z_str = list_line[1] + + if "T" in list_line or "B" in list_line: + topo_bath = list_line[2] + else: + topo_bath = "B" # Add new point to current section dist.append(float(dist_str)) z_list.append(float(z_str)) + topo_bath_list.append(topo_bath) + if self.has_layers: + layers_elev_list.append(layers_elev) # Add last section section = Section(section_id, section_pk, section_name) - if xa is not None and ya is not None: - section.set_axis(xa, ya) + if x_axis is not None and y_axis is not None: + section.set_axis(x_axis, y_axis) if self.has_ref: - section.set_points_from_xyz(x_list, y_list, z_list) + section.set_points_from_xyz(x_list, y_list, z_list, + topo_bath_list) else: - section.set_points_from_trans(dist, z_list) + section.set_points_from_trans(dist, z_list, topo_bath_list) + section.distances = dist + if self.has_layers: + section.add_layers_from_elevations(layers_elev_list) + self.layer_names = section.layer_names reach.add_section(section) def save(self, output_file_name): @@ -149,60 +220,79 @@ class MascaretGeoFile: if ref and not self.has_ref: raise MascaretException('Could not write `%s` format without\ - any geo-referenced data' % fformat) + any geo-referenced data' % fformat) - with open(output_file_name, 'w') as fileout: + with open(output_file_name, 'w', encoding="utf-8") as fileout: for _, reach in self.reaches.items(): for sec in reach: positions_str = '' if ref: # Get river_banks and `AXE` coordinates if necessary - xa, ya = sec.axis - positions_str += ' %f %f %f %f' %\ + x_axis, y_axis = sec.axis + positions_str += '%f %f %f %f' %\ (sec.x[0], sec.y[0], sec.x[-1], sec.y[-1]) - positions_str += ' AXE %f %f' % (xa, ya) + positions_str += ' AXE %f %f' % (x_axis, y_axis) # Write profile header fileout.write( - 'Profil %s %s %f%s\n' % + 'Profil %s %s %f %s\n' % (reach.name, sec.name, sec.pk, positions_str)) # Write points and layers if necessary if not ref and not layers: - for dist, z in zip(sec.distances, sec.z): - fileout.write('%f %f B\n' % (dist, z)) + for dist, z, topo_bath in zip(sec.distances, sec.z, + sec.topo_bath): + fileout.write('%f %f %s\n' % (dist, z, topo_bath)) elif ref and not layers: - for dist, x, y, z in zip(sec.distances, - sec.x, sec.y, sec.z): - fileout.write('%f %f B %f %f\n' % (dist, z, x, y)) + for dist, x, y, z, topo_bath in zip(sec.distances, + sec.x, + sec.y, + sec.z, + sec.topo_bath): + fileout.write('%f %f %s %f %f\n' + % (dist, z, topo_bath, x, y)) elif not ref and layers: - for i, (dist, z) in enumerate(zip(sec.distances, - sec.z)): + for i, (dist, z, topo_bath) in \ + enumerate(zip(sec.distances, sec.z, + sec.topo_bath)): if self.nlayers == 0: layers_str = '' else: layers_str = ' ' +\ - ' '.join( - [MascaretGeoFile.OUTPUT_FLOAT_FMT % zl - for zl in sec.layers_elev[:, i]]) - fileout.write('%f %f%s B\n' % - (dist, z, layers_str)) + ' '.join( + [MascaretGeoFile.OUTPUT_FLOAT_FMT % zl + for zl in sec.layers_elev[:, i]]) + + try: + fileout.write('%f %f %s %s\n' % + (dist, z, layers_str, topo_bath)) + except Exception as e: + raise e elif ref and layers: - for i, (dist, x, y, z) in enumerate(zip(sec.distances, - sec.x, sec.y, - sec.z)): + for i, (dist, x, y, z, topo_bat) in enumerate(zip(sec.distances, + sec.x, + sec.y, + sec.z, + sec.topo_bath)): + if self.nlayers == 0: layers_str = '' else: layers_str = ' ' + ' '\ - .join([MascaretGeoFile.OUTPUT_FLOAT_FMT % - zl for zl in sec.layers_elev[:, i]]) - fileout.write('%f %f%s B %f %f\n' % - (dist, z, layers_str, x, y)) + .join([MascaretGeoFile.OUTPUT_FLOAT_FMT % + zl for zl in sec.layers_elev[:, i]]) + + try: + fileout.write('%f %f %s %s %f %f\n' + % (dist, z, layers_str, + topo_bath, x, y)) + except Exception as e: + raise e + def __repr__(self): return 'MascaretGeoFile: %s' % self.file_name @@ -241,8 +331,8 @@ if __name__ == '__main__': 'examples', 'mascaret'), '*.geo') geo_files += recursive_glob(os.path.join(os.environ['HOMETEL'], 'examples', 'mascaret'), 'geometrie') - for file_name in sorted(geo_files): - geo_file = MascaretGeoFile(file_name, 'geo') + for filename in sorted(geo_files): + geo_file = MascaretGeoFile(filename, 'geo') print(geo_file.summary()) except MascaretException as e: print(str(e)) -- GitLab