Add discharge wells, linesinks

Make .renderer a classmethod instead of property

gistim/compute.py

@@ -1,7 +1,7 @@
 import json
 import pathlib
 from collections import defaultdict
-from typing import Dict, Union
+from typing import Any, Dict, Union
 
 import numpy as np
 import pandas as pd
@@ -112,6 +112,12 @@ def timml_head_observations(
     return pd.DataFrame(d)
 
 
+def timml_discharges(elements: Dict[str, Any]) -> Dict[str, pd.DataFrame]:
+    discharges = {}
+    for layername, content in elements.items():
+        first_element = content[0]
+
+
 def timml_headgrid(model: timml.Model, xmin, xmax, ymin, ymax, spacing) -> xr.DataArray:
     """
     Compute the headgrid of the TimML model, and store the results
@@ -146,6 +152,47 @@ def timml_headgrid(model: timml.Model, xmin, xmax, ymin, ymax, spacing) -> xr.Da
     )
 
 
+def extract_discharges(elements, nlayers):
+    tables = {}
+    for layername, content in elements.items():
+        sample = content[0]
+
+        if isinstance(sample, timml.WellBase):
+            table_rows = []
+            for well in content:
+                row = {f"discharge_layer{i}": q for i, q in enumerate(well.discharge())}
+                row["geometry"] = {
+                    "type": "Point",
+                    "coordinates": [well.xc[0], well.yc[0]],
+                }
+                table_rows.append(row)
+            tables[f"discharge-{layername}"] = pd.DataFrame.from_records(table_rows)
+
+        elif isinstance(sample, (timml.LineSinkDitchString, timml.HeadLineSinkString)):
+            table_rows = []
+            for linestring in content:
+                # Workaround for: https://github.com/mbakker7/timml/issues/85
+                discharges = np.zeros((linestring.nls, nlayers))
+                Q = (
+                    (linestring.parameters[:, 0] * linestring.dischargeinf())
+                    .reshape((linestring.nls, linestring.nlayers, linestring.order + 1))
+                    .sum(axis=2)
+                )
+                discharges[:, linestring.layers[0]] = Q
+
+                xy = linestring.xy
+                for q_layered, vertex0, vertex1 in zip(discharges, xy[:-1], xy[1:]):
+                    row = {f"discharge_layer{i}": q for i, q in enumerate(q_layered)}
+                    row["geometry"] = {
+                        "type": "LineString",
+                        "coordinates": [vertex0, vertex1],
+                    }
+                    table_rows.append(row)
+            tables[f"discharge-{layername}"] = pd.DataFrame.from_records(table_rows)
+
+    return tables
+
+
 def compute_steady(
     path: Union[pathlib.Path, str],
 ) -> None:
@@ -167,11 +214,10 @@ def compute_steady(
     write_ugrid(head, crs, path)
 
     # Compute observations and discharge, and write to geopackage.
-    tables = {}
+    tables = extract_discharges(elements, timml_model.aq.nlayers)
     if observations:
         for layername, content in observations.items():
             tables[layername] = timml_head_observations(timml_model, content["data"])
-    # TODO: Compute discharges...
 
     if tables:
         write_geopackage(tables, crs, path)

gistim/geopackage.py

@@ -37,7 +37,7 @@ def create_gpkg_contents(
     return pd.DataFrame(
         data={
             "table_name": table_names,
-            "data_type": ["features"],
+            "data_type": "features",
             "identifier": table_names,
             "description": "",
             "last_change": pd.Timestamp.now(),

plugin/qgistim/core/elements/building_pit.py

@@ -78,9 +78,9 @@ class BuildingPit(AssociatedElement):
     schema = BuildingPitSchema()
     assoc_schema = AssociatedBuildingPitSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.polygon_renderer(
             color=TRANSPARENT_RED, color_border=RED, width_border="0.75"
         )
 

plugin/qgistim/core/elements/circular_area_sink.py

@@ -63,9 +63,9 @@ class CircularAreaSink(TransientElement):
     )
     schema = CircularAreaSinkSchema()
 
-    @property
-    def renderer(self):
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls):
+        return cls.polygon_renderer(
             color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75"
         )
 

plugin/qgistim/core/elements/constant.py

@@ -24,9 +24,9 @@ class Constant(Element):
     )
     schema = ConstantSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.marker_renderer(color=RED, name="star", size="5")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.marker_renderer(color=RED, name="star", size="5")
 
     def process_timml_row(self, row, other=None):
         x, y = self.point_xy(row)

plugin/qgistim/core/elements/domain.py

@@ -32,12 +32,12 @@ def __init__(self, path: str, name: str):
         self.timml_name = f"timml {self.element_type}:Domain"
         self.ttim_name = "ttim Computation Times:Domain"
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
         """
         Results in transparent fill, with a medium thick black border line.
         """
-        return self.polygon_renderer(
+        return cls.polygon_renderer(
             color="255,0,0,0", color_border=BLACK, width_border="0.75"
         )
 

plugin/qgistim/core/elements/element.py

@@ -28,7 +28,7 @@
   and edit its data.
 
 Some elements require specific rendering in QGIS (e.g. no fill polygons), which
-are supplied by the `.renderer` property.
+are supplied by the `.renderer()` staticmethod.
 
 The coupling of separate tables (geometry table and time series table) is only
 explicit in the Dataset Tree. The only way of knowing that tables are
@@ -284,8 +284,8 @@ def polygon_renderer(**kwargs):
         symbol = QgsFillSymbol.createSimple(kwargs)
         return QgsSingleSymbolRenderer(symbol)
 
-    @property
-    def renderer(self):
+    @classmethod
+    def renderer(cls):
         return None
 
     def timml_layer_from_geopackage(self) -> QgsVectorLayer:

plugin/qgistim/core/elements/head_line_sink.py

@@ -67,9 +67,9 @@ class HeadLineSink(TransientElement):
     )
     schema = HeadLineSinkSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.line_renderer(color=BLUE, width="0.75")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.line_renderer(color=BLUE, width="0.75")
 
     def process_timml_row(self, row, other=None):
         return {

plugin/qgistim/core/elements/headwell.py

@@ -63,9 +63,9 @@ class HeadWell(TransientElement):
     )
     schema = HeadWellSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.marker_renderer(color=BLUE, size="3")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.marker_renderer(color=BLUE, size="3")
 
     def process_timml_row(self, row, other=None) -> Dict[str, Any]:
         x, y = self.point_xy(row)

plugin/qgistim/core/elements/impermeable_line_doublet.py

@@ -25,9 +25,9 @@ class ImpermeableLineDoublet(Element):
     }
     schema = ImpermeableLineDoubletSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.line_renderer(color=RED, width="0.75")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.line_renderer(color=RED, width="0.75")
 
     def process_timml_row(self, row, other=None):
         return {

plugin/qgistim/core/elements/leaky_building_pit.py

@@ -81,9 +81,9 @@ class LeakyBuildingPit(AssociatedElement):
     schema = LeakyBuildingPitSchema()
     assoc_schema = AssociatedLeakyBuildingPitchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.polygon_renderer(
             color=TRANSPARENT_RED,
             color_border=RED,
             width_border="0.75",

plugin/qgistim/core/elements/leaky_line_doublet.py

@@ -27,9 +27,9 @@ class LeakyLineDoublet(Element):
     }
     schema = LeakyLineDoubletSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.line_renderer(color=RED, width="0.75", outline_style="dash")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.line_renderer(color=RED, width="0.75", outline_style="dash")
 
     def process_timml_row(self, row, other=None):
         return {

plugin/qgistim/core/elements/line_sink_ditch.py

@@ -67,9 +67,9 @@ class LineSinkDitch(TransientElement):
     )
     schema = LineSinkDitchSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.line_renderer(color=GREEN, width="0.75")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.line_renderer(color=GREEN, width="0.75")
 
     def process_timml_row(self, row, other=None):
         return {

plugin/qgistim/core/elements/observation.py

@@ -35,9 +35,9 @@ class Observation(TransientElement):
     transient_columns = ("timeseries_id",)
     schema = ObservationSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.marker_renderer(color=LIGHT_BLUE, name="triangle", size="3")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.marker_renderer(color=LIGHT_BLUE, name="triangle", size="3")
 
     def process_timml_row(self, row, other=None):
         x, y = self.point_xy(row)

plugin/qgistim/core/elements/polygon_area_sink.py

@@ -29,9 +29,9 @@ class PolygonAreaSink(Element):
     }
     schema = PolygonAreaSinkSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.polygon_renderer(
             color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75"
         )
 

plugin/qgistim/core/elements/polygon_inhomogeneity.py

@@ -77,9 +77,9 @@ class PolygonInhomogeneity(AssociatedElement):
     schema = PolygonInhomogeneitySchema()
     assoc_schema = AssociatedPolygonInhomogeneitySchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.polygon_renderer(
             color=TRANSPARENT_GREY, color_border=GREY, width_border="0.75"
         )
 

plugin/qgistim/core/elements/polygon_semi_confined_top.py

@@ -33,9 +33,9 @@ class PolygonSemiConfinedTop(Element):
     }
     schema = PolygonSemiConfinedTopSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.polygon_renderer(
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.polygon_renderer(
             color=TRANSPARENT_BLUE, color_border=BLUE, width_border="0.75"
         )
 

plugin/qgistim/core/elements/well.py

@@ -78,9 +78,9 @@ class Well(TransientElement):
     )
     schema = WellSchema()
 
-    @property
-    def renderer(self) -> QgsSingleSymbolRenderer:
-        return self.marker_renderer(color=GREEN, size="3")
+    @classmethod
+    def renderer(cls) -> QgsSingleSymbolRenderer:
+        return cls.marker_renderer(color=GREEN, size="3")
 
     def process_timml_row(self, row, other=None) -> Dict[str, Any]:
         x, y = self.point_xy(row)