From 9d6fb8d0bb6c956215d6a409cdf14c176fbd424d Mon Sep 17 00:00:00 2001 From: Bastiaan Roos Date: Thu, 2 May 2019 19:44:33 +0200 Subject: [PATCH 1/2] fixes --- qt_models/area_tree.py | 4 +- qt_models/legger_tree.py | 2 +- sql_models/legger.py | 62 ++++++++++++++----------- sql_models/legger_views.py | 23 +++++----- utils/geom_collections/base.py | 19 +++++--- utils/legger_map_manager.py | 2 +- utils/map_layers.py | 4 +- utils/read_tdi_results.py | 82 +++++++++++++++++++++++++++++++--- utils/theoretical_profiles.py | 14 +++--- views/calculating_profiles.py | 69 +++++++++++++--------------- views/legger_network_widget.py | 5 +++ views/network_graph_widgets.py | 2 +- views/polder_selection.py | 57 +++++++++++++++-------- 13 files changed, 227 insertions(+), 118 deletions(-) diff --git a/qt_models/area_tree.py b/qt_models/area_tree.py index 19a9a18..b345ecf 100644 --- a/qt_models/area_tree.py +++ b/qt_models/area_tree.py @@ -10,8 +10,8 @@ {'field': 'selected', 'show': False, 'field_type': CHECKBOX_FIELD, 'column_width': 50, 'single_selection': True}, {'field': 'hover', 'show': False, 'field_type': CHECKBOX_FIELD, 'column_width': 50}, - {'field': 'weight', 'header': 'omvang', 'column_width': 50}, - {'field': 'distance', 'header': 'afstand', 'column_width': 50}, + {'field': 'weight', 'header': 'prioriteit', 'column_width': 50}, + {'field': 'distance', 'header': 'afstand', 'column_width': 50, 'show': False}, ) diff --git a/qt_models/legger_tree.py b/qt_models/legger_tree.py index 2133445..ce41e40 100644 --- a/qt_models/legger_tree.py +++ b/qt_models/legger_tree.py @@ -76,7 +76,7 @@ def __init__(self, data_dict, feature): 'selected_variant_id': 'geselecteerde_variant', 'begroeiingsvariant_id': 'begroeiingsvariant_id', 'selected_begroeiingsvariant_id': 'geselecteerde_begroeiingsvariant', - 'selected_remarks': 'selectie_opmerkingen', + 'selected_remarks': 'opmerkingen', } # set default values diff --git a/sql_models/legger.py b/sql_models/legger.py index 1ee6d03..0dc4e61 100644 --- a/sql_models/legger.py +++ b/sql_models/legger.py @@ -27,7 +27,7 @@ def get_or_create(session, model, defaults=None, **kwargs): class BegroeiingsVariant(Base): __tablename__ = 'begroeiingsvariant' - id = Column(Integer(), primary_key=True, autoincrement=True) + id = Column(Integer(), primary_key=True, autoincrement=True, index=True) naam = Column(String(20)) friction = Column(Float()) is_default = Column(Boolean, default=False) @@ -48,7 +48,7 @@ class Waterdeel(Base): __tablename__ = 'waterdeel' extend_existing = True - id = Column(Integer, primary_key=True) + id = Column(Integer, primary_key=True, index=True) shape_length = Column(Float) shape_area = Column(Float) geometry = Column("GEOMETRY", Geometry(geometry_type='MULTIPOLYGON', srid=28992)) @@ -62,7 +62,7 @@ class HydroObject(Base): __tablename__ = 'hydroobject' extend_existing = True - id = Column(Integer, primary_key=True) + id = Column(Integer, primary_key=True, index=True) geometry = Column("GEOMETRY", Geometry(geometry_type='MULTILINESTRING', srid=28992)) code = Column(String(50), index=True) categorieoppwaterlichaam = Column(Integer) @@ -72,8 +72,9 @@ class HydroObject(Base): flowline_id = Column(Integer) # link to 3di id score = Column(Float) begroeiingsvariant_id = Column(Integer, - ForeignKey(BegroeiingsVariant.__tablename__ + ".id")) - + ForeignKey(BegroeiingsVariant.__tablename__ + ".id"), + index=True) + opmerkingen = Column(String()) # shape_length = Column(Float) profielen = relationship("Profielen", @@ -103,7 +104,6 @@ class HydroObject(Base): # primaryjoin="Varianten.begroeiingsvariant_id == BegroeiingsVariant.id", uselist=False) - def __str__(self): return u'Hydro object {0}'.format( self.code) @@ -112,12 +112,13 @@ def __str__(self): class Profielen(Base): __tablename__ = 'profielen' - id = Column(Integer, primary_key=True) + id = Column(Integer, primary_key=True, index=True) proident = Column(String(24)) bron_profiel = Column(String(50)) pro_id = Column(Integer, index=True) hydro_id = Column(Integer, - ForeignKey(HydroObject.__tablename__ + ".id")) + ForeignKey(HydroObject.__tablename__ + ".id"), + index=True) # shape_lengte = Column(Float) hydro = relationship(HydroObject, back_populates="profielen") @@ -134,7 +135,7 @@ def __str__(self): class Profielpunten(Base): __tablename__ = 'profielpunten' - objectid = Column(Integer, primary_key=True, autoincrement=True) + objectid = Column(Integer, primary_key=True, autoincrement=True, index=True) pbp_id = Column(Integer) prw_id = Column(Integer) pbpident = Column(String(24)) @@ -143,7 +144,8 @@ class Profielpunten(Base): iws_hoogte = Column(Float) afstand = Column(Float) pro_pro_id = Column(Integer, - ForeignKey(Profielen.__tablename__ + '.pro_id')) + ForeignKey(Profielen.__tablename__ + '.pro_id'), + index=True) geometry = Column("GEOMETRY", Geometry(geometry_type='POINT', srid=28992)) profiel = relationship( @@ -158,7 +160,7 @@ def __str__(self): class Kenmerken(Base): __tablename__ = 'kenmerken' - id = Column(Integer, primary_key=True) + id = Column(Integer, primary_key=True, index=True) diepte = Column(Float) bron_diepte = Column(String(50)) bodemhoogte = Column(Float) @@ -170,8 +172,8 @@ class Kenmerken(Base): grondsoort = Column(String(50)) bron_grondsoort = Column(String(50)) hydro_id = Column(Integer, - ForeignKey(HydroObject.__tablename__ + ".id")) - + ForeignKey(HydroObject.__tablename__ + ".id"), + index=True) na_lengte = Column(Float) voor_lengte = Column(Float) new_category = Column(Integer) @@ -188,10 +190,12 @@ def __str__(self): class Varianten(Base): __tablename__ = 'varianten' - id = Column(String(), primary_key=True) + id = Column(String(), primary_key=True, index=True) begroeiingsvariant_id = Column(Integer, - ForeignKey(BegroeiingsVariant.__tablename__ + ".id")) - diepte = Column(Float) + ForeignKey(BegroeiingsVariant.__tablename__ + ".id"), + index=True) + diepte = Column(Float, + index=True) waterbreedte = Column(Float) bodembreedte = Column(Float) talud = Column(Float) @@ -199,7 +203,8 @@ class Varianten(Base): verhang_bos_bijkerk = Column(Float) opmerkingen = Column(String()) hydro_id = Column(Integer, - ForeignKey(HydroObject.__tablename__ + ".id")) + ForeignKey(HydroObject.__tablename__ + ".id"), + index=True) begroeiingsvariant = relationship(BegroeiingsVariant, # foreign_keys='begroeiingsvariant_id', @@ -230,11 +235,12 @@ class GeselecteerdeProfielen(Base): hydro_id = Column(Integer, ForeignKey(HydroObject.__tablename__ + ".id"), - primary_key=True) + primary_key=True, + index=True) variant_id = Column(String(), - ForeignKey(Varianten.__tablename__ + ".id")) + ForeignKey(Varianten.__tablename__ + ".id"), + index=True) selected_on = Column(DateTime, default=datetime.datetime.utcnow) - opmerkingen = Column(String()) hydro = relationship(HydroObject, back_populates="geselecteerd") @@ -248,8 +254,12 @@ class ProfielFiguren(Base): # object_id = Column(Integer, primary_key=True) hydro_id = Column('id_hydro', Integer, - ForeignKey(HydroObject.__tablename__ + ".id")) - profid = Column(String(16), ForeignKey(Varianten.__tablename__ + ".id"), primary_key=True, ) + ForeignKey(HydroObject.__tablename__ + ".id"), + index=True) + profid = Column(String(16), + ForeignKey(Varianten.__tablename__ + ".id"), + primary_key=True, + index=True) type_prof = Column(String(1)) coord = Column(String()) peil = Column(Float) @@ -282,7 +292,7 @@ class DuikerSifonHevel(Base): __tablename__ = 'duikersifonhevel' extend_existing = True - id = Column(Integer, primary_key=True) + id = Column(Integer, primary_key=True, index=True) geometry = Column("GEOMETRY", Geometry(geometry_type='MULTILINESTRING', srid=28992)) code = Column(String(50), index=True) categorie = Column(Integer) @@ -294,9 +304,9 @@ class DuikerSifonHevel(Base): vormkoker = Column(Float) # shape_lengte = Column(Float) - debiet = Column(Float) # extra? - channel_id = Column(Integer) # extra? - flowline_id = Column(Integer) # extra? + debiet = Column(Float) # extra + channel_id = Column(Integer) # extra + flowline_id = Column(Integer) # extra def __str__(self): return u'DuikerSifonHevel {0}'.format( diff --git a/sql_models/legger_views.py b/sql_models/legger_views.py index 9e47674..a4f892b 100644 --- a/sql_models/legger_views.py +++ b/sql_models/legger_views.py @@ -1,13 +1,13 @@ def create_legger_views(session): session.execute( """ - DROP VIEW IF EXISTS hydroobjects_kenmerken4; + DROP VIEW IF EXISTS hydroobjects_kenmerken; """ ) session.execute( """ - CREATE VIEW hydroobjects_kenmerken4 AS + CREATE VIEW hydroobjects_kenmerken AS SELECT h.id, code, @@ -27,7 +27,7 @@ def create_legger_views(session): geselecteerd_breedte, geselecteerde_variant, geselecteerde_begroeiingsvariant, - selectie_opmerkingen, + h.opmerkingen, CASE WHEN h.debiet >= 0 THEN "GEOMETRY" WHEN h.debiet THEN ST_REVERSE("GEOMETRY") @@ -64,8 +64,7 @@ def create_legger_views(session): v.diepte as geselecteerd_diepte, v.waterbreedte as geselecteerd_breedte, v.id as geselecteerde_variant, - v.begroeiingsvariant_id as geselecteerde_begroeiingsvariant, - g.opmerkingen as selectie_opmerkingen + v.begroeiingsvariant_id as geselecteerde_begroeiingsvariant FROM geselecteerd g, varianten v WHERE g.variant_id = v.id) as sel ON sel.hydro_id = h.id @@ -73,7 +72,7 @@ def create_legger_views(session): session.execute( """ - DELETE FROM views_geometry_columns WHERE view_name = 'hydroobjects_kenmerken3'; + DELETE FROM views_geometry_columns WHERE view_name = 'hydroobjects_kenmerken'; """ ) @@ -81,14 +80,14 @@ def create_legger_views(session): """ INSERT INTO views_geometry_columns (view_name, view_geometry, view_rowid, f_table_name, f_geometry_column, read_only) - VALUES('hydroobjects_kenmerken4', 'geometry', 'id', 'hydroobject', 'geometry', 1); + VALUES('hydroobjects_kenmerken', 'geometry', 'id', 'hydroobject', 'geometry', 1); """) session.execute( """ INSERT INTO views_geometry_columns (view_name, view_geometry, view_rowid, f_table_name, f_geometry_column, read_only) - VALUES('hydroobjects_kenmerken4', 'line', 'id', 'hydroobject', 'geometry', 1); + VALUES('hydroobjects_kenmerken', 'line', 'id', 'hydroobject', 'geometry', 1); """) session.commit() @@ -151,13 +150,13 @@ def create_legger_views(session): session.execute( """ - DROP VIEW IF EXISTS begroeiingsadvies4; + DROP VIEW IF EXISTS begroeiingsadvies; """ ) session.execute( """ - CREATE VIEW begroeiingsadvies4 AS + CREATE VIEW begroeiingsadvies AS SELECT bv.naam as advies_naam, a.* @@ -245,7 +244,7 @@ def create_legger_views(session): session.execute( """ - DELETE FROM views_geometry_columns WHERE view_name = 'begroeiingsadvies4'; + DELETE FROM views_geometry_columns WHERE view_name = 'begroeiingsadvies'; """ ) @@ -253,7 +252,7 @@ def create_legger_views(session): """ INSERT INTO views_geometry_columns (view_name, view_geometry, view_rowid, f_table_name, f_geometry_column, read_only) - VALUES('begroeiingsadvies4', 'geometry', 'id', 'hydroobject', 'geometry', 1); + VALUES('begroeiingsadvies', 'geometry', 'id', 'hydroobject', 'geometry', 1); """) diff --git a/utils/geom_collections/base.py b/utils/geom_collections/base.py index d38eced..7ac2edc 100644 --- a/utils/geom_collections/base.py +++ b/utils/geom_collections/base.py @@ -4,6 +4,7 @@ from collections import OrderedDict + # use fiona collection for 'normal' use @@ -101,7 +102,6 @@ def keys(self, start=0, stop=None, step=1, **kwds): qbbox = QgsRectangle(*bbox) selected.intersection_update(set(self._spatial_index.intersects(qbbox))) - if mask: # todo pass @@ -129,16 +129,25 @@ def writerecords(self, records): for record in records: record['id'] = nr self.ordered_dict[nr] = record - geom = shape(record['geometry']) + # rtree # self._spatial_index.insert(nr, geom) # QGIS: feature = QgsFeature() feature.setFeatureId(nr) - qgeom = QgsGeometry() - qgeom.fromWkb(geom.to_wkb()) - feature.setGeometry(qgeom) + try: + geom = shape(record['geometry']) + qgeom = QgsGeometry() + qgeom.fromWkb(geom.to_wkb()) + + except: + wkt = "{}({})".format( + record['geometry']['type'], + ",".join(["{} {}".format(*c) for c in record['geometry']['coordinates']])) + qgeom = QgsGeometry() + qgeom.fromWkt(wkt) + feature.setGeometry(qgeom) self._spatial_index.insertFeature(feature) nr += 1 diff --git a/utils/legger_map_manager.py b/utils/legger_map_manager.py index cd22174..bab0113 100644 --- a/utils/legger_map_manager.py +++ b/utils/legger_map_manager.py @@ -52,7 +52,7 @@ def get_layer(spatialite_path, table_name, geom_column=''): layer = get_layer( self.path_legger_db, - 'hydroobjects_kenmerken4', + 'hydroobjects_kenmerken', geometry_col) # todo: remove this filter when bidirectional islands are supported diff --git a/utils/map_layers.py b/utils/map_layers.py index 65acbe2..03b5a18 100644 --- a/utils/map_layers.py +++ b/utils/map_layers.py @@ -66,8 +66,8 @@ def add_layers_to_map(self): ]), ('tbv begroeiingsgraad', [ ('aanwijzen', 'hydroobject', 'begroeiingsvariant_id', 'begroeiingsvariant', 'geometry', None), - ('begroeiingsadvies', 'begroeiingsadvies4', 'advies_id', 'begroeiingsvariant', 'geometry', None), - ('begroeiingsvariant', 'begroeiingsadvies4', 'aangew_bv_id', 'begroeiingsvariant', + ('begroeiingsadvies', 'begroeiingsadvies', 'advies_id', 'begroeiingsvariant', 'geometry', None), + ('begroeiingsvariant', 'begroeiingsadvies', 'aangew_bv_id', 'begroeiingsvariant', 'geometry', None), # ('sterk min profiel', 'ruimte_view', 'ruim', 'min_max_line', 'geometry', None), # ('ruimte', 'ruimte_view', 'over_width', 'min_max_line', 'geometry', None), diff --git a/utils/read_tdi_results.py b/utils/read_tdi_results.py index 270abd5..7349cb6 100644 --- a/utils/read_tdi_results.py +++ b/utils/read_tdi_results.py @@ -10,11 +10,11 @@ from qgis.core import QgsGeometry, QgsLineStringV2, QgsPoint try: - from ThreeDiToolbox.datasource.netcdf_groundwater_h5py import NetcdfGroundwaterDataSourceH5py as NetcdfGroundwaterDataSource + from ThreeDiToolbox.datasource.netcdf_groundwater_h5py import \ + NetcdfGroundwaterDataSourceH5py as NetcdfGroundwaterDataSource except ImportError: from ThreeDiToolbox.datasource.netcdf_groundwater import NetcdfGroundwaterDataSource - from legger.utils.geom_collections.lines import LineCollection from legger.utils.geometries import LineString from legger.utils.geometries import shape @@ -22,7 +22,6 @@ from legger.sql_models.legger_database import LeggerDatabase from pyspatialite import dbapi2 as dbapi - log = logging.getLogger(__name__) @@ -142,6 +141,26 @@ def read_tdi_results(path_model_db, path_result_db, } }) + # if channel does not match, flowlines are used (relevant for culverts, weir.) + flowline_col = LineCollection() + + flowline_cursor = con_res.execute( + "SELECT " + "id, " + "ASGEOJSON(TRANSFORM(the_geom, 28992)) as geojson " + "FROM flowlines " + "WHERE type in ('v2_channel', 'v2_culvert', 'v2_orifice', 'v2_weir');" + ) + + for flowline in flowline_cursor.fetchall(): + flowline_col.write({ + 'geometry': json.loads(flowline['geojson']), + 'properties': { + 'id': flowline['id'], + 'q_end': qts[flowline['id'] - 1] + } + }) + hydro_cursor = con_legger.execute( 'SELECT ' 'id, ' @@ -154,7 +173,7 @@ def read_tdi_results(path_model_db, path_result_db, # loop over all hydroobjects and find flowline through link with channel for hydroobject in hydro_cursor.fetchall(): - if hydroobject['id'] in [198616, 198362, 659550]: + if hydroobject[0] in [157078]: a = 1 line = create_geom_line(json.loads(hydroobject['geojson'])['coordinates']) @@ -182,7 +201,7 @@ def read_tdi_results(path_model_db, path_result_db, points_on_line = [ line.interpolate((float(nr) / 10.0) * length) # normalized=True for nr in range(1, 10) - ] + ] # print(', '.join([str(point.x) for point in points_on_line])) @@ -274,6 +293,59 @@ def read_tdi_results(path_model_db, path_result_db, 'nr_candidates': len(candidates), 'score': score }) + else: + + flowline_candidates = [] + for flowline in flowline_col.filter(bbox=bbox): + geom_flowline = create_geom_line(flowline['geometry']['coordinates']) + distances = [ + geom_flowline.distance(point) + for point in points_on_line + ] + # check how many points are matching + nr_matches = sum([1 for distance in distances if distance <= max_link_distance]) + + if nr_matches >= match_criteria: + # lower score is better + score = sum([ + min(dist, max_link_distance) + for dist in distances + ]) / max_link_distance + flowline_candidates.append(( + score, + flowline, + geom_flowline + )) + + if len(flowline_candidates) > 0: + score, selected_fl, geom_flowline = sorted(flowline_candidates, key=lambda item: item[0])[0] + + # get orientation of hydroobject vs channel + if line.isMultipart(): + vertexes = line.asMultiPolyline() + d1, p1, v1 = geom_flowline.closestSegmentWithContext(vertexes[0][0]) + d2, p2, v2 = geom_flowline.closestSegmentWithContext(vertexes[-1][-1]) + else: + vertexes = line.asPolyline() + d1, p1, v1 = geom_flowline.closestSegmentWithContext(vertexes[0]) + d2, p2, v2 = geom_flowline.closestSegmentWithContext(vertexes[-1]) + + dist1 = geom_flowline.distanceToVertex(v1) - math.sqrt(geom_flowline.sqrDistToVertexAt(p1, v1)) + dist2 = geom_flowline.distanceToVertex(v2) - math.sqrt(geom_flowline.sqrDistToVertexAt(p2, v2)) + + if dist2 < dist1: + factor = -1 + else: + factor = 1 + + hydroobjects.append({ + 'id': hydroobject['id'], + 'qend': selected_fl['properties']['q_end'] * factor, + 'channel_id': None, + 'flowline_id': selected_fl['properties']['id'], + 'nr_candidates': len(flowline_candidates), + 'score': score + }) return hydroobjects diff --git a/utils/theoretical_profiles.py b/utils/theoretical_profiles.py index e436b59..211a01d 100644 --- a/utils/theoretical_profiles.py +++ b/utils/theoretical_profiles.py @@ -143,6 +143,7 @@ def calc_profile_max_ditch_width(object_id, normative_flow, length, slope, max_d in a new dataframe. """ # todo: question. why both bos and bijkerk and manning? + # manning not used for begroeiingsvariants # initial values water_depth = minimal_waterdepth ditch_bottom_width = max_ditch_width - 2 * slope * water_depth @@ -218,7 +219,7 @@ def calc_profile_variants(hydro_objects_satisfy, gradient_norm, friction_bos_bij # todo: waarom round(.., 1)? # First two empty tables: # 1st one with hydro objects that shows the amount of table variants pssoible. - options_table = DataFrame(data=hydro_objects_satisfy.object_id, columns=['object_id', 'possibilities_count']) + # options_table = DataFrame(data=hydro_objects_satisfy.object_id, columns=['object_id', 'possibilities_count']) # 2nd one a table where variants are saved. variants_table = DataFrame(columns=['object_id', 'object_waterdepth_id', 'slope', @@ -309,7 +310,7 @@ def calc_profile_variants(hydro_objects_satisfy, gradient_norm, friction_bos_bij variants_table = variants_table.append(df_temp) count = 1 - options_table.possibilities_count[options_table.object_id == options_table.object_id[i]] = count + # options_table.possibilities_count[options_table.object_id == options_table.object_id[i]] = count # produces a table with the options per hydro object. This is not part of the return, because it is not used. variants_table = variants_table.reset_index(drop=True) @@ -401,14 +402,13 @@ def create_theoretical_profiles(legger_db_filepath, gradient_norm, bv): # Calculate a profile profile = calc_profile_max_ditch_width(object_id, normative_flow, length, slope, max_ditch_width, - gradient_norm, friction_bos_bijkerk=Kb) # todo: Kb + gradient_norm, friction_bos_bijkerk=Kb) # Add the profile to the previous made table where the results are stored profile_max_ditch_width = profile_max_ditch_width.append(profile) - # todo: this function not one indent less? - # When all the results are stored in the table, re-index the table. - profile_max_ditch_width = profile_max_ditch_width.reset_index(drop=True) + # When all the results are stored in the table, re-index the table. + profile_max_ditch_width = profile_max_ditch_width.reset_index(drop=True) log.info("Finished 3: Successfully calculated profiles based on max ditch width\n") @@ -444,7 +444,7 @@ def create_theoretical_profiles(legger_db_filepath, gradient_norm, bv): hydro_objects_unsatisfy = profile_max_ditch_width[profile_max_ditch_width['gradient_bos_bijkerk'] > gradient_norm] - profile_variants = calc_profile_variants(hydro_objects_satisfy, bv.friction) + profile_variants = calc_profile_variants(hydro_objects_satisfy, gradient_norm, bv.friction) log.info("Finished 7: variants created\n") diff --git a/views/calculating_profiles.py b/views/calculating_profiles.py index 1288543..c667742 100644 --- a/views/calculating_profiles.py +++ b/views/calculating_profiles.py @@ -13,7 +13,6 @@ from PyQt4.QtGui import QComboBox, QWidget from legger.sql_models.legger import BegroeiingsVariant, HydroObject, get_or_create from legger.sql_models.legger_database import LeggerDatabase -from legger.sql_models.legger_views import create_legger_views from legger.utils.profile_match_a import doe_profinprof, maaktabellen from legger.utils.read_tdi_results import (get_timestamps, read_tdi_culvert_results, read_tdi_results, write_tdi_culvert_results_to_db, write_tdi_results_to_db) @@ -109,7 +108,7 @@ def __init__( self.last_surge_text = "kies opstuwingsnorm" # fill surge combobox - surge_choices = [self.last_surge_text] + ['%s' % s for s in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]] + surge_choices = ['%s' % s for s in [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]] self.surge_combo_box.insertItems(0, surge_choices) self.surge_combo_box.setCurrentIndex(3) @@ -185,12 +184,14 @@ def execute_step1(self): session = db.get_session() session.query(HydroObject) - timestamp = 'laatst' if self.timestep == 0 else '{0} op {1:.0f} s'.format(self.timestep + 1, - self.timestamps[self.timestep]) + timestamp = 'laatst' if self.timestep == 0 else '{0} op {1:.0f} s'.format( + self.timestep + 1, + self.timestamps[self.timestep]) + self.feedbackmessage = 'Neem tijdstap {0}'.format(timestamp) - # try: - if True: + result = None + try: # read 3di channel results result = read_tdi_results( self.path_model_db, @@ -201,26 +202,25 @@ def execute_step1(self): ) self.feedbackmessage += "\nDatabases zijn gekoppeld." - # except Exception, e: - # self.feedbackmessage += "\nDatabases zijn niet gekoppeld. melding: {0}\n".format(e.message) - # finally: - # self.feedbacktext.setText(self.feedbackmessage) - - try: - # write 3di channel result to legger spatialite - write_tdi_results_to_db( - result, - self.polder_datasource) - - con_legger = dbapi.connect(self.polder_datasource) - create_legger_views(con_legger) - - self.feedbackmessage = self.feedbackmessage + "\n3Di resultaten weggeschreven naar polder database." - except: - self.feedbackmessage = self.feedbackmessage + "\nFout in wegschrijven 3Di resultaten naar polder database." + except Exception, e: + self.feedbackmessage += "\nDatabases zijn niet gekoppeld. melding: {0}\n".format(e.message) finally: self.feedbacktext.setText(self.feedbackmessage) + if result is not None: + try: + # write 3di channel result to legger spatialite + write_tdi_results_to_db( + result, + self.polder_datasource) + + self.feedbackmessage = self.feedbackmessage + "\n3Di resultaten weggeschreven naar polder database." + except: + self.feedbackmessage = self.feedbackmessage + "\nFout in wegschrijven 3Di resultaten naar polder database." + finally: + self.feedbacktext.setText(self.feedbackmessage) + + culv_results = None try: # read 3di culvert results culv_results = read_tdi_culvert_results( @@ -236,15 +236,15 @@ def execute_step1(self): finally: self.feedbacktext.setText(self.feedbackmessage) - try: - # write 3di culvert results to legger spatialite - write_tdi_culvert_results_to_db(culv_results, - self.polder_datasource) - self.feedbackmessage = self.feedbackmessage + "\n3Di culvert resultaten weggeschreven." - except: - self.feedbackmessage = self.feedbackmessage + "\nFout, 3Di culvert resultaten niet weggeschreven." - finally: - self.feedbacktext.setText(self.feedbackmessage) + if culv_results is not None: + try: + # write 3di culvert results to legger spatialite + write_tdi_culvert_results_to_db(culv_results, self.polder_datasource) + self.feedbackmessage = self.feedbackmessage + "\n3Di culvert resultaten weggeschreven." + except: + self.feedbackmessage = self.feedbackmessage + "\nFout, 3Di culvert resultaten niet weggeschreven." + finally: + self.feedbacktext.setText(self.feedbackmessage) def explain_step2(self): """ @@ -279,11 +279,6 @@ def execute_step2(self): # do one query, don't know what the reason was for this... session = db.get_session() - # delete existing variants - # session.execute("Verwijder van varianten") - # session.execute("Verwijder begroeiingsvariant") - session.commit() - get_or_create(session, BegroeiingsVariant, naam='standaard', defaults={'friction': Kb, 'is_default': True}) diff --git a/views/legger_network_widget.py b/views/legger_network_widget.py index 9553c24..135bdc2 100644 --- a/views/legger_network_widget.py +++ b/views/legger_network_widget.py @@ -568,6 +568,8 @@ def data_changed_variant(self, index): if self.legger_model.ep: traject = self.legger_model.ep.up(self.legger_model.selected) traject.reverse() + if len(traject) > 0: + traject.pop(0) else: messagebar_message( 'Traject nodig', @@ -610,6 +612,9 @@ def data_changed_variant(self, index): if self.legger_model.ep: traject = self.legger_model.ep.up(self.legger_model.selected) traject.reverse() + if len(traject) > 0: + traject.pop(0) + else: messagebar_message( 'Traject nodig', diff --git a/views/network_graph_widgets.py b/views/network_graph_widgets.py index 8e6de7e..8b631d7 100644 --- a/views/network_graph_widgets.py +++ b/views/network_graph_widgets.py @@ -554,7 +554,7 @@ def data_changed_legger(self, index): def _get_data(self): if self.legger_model.ep is None: - return [] + return OrderedDict([]) up = self.legger_model.ep.up(end=self.legger_model.sp) out = [] before = up[-1] diff --git a/views/polder_selection.py b/views/polder_selection.py index 79e4466..a96f915 100644 --- a/views/polder_selection.py +++ b/views/polder_selection.py @@ -1,14 +1,16 @@ # -*- coding: utf-8 -*- from __future__ import division -import datetime import logging import os from PyQt4 import QtCore, QtGui from PyQt4.QtCore import QSettings, pyqtSignal from PyQt4.QtGui import QFileDialog, QWidget +from legger.sql_models.legger_views import create_legger_views from legger.utils.read_data_and_make_leggerdatabase import CreateLeggerSpatialite +from legger.utils.user_message import messagebar_message +from pyspatialite import dbapi2 as dbapi log = logging.getLogger(__name__) @@ -20,6 +22,8 @@ def _fromUtf8(s): try: _encoding = QtGui.QApplication.UnicodeUTF8 + + def _translate(context, text, disambig): return QtGui.QApplication.translate(context, text, disambig, _encoding) except AttributeError: @@ -27,7 +31,7 @@ def _translate(context, text, disambig): return QtGui.QApplication.translate(context, text, disambig) -class PolderSelectionWidget(QWidget):#, FORM_CLASS): +class PolderSelectionWidget(QWidget): # , FORM_CLASS): """Dialog for selecting model (spatialite and result files netCDFs)""" closingDialog = pyqtSignal() @@ -51,7 +55,8 @@ def __init__( self.var_text_DAMO.setText("...") self.var_text_HDB.setText("...") - self.var_text_leggerdatabase.setText(self.root_tool.polder_datasource) # De tekst verwijst naar de tekst in de root_tool totdat deze geupdated wordt. + self.var_text_leggerdatabase.setText( + self.root_tool.polder_datasource) # De tekst verwijst naar de tekst in de root_tool totdat deze geupdated wordt. def closeEvent(self, event): """ @@ -66,17 +71,18 @@ def explain_leggerdatabase(self): self.msg_upper_row = QtGui.QMessageBox(self) self.msg_upper_row.setIcon(QtGui.QMessageBox.Information) self.msg_upper_row.setText("Het selecteren van een leggerdatabase") - self.msg_upper_row.setInformativeText("Voor de toewijzing van leggerprofielen wordt een aparte 'leggerdatabase' " - "gemaakt. Deze database is een aparte .sqlite bestand waar data uit " - "DAMO en de Hydrologendatabase (HDB) gecombineerd wordt als randvoorwaarden " - "voor de leggerprofielen, zoals breedte en talud per hydro-object.\n" - "Wanneer een nieuwe leggerdatabase gemaakt moet worden, selecteer dan bij " - "voorkeur de DAMO en HDB die ook voor de opbouw van het 3di model zijn " - "gebruikt.\n" - "Is er al een 'leggerdatabase' aangemaakt, sla deze stap dan over en zorg " - "dat dit bestand (met als extentie .sqlite) in de tweede stap geselecteerd " - "wordt. Let wel op: opnieuw uitgevoerde stappen en leggerkeuzes zullen " - "bestaande data overschrijven.") + self.msg_upper_row.setInformativeText( + "Voor de toewijzing van leggerprofielen wordt een aparte 'leggerdatabase' " + "gemaakt. Deze database is een aparte .sqlite bestand waar data uit " + "DAMO en de Hydrologendatabase (HDB) gecombineerd wordt als randvoorwaarden " + "voor de leggerprofielen, zoals breedte en talud per hydro-object.\n" + "Wanneer een nieuwe leggerdatabase gemaakt moet worden, selecteer dan bij " + "voorkeur de DAMO en HDB die ook voor de opbouw van het 3di model zijn " + "gebruikt.\n" + "Is er al een 'leggerdatabase' aangemaakt, sla deze stap dan over en zorg " + "dat dit bestand (met als extentie .sqlite) in de tweede stap geselecteerd " + "wordt. Let wel op: opnieuw uitgevoerde stappen en leggerkeuzes zullen " + "bestaande data overschrijven.") self.box_explanation.addWidget(self.msg_upper_row) def select_DAMO(self): @@ -156,7 +162,6 @@ def select_spatialite(self): def create_spatialite_database(self): # todo: # feedback if input is incorrect. - # feedback of proces en of het gelukt is settings = QSettings('leggertool', 'filepaths') try: @@ -183,8 +188,13 @@ def create_spatialite_database(self): filepath_DAMO = self.var_text_DAMO.text() filepath_HDB = self.var_text_HDB.text() - if not filepath_DAMO: - # todo: raise error + if not os.path.exists(filepath_DAMO): + messagebar_message( + 'Aanmaken leggerdatabase mislukt', + 'Opgegeven DAMO database bestaat niet', + level=2, + duration=10) + raise Exception('Geselecteerde DAMO database mist') legger_class = CreateLeggerSpatialite( @@ -198,6 +208,16 @@ def create_spatialite_database(self): # set root_tool as last, because this triggers other actions self.root_tool.polder_datasource = database + # create views + con_legger = dbapi.connect(self.root_tool.polder_datasource) + create_legger_views(con_legger) + + messagebar_message( + 'Aanmaken leggerdatabase', + 'Aanmaken leggerdatabase is gelukt', + level=3, + duration=10) + def setup_ui(self): self.setMinimumWidth(700) @@ -281,7 +301,6 @@ def setup_ui(self): self.hbox_LDB.addWidget(self.load_leggerdatabase_button) self.box_leggerdatabase_input.addLayout(self.hbox_LDB) - # Create groupbox and add H or VBoxes to it self.groupBox_explanation = QtGui.QGroupBox(self) self.groupBox_explanation.setTitle("Uitleg") @@ -316,4 +335,4 @@ def retranslateUi(self, Dialog): self.load_HDB_dump_button.setText(_translate("Dialog", "Selecteer HDB gdb", None)) self.create_leggerdatabase_button.setText(_translate("Dialog", "Aanmaken database", None)) self.load_leggerdatabase_button.setText(_translate("Dialog", "Selecteer leggerdb", None)) - self.cancel_button.setText(_translate("Dialog", "Sluiten", None)) \ No newline at end of file + self.cancel_button.setText(_translate("Dialog", "Sluiten", None)) From ad7bab522e75dea82d2834514dfd577c49c7e075 Mon Sep 17 00:00:00 2001 From: Bastiaan Roos Date: Thu, 2 May 2019 19:44:46 +0200 Subject: [PATCH 2/2] remove unused code --- utils/profile_match.py | 507 ----------------------------------------- 1 file changed, 507 deletions(-) delete mode 100644 utils/profile_match.py diff --git a/utils/profile_match.py b/utils/profile_match.py deleted file mode 100644 index f84bd0f..0000000 --- a/utils/profile_match.py +++ /dev/null @@ -1,507 +0,0 @@ -# profile_match.py -# onderdeel van leggertool -# plaatst theoretische profielen zo goed mogelijk binnen gemeten profielen - -from pyspatialite import dbapi2 as sql -import shapely -import shapely.geometry -import future_builtins -import matplotlib -matplotlib.use('AGG') -from matplotlib.pyplot import savefig -import matplotlib.pyplot as plt - -from descartes import PolygonPatch -import logging - -logger = logging.getLogger('legger.utils.profile_match') - - -def mk_pro_x_hy_kompas(cur, straal=0.5, aantalstappen=90, srid=28992): - """ Voor alle gemeten profielen wordt het snijpunt met het bijbehorende hydroobject gemaakt en het azimut - (kompasrichting) van het hydroobject ter plekke. - Invoer: cur = een cursor naar de database met gemetenprofielen, hydroobjecten en theoretische profielen - straal = de straal van een semi-cirkel met als middelpunt het snijpunt van gemetenprofiel en hydroobject - aantalstappen = het aantal lijnstukken van de semi-cirkel - srid = het nummer van het geografische referentiestelsel; 28992 = Rijksdriehoekmeting - Uitvoer: gevulde tabel pro_x_hy_kompas in de database - - Op de loodlijn op het lijnstuk van het hydroobject ter plekke van het snijpunt, worden later het gemeten - profiel en de theoretische profielen geprojecteerd. - Het snijpunt met de kompasrichting wordt opgeslagen in tabel pro_x_hy_kompas in de database. - - Het snijpunt wordt bepaald met de spatialite functie Intersection(profiellijn,hydroobject) - De kommpasrichting wordt bepaald met de spatialite functie Azimuth(punt1,punt2) - punt1 en punt2 worden bepaald door Intersection(profiellijn, cirkeltje_om_snijpunt) - cirkeltje_om_snijpunt wordt gemaakt door Buffer(snijpunt, straal, aantalstappen) - straal is een variabele evenals aantal stappen""" - - cur.execute('drop table if exists pro_x_hy_kompas') - cur.execute('create table pro_x_hy_kompas (pro_id bigint primary key, ovk_ovk_id bigint, kompas float,' - 'CONSTRAINT fk_pro FOREIGN KEY (pro_id) REFERENCES pro(pro_id), ' - 'CONSTRAINT fk_hy FOREIGN KEY (ovk_ovk_id) REFERENCES hydroobject(id)) ') - cur.execute('select DiscardGeometryColumn("pro_x_hy_kompas","geometry")') - cur.execute('select AddGeometryColumn("pro_x_hy_kompas", "geometry", %d, "POINT")' % srid) - cur.execute('insert into pro_x_hy_kompas (pro_id, ovk_ovk_id, geometry, kompas)' - ' select pro_id, ovk_ovk_id, Intersection(pro.GEOMETRY, hydroobject.GEOMETRY),' - 'Azimuth(' - 'PointN(Intersection(pro.GEOMETRY,Buffer(Intersection(pro.GEOMETRY, hydroobject.GEOMETRY),%f,%d)),1), ' - 'PointN(Intersection(pro.GEOMETRY,Buffer(Intersection(pro.GEOMETRY, hydroobject.GEOMETRY),%f,%d)),2)) ' - 'from pro inner join hydroobject on ' - '(pro.ovk_ovk_id=hydroobject.id)' % (straal, aantalstappen, straal, aantalstappen)) - return - - -def peilperprofiel(cur, peilcriterium="min", debug=0): - """ Haal per gemeten profiel het heersende peil op. - Invoer: cur = een cursor naar de database met gemetenprofielen, hydroobjecten en theoretische profielen - peilcriterium = vlag met waarden min of max om resp.de minimale of maximale waterhoogte te kiezen - Uitvoer: een dictionary met als sleutel het id van het profiel en als waarden het id van het hydroobject en peil - - In versie 0 wordt domweg op grond van de administratieve joins de waterhoogte uit de tabel streefpeilen gebruikt - Dit kan later zo nodig verfijnt worden met een spatial join - In het testgebied Geestmerambacht levert de gebruikte query 43 null waarden op voor waterhoogte - steekproeven geven aan dat dit vooral komt doordat profielen geselecteerd zijn met een buffer - rond het gebied """ - if peilcriterium != 'min' and peilcriterium != 'max': - peilcriterium = 'min' - if debug: - logger.debug("PEILCRITERIUM AANGEPAST NAAR %s", peilcriterium) - q = '''select pro.pro_id, pro.ovk_ovk_id, %s(streefpeil.waterhoogte) from - pro left outer join hydroobject on (pro.ovk_ovk_id = hydroobject.id) - left outer join peilgebiedpraktijk on (hydroobject.ws_in_peilgebied = peilgebiedpraktijk.code) - left outer join streefpeil on (peilgebiedpraktijk.id=streefpeil.peilgebiedpraktijkid) - group by pro.pro_id, pro.ovk_ovk_id - ''' % peilcriterium - prof = {} - for r in cur.execute(q): - prof[r[0]] = (r[1], r[2]) - q = '''insert into hyob_voorkeurpeil select hydroobject.id, %s(streefpeil.waterhoogte) from - hydroobject left outer join peilgebiedpraktijk on (hydroobject.ws_in_peilgebied = peilgebiedpraktijk.code) - left outer join streefpeil on (peilgebiedpraktijk.id=streefpeil.peilgebiedpraktijkid) - group by hydroobject.id''' % peilcriterium - cur.execute(q) - if debug: - logger.debug("aantal gemeten profielen in een hydro_object met een peil: %d ", len(prof)) - return prof - - -def haal_meetprofielen1(cur, profielsoort="Z1", peilcriterium="min", debug=0): - """ Haal de gemeten profieelpunten op uit de database voor de profielsoort vastebodem (Z1) - Invoer: cur = een cursor naar de database met gemetenprofielen, hydroobjecten en theoretische profielen - profielsoort = de code voor de harde bodem - peilcriterium = vlag met waarden min of max om resp.de minimale of maximale waterhoogte te kiezen - Uitvoer: een dictionary met als sleutel de profielids van gemeten profielen - met per profielid: - het hydroid (hydro object id) - het peil - de punten van het gemeten profiel + extra begin- en eindpunt 100m hoger """ - prof = {} - peilvanprofiel = peilperprofiel(cur, peilcriterium, debug) # per profielid het hydroobject_id en het peil - #q = '''select "c"."pro_pro_id", "b"."iws_volgnr", X("a"."GEOMETRY"), Y("a"."GEOMETRY"), - # "b"."iws_hoogte", a."OGC_FID" - # from "profielpunten" as "a" - # join "pbp" as "b" on ("a"."OGC_FID" = "b"."OGC_FID") - # join "prw" as "c" on ("b"."OGC_FID" = "c"."OGC_FID") - # where "c"."pro_pro_id" = %d and "c"."osmomsch" = "%s" - # order by "b"."iws_volgnr" - # ''' - #q = '''select pro.pro_id, pbp.iws_volgnr, - #X(profielpunten.GEOMETRY), Y(profielpunten.GEOMETRY), - #pbp.iws_hoogte, pro.ovk_ovk_id, pro.proident - #from pro inner join prw on (pro.pro_id=prw.pro_pro_id) - #inner join pbp on (prw.prw_id=pbp.prw_prw_id) - #inner join profielpunten on (pbp.pbp_id=profielpunten.pbp_pbp_id) - #where pro.pro_id = %d and prw.osmomsch="%s" - #order by pbp.iws_volgnr''' - q = '''select pro.pro_id, pbp.iws_volgnr, - X(profielpunten.GEOMETRY), Y(profielpunten.GEOMETRY), - pbp.iws_hoogte, pro.ovk_ovk_id, pro.proident - from pro inner join prw on (pro.pro_id=prw.pro_pro_id) - inner join pbp on (prw.OGC_FID=pbp.OGC_FID) - inner join profielpunten on (pbp.OGC_FID=profielpunten.OGC_FID) - where pro.pro_id = %d and prw.osmomsch="%s" - order by pbp.iws_volgnr''' - for proid in peilvanprofiel: - if peilvanprofiel[proid][1] is not None: # er is een peil ! - prof[proid] = {} # veronderstelling 1 profiel per hydroobject is fout!!!! nu gewoon proid invullen - prof[proid]['hydroid'] = peilvanprofiel[proid][0] - prof[proid]['peil'] = peilvanprofiel[proid][1] - s = 0 - for r in cur.execute(q % (proid, profielsoort)): - if s == 0: # beginconditie, eerste punt van profiel 1000m hoger, zodat altijd boven peil gestart wordt - prof[proid]['orig'] = [[r[2], r[3], max(r[4], peilvanprofiel[proid][1]) + 1000.0, 0]] - s += 1 - prof[proid]["orig"].append([r[2], r[3], r[4], r[5]]) - # laatste punt van profiel, 1000 m hoger - prof[proid]["orig"].append([r[2], r[3], max(r[4], peilvanprofiel[proid][1]) + 1000.0, 0]) - prof[proid]['proident'] = r[6] - if debug: - logger.debug('aantal profielen in hydro-objecten met een peil: %d', len(prof)) - return prof - - -def interpoleerafstand(l, r, p): - """interpoleer de afstand op grond van de hoogtes - Invoer: l = linker list met afstand, x, y en z - r = rechter list met [a, x, y, z] - p = z waarde tussen l[3] en r[3] in - - Uitvoer: tuple van afstand a en hoogte z - """ - factor = (p - l[3])/(r[3]-l[3]) - return l[0] + factor * (r[0] - l[0]), p - - -def projecteerprofielen(prof, projectie="eindpunt", debug=0): - """ Verrijk profielen met de projectie op een rechte lijn - Invoer: prof = dictionary met gemeten profielen; punten staan onder key "orig" - projectie = keuze soort projectie van de xy punten van het gemeten profiel; waarden: - eindpunt: op de rechte tussen de eindpunten van het gemeten profiel - loodlijn: op de loodlijn op het lijnstuk van het hydroObject tpv het kruispunt - van het lijnstuk van het hydroobject met de gemeten profiellijn (pro) - Uitvoer: prof verrijkt met de key "proj" met daarin een list van lists van afstand-geprojecteerd, - x-geprojecteerd, y-geprojecteerd en diepte - In eerste instantie is alleen eindpunt geimplementeerd! - """ - if projectie == 'eindpunt': - for proid in prof: - lijn = shapely.geometry.LineString([(prof[proid]['orig'][0][0], prof[proid]['orig'][0][1]), - (prof[proid]['orig'][-1][0], prof[proid]['orig'][-1][1])]) - prof[proid]['proj'] = [] - for p in prof[proid]['orig']: - afstand = lijn.project(shapely.geometry.Point((p[0], p[1]))) - pr = lijn.interpolate(afstand) - prof[proid]['proj'].append([afstand, pr.x, pr.y, p[2], p[3]]) - else: - if debug: - logger.debug("PROJECTIECRITERIUM NIET GELDIG %s, GEEN PROJECTIE!", projectie) - return prof - - -def verrijkgemprof(cur, prof): - """Verrijk de tabel profielpunten met de afstanden zoals die in projecteerprofielen berekend zijn""" - q = 'update profielpunten set afstand = %f where OGC_FID= %d' - for proid in prof: - for p in prof[proid]['proj']: - if p[4] !=0 : - cur.execute(q % (p[0], p[4])) - return - -def mkmogelijktheoprofiel(talud, waterdiepte, bodembreedte, peil): - """Maak een shapely polygoon van het theoretisch profiel """ - return shapely.geometry.Polygon([(0, peil), (talud * waterdiepte, peil - waterdiepte), - (talud * waterdiepte + bodembreedte, peil - waterdiepte), - (talud * waterdiepte + bodembreedte + talud * waterdiepte, peil)]) - - -def mkrechthoekondertheoprofiel(rhlb, rhrb): - """Maak een shapely polygoon van een rechthoek onder het theoretisch profiel""" - return shapely.geometry.Polygon([rhlb, (rhlb[0], rhlb[1]-1000.0) , (rhrb[0], rhrb[1]-1000.0), rhrb]) - - -def grootste(xy): - index = 0 - m = 0 - for i in xy: - if len(xy[i]) > 2: - if (xy[i][-1][0] - xy[i][0][0]) > m: - index = i - m = xy[i][-1][0] - xy[i][0][0] - return index - - -def mkgemprof(axyzlist, peil): - """ Maak een shapely polygoon van het gemeten profiel (afstanden en diepte), doorsnijden met het peil - helaas is split pas vanaf shapely versie 1.6 aanwezig (qgis 2.18 heeft shapely 1.2) daarom - - Invoer: list van lists met afstand, x, y en z - peil - Uitvoer: shapely polygon - """ - xy = {} - tel = 0 - xy[tel] = [] - links = axyzlist[0] - positie = 'b' - for c in axyzlist[1:]: - if positie == 'b': - if c[3] < peil: - xy[tel].append(interpoleerafstand(links, c, peil)) - xy[tel].append((c[0], c[3])) - positie = "o" - elif c[3] == peil: - xy[tel].append((c[0], c[3])) - positie = "o" - else: - if c[3] > peil: - xy[tel].append(interpoleerafstand(links, c, peil)) - positie = "b" - tel += 1 - xy[tel] = [] - elif c[3] == peil: - xy[tel].append((c[0], c[3])) - positie = "b" - tel += 1 - xy[tel] = [] - else: - xy[tel].append((c[0], c[3])) - links = c - if tel > 0: - tel = grootste(xy) - return shapely.geometry.Polygon(xy[tel]) - - -def prof_in_prof(profgem, proftheo, aantstap=100, delta=0.001, obdiepte=0.3, debug=0): - """Het theoretisch profiel wordt in stapjes verschoven over het gemeten profiel (te beginnen links van het - gemeten profiel zonder overlap, tot en met rechts van het gemeten profiel zonder overlap). Indien het - theoretisch profiel nergens past binnen het gemeten profiel is de plek met het maximale oppervlak van de - intersectie van het theoretisch met het gemeten profiel de optimale plek voor het theoretisch profiel, - tenzij er een traject is met een gelijk maximaal oppervlak; in dat geval wordt het midden van dit traject de - optimale plek voor het theoretisch profiel. - Wanneer het oppervlak van deze intersectie gelijk is aan het oppervlak van het theoretisch profiel, dan - past het theoretisch profiel volkomen binnen het gemeten profiel. In dat geval wordt het midden van het - traject waarvoor het intersectie oppervlak gelijk is aan het oppervlak van het theoretisch profiel, de - optimale plek voor het theoretisch profiel - Invoer: profgem: het shapely polygoon van het gemeten profiel - proftheo: het shapely polygoon van het theoretisch profiel - aantstap: het aantal stappen dat gebruikt wordt voor de verschuiving van het theoretisch profiel - delta: het acceptabele verschil om vast te stellen of twee oppervlakken gelijk zijn - obdiepte: de diepte waarop de beschikbare overbreedte bepaald wordt - Uitvoer: fit: de fractie van het oppervlak van het theoretisch profiel dat past binnen het gemeten - optimaal: de optimale verschuiving (afstand) van het theoretisch profiel gemeten vanaf - de start van het gemeten profiel - fractie: 1 - de fractie van het gemeten profiel dat bedekt wordt door het theoretisch profiel - overdiepte: de gemiddelde afstand onder rechte stuk van het theoretisch profiel tot het gemeten profiel - linksover: de afstand tussen het gemeten profiel en het theoretisch profiel op obdiepte links - rechtsover: de afstand tussen het gemeten profiel en het theoretisch profiel op obdiepte rechts""" - waterbreedte_gemeten = profgem.bounds[2] - profgem.bounds[0] - waterbreedte_theo = proftheo.bounds[2] - proftheo.bounds[0] - waterdiepte_theo = proftheo.bounds[3] - proftheo.bounds[1] - rhlb = proftheo.exterior.coords[1] # de linkerbovenhoek resp rechterbovenhoek van een rechthoek onder - rhrb = proftheo.exterior.coords[2] # het rechte stuk van het theoretisch profiel (tbv overdiepte) - oblijn = shapely.geometry.LineString([(0,profgem.bounds[3]-obdiepte), - (waterbreedte_gemeten + 2 * waterbreedte_theo, profgem.bounds[3]-obdiepte)]) - clijn = profgem.intersection(oblijn) # oblijn tbv overbreedte, clijn controle lijnstuk tbv overbreedte - gemprof = shapely.affinity.translate(profgem, -profgem.bounds[0], 0.0, 0.0) # op nul meter laten beginnen - profzoek = shapely.affinity.translate(proftheo, -waterbreedte_theo, 0.0, 0.0) # verschuif naar nul overlap - stap = (waterbreedte_gemeten + waterbreedte_theo + waterbreedte_theo) / aantstap - zoekopp = profzoek.area # het oppervlak van het theoretisch profiel - maxopp = -9.9 # het maximale oppervlak van de intersectie - traject = '' # in traject komt per stap een 0 of 1 (1 indien maxopp == zoekopp) - optimaal = -waterbreedte_theo # dit is de start van het theoretisch profiel, overlap is nul! - fit = 0 # de verhouding maxopp / zoekopp (een goodness of fit) - if debug: - logger.debug("wb_gem: %.2f, wb_theo: %.2f; stap: %.2f, zoekopp: %.3f", - waterbreedte_gemeten, waterbreedte_theo, stap, zoekopp) - for i in xrange(aantstap): - profzoek = shapely.affinity.translate(profzoek, stap, 0.0, 0.0) - inter = gemprof.intersection(profzoek) - if abs(zoekopp-inter.area) < delta: # opp intersectie == zoekopp dus past profzoek volledig in gemprof - optimaal = profzoek.bounds[0] - maxopp = inter.area - traject += '1' - if debug: - logger.debug('Volledig: optimaal: %f; maxopp: %f', optimaal, maxopp) - else: - if inter.area == maxopp: - traject += '2' - elif inter.area > maxopp: # opp intersectie groter dan voorgaand oppervlak - optimaal = profzoek.bounds[0] - maxopp = inter.area - traject = traject.replace('2', '0') # evt oud traject met kleiner oppervlak weghalen - traject += '2' - if debug: - logger.debug('Niet vol: optimaal: %f; maxopp: %f', optimaal, maxopp) - else: - traject += '0' - fit = maxopp / zoekopp # de best fit - fractie = (gemprof.area-maxopp) / gemprof.area # de overblijvende fractie oppervlak van het gemetenprofiel - if debug: - logger.debug('Na loop: Fit: %.3f; optimaal: %f', fit, optimaal) - logger.debug(traject) - - if traject.find('1') > 0: # er zijn 1 tekens in traject - traject = traject.replace('2', '0') # evt oud traject met kleiner oppervlak weghalen - zoek = max(traject.split('0')) # in zoek komt de eerste langste reeks 1 tekens in traject - astap = traject.find(zoek) + len(zoek)/2.0 + 1 - optimaal = astap * stap - optimaal -= waterbreedte_theo # corrigeer voor de verschuiving van het theoretisch profiel - elif traject.find('2') > 0: # er zijn trajecten met een kleiner oppervlak - zoek = max(traject.split('0')) # in zoek komt de eerste langste reeks 1 tekens in traject - astap = traject.find(zoek) + len(zoek)/2.0 + 1 - optimaal = astap * stap - optimaal -= waterbreedte_theo # corrigeer voor de verschuiving van het theoretisch profiel - optimaal += profgem.bounds[0] # corrigeer voor de verschuiving van het gemeten profiel - - roth = shapely.affinity.translate(mkrechthoekondertheoprofiel(rhlb, rhrb), optimaal, 0.0, 0.0) - poloth = profgem.intersection(roth) - overdiepte = poloth.area/(rhrb[0]-rhlb[0]) # oppervlak gedeeld door breedte geeft diepte - - linksover = 0.0 - rechtsover = 0.0 - restbak = profgem.difference(shapely.affinity.translate(proftheo,optimaal, 0.0, 0.0)) - hlijn = restbak.intersection(oblijn) # restbak is restant gemeten profiel min theor. profiel - if obdiepte < (profgem.bounds[3] - profgem.bounds[1]): # de obdiepte is hoger dan de bodem van het gemetenprofiel - try: - xlinks = clijn.coords[0][0] - xrechts = clijn.coords[1][0] - - try: - if len(hlijn) == 2: # is hlijn een MultiLineString van twee LineStrings - if (hlijn[0].coords[0][0] == xlinks) and (hlijn[1].coords[1][0] == xrechts): - linksover = hlijn[0].length - rechtsover = hlijn[1].length - elif len(hlijn) == 3: # is hlijn een MultiLineString van drie LineStrings - if debug: - logger.debug("hlijn 3 stuks, xl, xr, hl", xlinks, xrechts, hlijn) - if (hlijn[0].coords[1][0] == xlinks) and (hlijn[2].coords[1][0] == xrechts): - linksover = hlijn[1].length - rechtsover = hlijn[2].length - else: - logger.debug(hlijn) - pass - except: - if debug: - logger.debug("hlijn is geen MultiLindeString") - pass # hlijn is geen MultiLineString - except: - if debug: - logger.debug("clijn is geen LineString") - pass # clijn is geen LineString (mogelijk een MultiLineString) - return fit, optimaal, fractie, overdiepte, linksover, rechtsover - - -def altertable(cur, tabelnaam, veldnaam, veldtype): - """"Primitieve alter table voor float, integer, double met controle of het veld al bestaat""" - gevonden = 0 - for r in cur.execute('pragma table_info("%s")' % tabelnaam): - if r[1] == veldnaam: - gevonden = 1 - if not gevonden: - cur.execute('alter table "%s" add column "%s" "%s"' % (tabelnaam, veldnaam, veldtype) ) - return - - -def maaktabellen(cur): - """"Maak de tabellen met verrijkte platgeslagen profielen tbv presentatie - presentatie: klik op de kaart nabij een hydroobject en een profiel (selecteer dichtstbijzijnde, - een hydroobject kan meer gemeten profielen hebben!! - tabel hyob_voorkeurpeil geeft op grond van het id van het hydroobject het gekozen peil (kan natuurlijk - ook als view, maar dit zal sneller zijn, geen idee of het van belang is) - tabel profielfiguren is een platte tabel met alle info voor figuren met gemeten en theoretische profielen - met infor over fit, overdiepte enz enz. - - Aanpassing van bestaande tabellen: - hydroobject met voorkeurpeil - profielpunten met afstand - """ - altertable(cur, "hydroobject", "voorkeurpeil", "float") - altertable(cur, "profielpunten", "afstand", "float") - - cur.execute('drop table if exists hyob_voorkeurpeil') - cur.execute('create table hyob_voorkeurpeil (id integer primary key, voorkeurpeil float)') - cur.execute('drop table if exists profielfiguren') - cur.execute('drop index if exists profielfiguren0') - cur.execute('drop index if exists profielfiguren1') - cur.execute('create table profielfiguren(id_hydro integer, profid varchar(16), type_prof char(1), coord text, ' - 'peil float, t_talud float, t_waterdiepte float, t_bodembreedte float, t_fit float, t_afst float, ' - 'g_rest float, t_overdiepte float, t_overbreedte_l float, t_overbreedte_r float)') - cur.execute('vacuum') - return - - -def controlefig(gemprof, theoprof, afstand, fit, fractie, overdiepte, overlinks, overrechts, hydroid, profid, - talud, diepte, breedte, peil): - fnm = 'cf/%d_%d_%d_%.2f_%.2f' % (profid, hydroid, talud, diepte, breedte) - fnm = fnm.replace('.', ',') - txt1 = 'Hydroid: %d, profid: %d, peil: %.2f; Talud: %d; Waterdiepte: %.2f; Bodembreedte: %2f.' % \ - (hydroid, profid, peil, talud, diepte, breedte) - txt2 = 'Fit: %.2f; Fractie %.2f, Overdiepte %.3f; Overbreedte links: %.1f; Overbreedte rechts: %.1f.' % \ - (fit, fractie, overdiepte, overlinks, overrechts) - blue = '#6699cc' - orange = '#cc9933' - fig, ax = plt.subplots() - fig.text(0.95, 0.15, txt1, - fontsize=7, color='black', - ha='right', va='bottom', alpha=0.9) - fig.text(0.95, 0.10, txt2, - fontsize=7, color='black', - ha='right', va='bottom', alpha=0.9) - p2 = PolygonPatch(gemprof, fc=blue, ec=blue, alpha=0.5) - ax.add_patch(p2) - p1 = PolygonPatch(shapely.affinity.translate(theoprof, afstand, 0.0, 0.0), fc=orange, ec=orange, alpha=0.5) - ax.add_patch(p1) - ax.axis('scaled') - savefig(fnm, dpi=200) - plt.close() - return - -def doe_profinprof(cur0, cur1, aantalstappen=200, precisie=0.0001, codevastebodem="Z1", peilcriterium="min", - projectiecriterium="eindpunt", obdiepte=0.3, debug=0): - try: - maaktabellen(cur0) - - gemetenprofielen = haal_meetprofielen1(cur0, codevastebodem, peilcriterium, debug) - gemetenprofielen = projecteerprofielen(gemetenprofielen, projectiecriterium, debug) - verrijkgemprof(cur0, gemetenprofielen) - - # alleen theoretische profielen die liggen in hydro-objecten waar ook gemeten profielen zijn ophalen - q = """select ROWID, id, talud, waterdiepte, bodembreedte from profiel_varianten where hydro_object_id='%s'""" - qm = """insert into profielfiguren (id_hydro, profid, type_prof, coord, peil) values(%d, "%s", "m", "%s", %f)""" - qt = """insert into profielfiguren (id_hydro, profid, type_prof, coord, peil, t_talud, t_waterdiepte, t_bodembreedte, - t_fit, t_afst, g_rest, t_overdiepte, t_overbreedte_l, t_overbreedte_r) values - (%d, "%s", "t", "%s", %f, %f, %f, %f, %f, %f, %f, %f, %f, %f)""" - for profielid in gemetenprofielen.keys(): - # mkgemprof aanroepen met list van lists met afstand, x, y, z (geprojecteerd); en het peil; - # levert een shapely polygoon - gemprofshapely = mkgemprof(gemetenprofielen[profielid]['proj'], gemetenprofielen[profielid]['peil']) - - h = qm % (gemetenprofielen[profielid]['hydroid'], gemetenprofielen[profielid]['proident'], - gemprofshapely.wkt, gemetenprofielen[profielid]['peil']) - cur1.execute(h) - for theo_data in cur0.execute(q % gemetenprofielen[profielid]['hydroid']): - # mkmogelijkprofiel aanroepen met talud, waterdiepte, bodembreedte en peil, levert een shapely polygon - theoprofshapely = mkmogelijktheoprofiel(theo_data[2], theo_data[3], theo_data[4], - gemetenprofielen[profielid]['peil']) - # prof_in_prof aanroepen met gemetenprofiel, theoretisch profiel aantal stappen en aanvaardbaar verschil - # (profielen bestaan uit shapely polygons) - fit, afstand, fractie, overdiepte, overlinks, overrechts = \ - prof_in_prof(gemprofshapely, theoprofshapely, aantalstappen, precisie, obdiepte, debug) - cur1.execute(qt % (gemetenprofielen[profielid]['hydroid'], theo_data[1], - shapely.affinity.translate(theoprofshapely, afstand, 0.0, 0.0).wkt, - gemetenprofielen[profielid]['peil'], theo_data[2], theo_data[3], theo_data[4], - fit, afstand, fractie, overdiepte, overlinks, overrechts)) - if debug: - controlefig(gemprofshapely, theoprofshapely, afstand, fit, fractie, overdiepte, overlinks, - overrechts, gemetenprofielen[profielid]['hydroid'], profielid, - theo_data[2], theo_data[3], theo_data[4], gemetenprofielen[profielid]['peil']) - cur0.execute('create index profielfiguren0 on profielfiguren(id_hydro)') - cur0.execute('create index profielfiguren1 on profielfiguren(profid)') - cur0.execute('vacuum') - resultaat = "klaar" - except: - resultaat = "Niet correct, run evt opnieuw met debug = 1 om te onderzoeken" - return resultaat - - -# start van het programma -# Variabelen -aantalstappen = 200 # Het aantal stappen dat gebruikt wordt door prof_in_prof (stapgrootte is: - # (breedte gemeten profiel + 2 * breedte theoretisch profiel) / aantalstappen -precisie = 0.0001 # De afwijking die aanvaardbaar is voor de vaststelling van gelijkheid van oppervlakken -codevastebodem = "Z1" # de code in de profieldata voor de vaste bodem -peilcriterium = "min" # het criterium om het peil voor een hydroobject te kiezen, geldige waarden min en max -projectiecriterium = 'eindpunt' # het criterium voor de projectie van de profielpunten, geldige waarden eindpunt -obdiepte = 0.3 # de diepte waarop de overbreedte bepaald moet worden -debug = 0 # vlag om extra output te genereren, geldige waarden 0 of 1 - -con = sql.connect('../tests/data/test_spatialite_with_theoretical_profiles_mp_20171206.sqlite') -cur0 = con.cursor() -cur1 = con.cursor() - -resultaat = doe_profinprof(cur0, cur1, aantalstappen, precisie, codevastebodem, peilcriterium, - projectiecriterium, obdiepte, debug) -logger.debug(resultaat) -con.commit() -con.close() \ No newline at end of file