diff --git a/gistim/__init__.py b/gistim/__init__.py index 83d8b71..413d7b5 100644 --- a/gistim/__init__.py +++ b/gistim/__init__.py @@ -1,13 +1,7 @@ -import pathlib -from typing import Dict, Union - import pkg_resources -from gistim import timml_elements, ttim_elements -from gistim.common import gridspec, model_specification -from gistim.compute import compute_steady, compute_transient -from gistim.data_extraction import as_aquifer_aquitard, layer_statistics -from gistim.ugrid import to_ugrid2d +import gistim.compute +import gistim.data_extraction # version try: @@ -15,54 +9,3 @@ except pkg_resources.DistributionNotFound: # package is not installed pass - - -def convert_to_script(inpath: str, outpath: str) -> None: - timml_spec, ttim_spec = model_specification(inpath, {}) - timml_script = timml_elements.convert_to_script(timml_spec) - try: - ttim_script = ttim_elements.convert_to_script(ttim_spec) - except Exception: - ttim_script = "" - - with open(outpath, "w") as f: - f.write(timml_script) - f.write("\n") - f.write(ttim_script) - - -def compute( - inpath: Union[pathlib.Path, str], - outpath: Union[pathlib.Path, str], - mode: str, - cellsize: float, - active_elements: Dict[str, bool], -) -> None: - """ - Compute the results of TimML model. - - The model is fully specified by the GeoPacakge dataset in the path. - - The extent of the head grids is read from a vector layer in the - GeoPackage file. - - Parameters - ---------- - path: Union[pathlib.Path, str] - Path to the GeoPackage file containing the full model input. - cellsize: float - Grid cell size of the computed output - - Returns - ------- - None - The result is written to a netCDF file. Its name is generated from - the geopackage name, and the requested grid cell size. - """ - if mode == "steady-state": - compute_steady(inpath, outpath, cellsize, active_elements) - elif mode == "transient": - compute_transient(inpath, outpath, cellsize, active_elements) - else: - raise ValueError(f'Invalid mode: {mode}: should be "steady" or "transient".') - return diff --git a/gistim/__main__.py b/gistim/__main__.py index 9fc1646..2e12e0a 100644 --- a/gistim/__main__.py +++ b/gistim/__main__.py @@ -81,21 +81,12 @@ def handle(line) -> None: # print(json.dumps(data, indent=4)) operation = data.pop("operation") - if operation == "compute": - gistim.compute( - inpath=data["inpath"], - outpath=data["outpath"], - cellsize=data["cellsize"], - mode=data["mode"], - active_elements=data["active_elements"], + gistim.compute.compute( + path=data["path"], + transient=data["transient"], ) - response = "Computation of {inpath} to {outpath}".format(**data) - elif operation == "convert": - inpath = data["inpath"] - outpath = data["outpath"] - gistim.convert_to_script(inpath, outpath) - response = "Conversion of {inpath} to {outpath}".format(**data) + response = "Computation of {path}".format(**data) elif operation == "extract": inpath = data["inpath"] outpath = data["outpath"] @@ -148,28 +139,8 @@ def extract(args) -> None: gistim.data_extraction.netcdf_to_table(inpath, outpath, wkt_geometry) -def convert(args) -> None: - """ - Convert a Geopackage into a Python script. - """ - inpath = args.inpath[0] - outpath = args.outpath[0] - gistim.convert_to_script(inpath, outpath) - - def compute(args) -> None: - jsonpath = args.jsonpath[0] - - with open(jsonpath, "r") as f: - data = json.loads(f.read()) - - gistim.compute( - inpath=data["inpath"], - outpath=data["outpath"], - cellsize=data["cellsize"], - mode=data["mode"], - active_elements=data["active_elements"], - ) + gistim.compute.compute(path=args.path[0], transient=args.transient[0]) return @@ -192,12 +163,10 @@ def compute(args) -> None: parser_extract.add_argument("outpath", type=str, nargs=1, help="outpath") parser_extract.add_argument("wkt", type=str, nargs=1, help="wkt") - parser_convert.set_defaults(func=convert) - parser_convert.add_argument("inpath", type=str, nargs=1, help="inpath") - parser_convert.add_argument("outpath", type=str, nargs=1, help="outpath") - parser_compute.set_defaults(func=compute) - parser_compute.add_argument("jsonpath", type=str, nargs=1, help="jsonpath") + parser_compute.add_argument("path", type=str, nargs=1, help="path to JSON file") + parser_compute.add_argument("--transient", action=argparse.BooleanOptionalAction) + parser.set_defaults(transient=False) # Parse and call the appropriate function args = parser.parse_args() diff --git a/gistim/common.py b/gistim/common.py deleted file mode 100644 index 4f918ee..0000000 --- a/gistim/common.py +++ /dev/null @@ -1,379 +0,0 @@ -""" -Common utilities used by conversions from GeoPackage (GPKG) layers to TimML and -TTim elements. - -The relation between a GPKG layer and the elements is not perfect. A GPGK layer -consists of a table, with optionally an associated geometry for every row. This -matches one to one for elements such as HeadLineSinkStrings, Wells, etc: for -these elements, one row equals one element. - -For elements such as PolygonImhomogenities, this is not the case. Every geometry -(a polygon) requires a table of its own. These tables are stored in associated -tables; their association is by name. - -For transient (TTim) elements, the same is true: elements require an additional -table for their timeseries data, which should require repeating the geometry -for every time step. In this package, we assume that any ttim element is -accompanied by a timml element; the QGIS plugin always sets up the layer ih -that manner. - -When processing a GPKG, we first parse the names, and group the different -tables together in the ElementSpecifications below. The ``timml_elements`` and -``ttim_elements`` then convert these grouped tables into ``timml`` and ``ttim`` -models. -""" -import numbers -import pathlib -import re -from collections import defaultdict -from functools import partial -from typing import Any, Dict, NamedTuple, Tuple, Union - -import fiona -import geopandas as gpd -import numpy as np -import pandas as pd - -FloatArray = np.ndarray - - -def filter_scalar_nan(value: Any) -> Any: - if isinstance(value, numbers.Real) and np.isnan(value): - return None - else: - return value - - -def filter_nan(row: Union[Dict, pd.Series]): - """ - Newer versions of geopandas return NaN rather than None for columns. - TimML & TTim expect None for optional values. - """ - return {k: filter_scalar_nan(v) for k, v in row.items()} - - -class ElementSpecification(NamedTuple): - elementtype: str - active: bool - dataframe: gpd.GeoDataFrame - associated_dataframe: gpd.GeoDataFrame - - -class TransientElementSpecification(NamedTuple): - elementtype: str - active: bool - dataframe: gpd.GeoDataFrame - steady_spec: ElementSpecification - - -class TimmlModelSpecification(NamedTuple): - aquifer: gpd.GeoDataFrame - elements: Dict[str, ElementSpecification] - domain: gpd.GeoDataFrame - - -class TtimModelSpecification(NamedTuple): - aquifer: gpd.GeoDataFrame - temporal_settings: gpd.GeoDataFrame - elements: Dict[str, ElementSpecification] - domain: gpd.GeoDataFrame - output_times: FloatArray - - -# Extract coordinates from geodataframe -# ------------------------------------- -def point_coordinates(dataframe) -> Tuple[FloatArray, FloatArray]: - return dataframe["geometry"].x, dataframe["geometry"].y - - -def remove_zero_length(coords: FloatArray): - dx_dy = np.diff(coords, axis=0) - notzero = (dx_dy != 0).any(axis=1) - keep = np.full(len(coords), True) - keep[1:] = notzero - return coords[keep] - - -def linestring_coordinates(row) -> FloatArray: - return remove_zero_length(np.array(row["geometry"].coords)) - - -def polygon_coordinates(row) -> FloatArray: - return remove_zero_length(np.array(row["geometry"].exterior.coords)) - - -# Parse GPKG content to Tim input -# ------------------------------- -def aquifer_data( - dataframe: gpd.GeoDataFrame, transient: bool = False -) -> Dict[str, Any]: - """ - Convert a table created by the QGIS plugin to the layer configuration and - keywords arguments as expected by TimML or TTim. - """ - # Make sure the layers are in the right order. - dataframe = dataframe.sort_values(by="layer").set_index("layer") - nlayer = len(dataframe) - # Deal with optional semi-confined top layer. - hstar = dataframe.loc[0, "semiconf_head"] - semi = pd.notnull(hstar) - kaq = dataframe["aquifer_k"].values - - if semi: - c = dataframe["aquitard_c"].values - porosity = np.empty(nlayer * 2) - z = np.empty(nlayer * 2 + 1) - z[0] = dataframe.loc[0, "semiconf_top"] - z[1::2] = dataframe["aquifer_top"].values - z[2::2] = dataframe["aquifer_bottom"].values - porosity[::2] = dataframe["aquitard_npor"].values - porosity[1::2] = dataframe["aquifer_npor"].values - topboundary = "semi" - storage_aquifer = dataframe["aquifer_s"].values - storage_aquitard = dataframe["aquitard_s"].values - else: - c = dataframe["aquitard_c"].values[1:] - z = np.empty(nlayer * 2) - z[::2] = dataframe["aquifer_top"].values - z[1::2] = dataframe["aquifer_bottom"].values - porosity = np.empty(nlayer * 2 - 1) - porosity[::2] = dataframe["aquifer_npor"].values - porosity[1::2] = dataframe["aquitard_npor"].values[1:] - topboundary = "conf" - storage_aquifer = dataframe["aquifer_s"].values - storage_aquitard = dataframe["aquifer_s"].values[1:] - - d = { - "kaq": kaq, - "z": z, - "c": c, - "topboundary": topboundary, - } - if transient: - d["Sll"] = storage_aquitard - d["Saq"] = storage_aquifer - # TODO: for now, assume there is always at least one specific yield - # this is the aquifer if conf, the aquitard if semi - d["phreatictop"] = True - else: - d["npor"] = porosity - d["hstar"] = hstar - - # For inhomogeneities - if "rate" in dataframe: - d["N"] = dataframe.loc[0, "rate"] - - filtered = {} - for k, value in d.items(): - if isinstance(value, np.ndarray): - filtered[k] = [filter_scalar_nan(v) for v in value] - else: - filtered[k] = filter_scalar_nan(value) - - return filtered - - -def parse_name(layername: str) -> Tuple[str, str, str]: - """ - Based on the layer name find out: - - * whether it's a timml or ttim element; - * which element type it is; - * what the user provided name is. - - For example: - parse_name("timml Headwell: drainage") -> ("timml", "Head Well", "drainage") - - Some grouping of tables occurs here. - """ - prefix, name = layername.split(":") - element_type = re.split("timml |ttim ", prefix)[1] - mapping = { - "Computation Times": "Domain", - "Temporal Settings": "Aquifer", - "Polygon Inhomogeneity Properties": "Polygon Inhomogeneity", - "Building Pit Properties": "Building Pit", - } - element_type = mapping.get(element_type, element_type) - if "timml" in prefix: - if "Properties" in prefix: - tim_type = "timml_assoc" - else: - tim_type = "timml" - elif "ttim" in prefix: - tim_type = "ttim" - else: - raise ValueError(f"Neither timml nor ttim in layername: {layername}") - return tim_type, element_type, name - - -def model_specification( - path: Union[str, pathlib.Path], active_elements: Dict[str, bool] -) -> Tuple[TimmlModelSpecification, TtimModelSpecification]: - """ - Group the different layers of a GPKG into model specifications for timml - and ttim. The grouping occurs solely on the basis of layer names. - """ - # Start by listing all layers - gpkg_names = fiona.listlayers(path) - # Group all these names together (using a defaultdict) - dd = defaultdict - grouped_names = dd(partial(dd, partial(dd, list))) - for layername in gpkg_names: - tim_type, element_type, name = parse_name(layername) - grouped_names[element_type][name][tim_type] = layername - - # Grab the names of the required elements and load the data into - # geodataframes. - aquifer_entry = grouped_names.pop("Aquifer")["Aquifer"] - aquifer = gpd.read_file(path, layer=aquifer_entry["timml"]) - temporal_settings = gpd.read_file(path, layer=aquifer_entry["ttim"]) - domain_entry = grouped_names.pop("Domain")["Domain"] - domain = gpd.read_file(path, layer=domain_entry["timml"]) - - if len(domain.index) == 0: - raise ValueError("Domain not defined") - - output_times = gpd.read_file(path, layer=domain_entry["ttim"]) - - # Load the data all other elements into geodataframes. - ttim_elements = {} - timml_elements = {} - for element_type, element_group in grouped_names.items(): - for name, group in element_group.items(): - timml_name = group["timml"] - timml_assoc_name = group.get("timml_assoc", None) - ttim_name = group.get("ttim", timml_name) - - timml_df = gpd.read_file(path, layer=timml_name) - timml_assoc_df = ( - gpd.read_file(path, layer=timml_assoc_name) - if timml_assoc_name is not None - else None - ) - ttim_df = ( - gpd.read_file(path, layer=ttim_name) if ttim_name is not None else None - ) - - # Use the background aquifer for a semi-confined top - if element_type in ("Polygon Area Sink", "Polygon Semi-Confined Top"): - timml_assoc_df = aquifer - - timml_spec = ElementSpecification( - elementtype=element_type, - active=active_elements.get(timml_name, False), - dataframe=timml_df, - associated_dataframe=timml_assoc_df, - ) - - ttim_spec = TransientElementSpecification( - elementtype=element_type, - active=active_elements.get(ttim_name, False), - dataframe=ttim_df, - steady_spec=timml_spec, - ) - timml_elements[timml_name] = timml_spec - ttim_elements[ttim_name] = ttim_spec - - return ( - TimmlModelSpecification(aquifer, timml_elements, domain), - TtimModelSpecification( - aquifer, - temporal_settings, - ttim_elements, - domain, - np.sort(output_times["time"].values), - ), - ) - - -# Three helpers for conversion to Python scripts -# ---------------------------------------------- -def dict_to_kwargs_code(data: dict) -> str: - strings = [] - for key, value in data.items(): - if isinstance(value, np.ndarray): - value = value.tolist() - elif isinstance(value, str) and key not in ("model", "timmlmodel"): - value = f'"{value}"' - strings.append(f"{key}={value}") - return ",".join(strings) - - -def sanitized(name: str) -> str: - return name.split(":")[-1].replace(" ", "_") - - -def headgrid_code(domain: gpd.GeoDataFrame) -> Tuple[str, str]: - xmin, ymin, xmax, ymax = domain.bounds.iloc[0] - dy = (ymax - ymin) / 50.0 - if dy > 500.0: - dy = round(dy / 500.0) * 500.0 - elif dy > 50.0: - dy = round(dy / 50.0) * 50.0 - elif dy > 5.0: # round to five - dy = round(dy / 5.0) * 5.0 - elif dy > 1.0: - dy = round(dy) - (xmin, xmax, ymin, ymax) = round_extent((xmin, xmax, ymin, ymax), dy) - xmin += 0.5 * dy - xmax += 0.5 * dy - ymax -= 0.5 * dy - xmin -= 0.5 * dy - xg = f"np.arange({xmin}, {xmax}, {dy})" - yg = f"np.arange({ymax}, {ymin}, -{dy})" - return xg, yg - - -# Output methods -# -------------- -def round_extent(extent: Tuple[float], cellsize: float) -> Tuple[float]: - """ - Increases the extent until all sides lie on a coordinate - divisible by cellsize. - - Parameters - ---------- - extent: Tuple[float] - xmin, xmax, ymin, ymax - cellsize: float - Desired cell size of the output head grids - - Returns - ------- - extent: Tuple[float] - xmin, xmax, ymin, ymax - """ - xmin, xmax, ymin, ymax = extent - xmin = np.floor(xmin / cellsize) * cellsize - ymin = np.floor(ymin / cellsize) * cellsize - xmax = np.ceil(xmax / cellsize) * cellsize - ymax = np.ceil(ymax / cellsize) * cellsize - return xmin, xmax, ymin, ymax - - -def gridspec( - path: Union[pathlib.Path, str], cellsize: float -) -> Tuple[Tuple[float], Any]: - """ - Infer the grid specification from the geopackage ``timmlDomain`` layer and - the provided cellsize. - - Parameters - ---------- - path: Union[pathlib.Path, str] - Path to the GeoPackage file. - cellsize: float - Desired cell size of the output head grids - - Returns - ------- - extent: tuple[float] - xmin, xmax, ymin, ymax - crs: Any - Coordinate Reference System - """ - domain = gpd.read_file(path, layer="timml Domain:Domain") - xmin, ymin, xmax, ymax = domain.bounds.iloc[0] - extent = (xmin, xmax, ymin, ymax) - return round_extent(extent, cellsize), domain.crs diff --git a/gistim/compute.py b/gistim/compute.py index 6f87ca3..9f12275 100644 --- a/gistim/compute.py +++ b/gistim/compute.py @@ -1,120 +1,330 @@ +import json import pathlib -from typing import Dict, Union +from collections import defaultdict +from functools import singledispatch +from typing import Any, Dict, List, Union -import geopandas as gpd import numpy as np -import rioxarray # noqa # pylint: disable=unused-import +import pandas as pd +import timml +import ttim import xarray as xr -import gistim +from gistim.geopackage import CoordinateReferenceSystem, write_geopackage +from gistim.netcdf import write_raster, write_ugrid +TIMML_MAPPING = { + "Constant": timml.Constant, + "Uflow": timml.Uflow, + "CircAreaSink": timml.CircAreaSink, + "Well": timml.Well, + "HeadWell": timml.HeadWell, + "PolygonInhomMaq": timml.PolygonInhomMaq, + "HeadLineSinkString": timml.HeadLineSinkString, + "LineSinkDitchString": timml.LineSinkDitchString, + "LeakyLineDoubletString": timml.LeakyLineDoubletString, + "ImpLineDoubletString": timml.ImpLineDoubletString, + "BuildingPit": timml.BuildingPit, + "LeakyBuildingPit": timml.LeakyBuildingPit, +} +TTIM_MAPPING = { + "CircAreaSink": ttim.CircAreaSink, + "Well": ttim.Well, + "HeadWell": ttim.HeadWell, + "HeadLineSinkString": ttim.HeadLineSinkString, + "LineSinkDitchString": ttim.LineSinkDitchString, + "LeakyLineDoubletString": ttim.LeakyLineDoubletString, +} -def write_raster( - head: xr.DataArray, - crs: int, - outpath: Union[pathlib.Path, str], -) -> None: - out = head.rio.write_crs(crs).rio.write_coordinate_system() - out.to_netcdf(outpath.with_suffix(".nc")) - return +def initialize_elements(model, mapping, data): + elements = defaultdict(list) + for name, entry in data.items(): + klass = mapping[entry["type"]] + for kwargs in entry["data"]: + element = klass(model=model, **kwargs) + elements[name].append(element) + return elements -def write_ugrid( - head: xr.DataArray, - crs: int, - outpath: Union[pathlib.Path, str], -) -> None: - ugrid_head = gistim.to_ugrid2d(head) - ugrid_head["projected_coordinate_system"] = xr.DataArray( - data=np.int32(0), - attrs={"epsg": np.int32(crs.to_epsg())}, + +def initialize_timml(data): + aquifer = data.pop("ModelMaq") + timml_model = timml.ModelMaq(**aquifer) + elements = initialize_elements(timml_model, TIMML_MAPPING, data) + return timml_model, elements + + +def initialize_ttim(data, timml_model): + aquifer = data.pop("ModelMaq") + ttim_model = ttim.ModelMaq(**aquifer, timmlmodel=timml_model) + elements = initialize_elements(ttim_model, TTIM_MAPPING, data) + return ttim_model, elements + + +@singledispatch +def headgrid(model, **kwargs): + raise TypeError("Expected timml or ttim model") + + +@headgrid.register +def _( + model: timml.Model, + xmin: float, + xmax: float, + ymin: float, + ymax: float, + spacing: float, + **_, +) -> xr.DataArray: + """ + Compute the headgrid of the TimML model, and store the results + in an xarray DataArray with the appropriate dimensions. + + Parameters + ---------- + model: timml.Model + Solved model to get heads from + data: Dict[str, Any] + + Returns + ------- + head: xr.DataArray + DataArray with dimensions ``("layer", "y", "x")``. + """ + x = np.arange(xmin, xmax, spacing) + 0.5 * spacing + # In geospatial rasters, y is DECREASING with row number + y = np.arange(ymax, ymin, -spacing) - 0.5 * spacing + head = model.headgrid(xg=x, yg=y) + nlayer = model.aq.find_aquifer_data(x[0], y[0]).naq + layer = [i for i in range(nlayer)] + return xr.DataArray( + data=head, + name="head", + coords={"layer": layer, "y": y, "x": x}, + dims=("layer", "y", "x"), ) - ugrid_head.to_netcdf(outpath.with_suffix(".ugrid.nc")) - return -def write_vector( - gdf_head: gpd.GeoDataFrame, - crs: int, - outpath: Union[pathlib.Path, str], - layername: str, -) -> None: - if len(gdf_head.index) > 0: - gdf_head = gdf_head.set_crs(crs) - gdf_head.to_file( - outpath.with_suffix(".output.gpkg"), - driver="GPKG", - layer=layername, +@headgrid.register +def _( + model: ttim.ModelMaq, + xmin: float, + xmax: float, + ymin: float, + ymax: float, + spacing: float, + reference_date: str, + time: List[float], +) -> Union[None, xr.DataArray]: + if time is None: + return None + + # Get coordinates ready + x = np.arange(xmin, xmax, spacing) + 0.5 * spacing + # In geospatial rasters, y is DECREASING with row number + y = np.arange(ymax, ymin, -spacing) - 0.5 * spacing + nlayer = model.aq.find_aquifer_data(x[0], y[0]).naq + + if 0.0 in time: + steady_head = model.timmlmodel.headgrid(xg=x, yg=y)[:, np.newaxis, :, :] + transient_head = model.headgrid(xg=x, yg=y, t=time[1:]) + head = np.hstack((steady_head, transient_head)) + else: + head = model.headgrid(xg=x, yg=y, t=time) + + # Other coordinates + layer = [i for i in range(nlayer)] + time = pd.to_datetime(reference_date) + pd.to_timedelta(time, "D") + return xr.DataArray( + data=head, + name="head", + coords={"layer": layer, "time": time, "y": y, "x": x}, + dims=("layer", "time", "y", "x"), + ) + + +@singledispatch +def head_observations(model, observations): + raise TypeError("Expected timml or ttim model") + + +@head_observations.register +def _( + model: timml.Model, + observations: Dict, + **_, +) -> Dict[str, pd.DataFrame]: + d = {"geometry": [], "label": []} + heads = [] + for kwargs in observations: + x = kwargs["x"] + y = kwargs["y"] + heads.append(model.head(x=x, y=y)) + d["geometry"].append({"type": "Point", "coordinates": [x, y]}) + d["label"].append(kwargs["label"]) + for i, layerhead in enumerate(np.vstack(heads).T): + d[f"head_layer{i}"] = layerhead + return pd.DataFrame(d) + + +@head_observations.register +def _( + model: ttim.ModelMaq, observations: Dict, reference_date: pd.Timestamp +) -> Dict[str, pd.DataFrame]: + d = { + "geometry": [], + "datetime_start": [], + "datetime_end": [], + "label": [], + "observation_id": [], + } + heads = [] + reference_date = pd.to_datetime(reference_date, utc=False) + for observation_id, kwargs in enumerate(observations): + x = kwargs["x"] + y = kwargs["y"] + t = kwargs["t"] + n_time = len(t) + datetime = reference_date + pd.to_timedelta([0] + t, "day") + d["geometry"].extend([{"type": "Point", "coordinates": [x, y]}] * n_time) + d["datetime_start"].extend(datetime[:-1]) + d["datetime_end"].extend(datetime[1:] - pd.to_timedelta(1, "minute")) + d["label"].extend([kwargs["label"]] * n_time) + d["observation_id"].extend([observation_id] * n_time) + heads.append(model.head(x=x, y=y, t=t)) + + for i, layerhead in enumerate(np.hstack(heads)): + d[f"head_layer{i}"] = layerhead + + df = pd.DataFrame(d) + return df + + +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 write_output( + model: Union[timml.Model, ttim.ModelMaq], + elements: Dict[str, Any], + data: Dict[str, Any], + path: Union[pathlib.Path, str], +): + crs = CoordinateReferenceSystem(**data["crs"]) + output_options = data["output_options"] + observations = data["observations"] + reference_date = pd.to_datetime(data.get("reference_date")) + + # Compute gridded head data and write to netCDF. + head = None + if output_options["raster"] or output_options["mesh"]: + head = headgrid(model, **data["headgrid"], reference_date=reference_date) + + if head is not None: + if output_options["raster"]: + write_raster(head, crs, path) + if output_options["mesh"]: + write_ugrid(head, crs, path) + + # Compute observations and discharge, and write to geopackage. + if output_options["discharge"]: + tables = extract_discharges( + elements, model.aq.nlayers, reference_date=reference_date ) + else: + tables = {} + + if output_options["head_observations"] and observations: + for layername, content in observations.items(): + tables[layername] = head_observations( + model, content["data"], reference_date=reference_date + ) + + write_geopackage(tables, crs, path) return def compute_steady( - inpath: Union[pathlib.Path, str], - outpath: Union[pathlib.Path, str], - cellsize: float, - active_elements: Dict[str, bool], + path: Union[pathlib.Path, str], ) -> None: - inpath = pathlib.Path(inpath) - outpath = pathlib.Path(outpath) + with open(path, "r") as f: + data = json.load(f) - timml_spec, _ = gistim.model_specification(inpath, active_elements) - timml_model, _, observations = gistim.timml_elements.initialize_model(timml_spec) + timml_model, elements = initialize_timml(data["timml"]) timml_model.solve() - extent, crs = gistim.gridspec(inpath, cellsize) - gdfs_obs = gistim.timml_elements.head_observations(timml_model, observations) - head = gistim.timml_elements.headgrid(timml_model, extent, cellsize) - - write_raster(head, crs, outpath) - write_ugrid(head, crs, outpath) - - for name, gdf in gdfs_obs.items(): - write_vector(gdf, crs, outpath, layername=name) - + write_output( + timml_model, + elements, + data, + path, + ) return def compute_transient( - inpath: Union[pathlib.Path, str], - outpath: Union[pathlib.Path, str], - cellsize: float, - active_elements: Dict[str, bool], + path: Union[pathlib.Path, str], ) -> None: - inpath = pathlib.Path(inpath) - outpath = pathlib.Path(outpath) + with open(path, "r") as f: + data = json.load(f) - timml_spec, ttim_spec = gistim.model_specification(inpath, active_elements) - timml_model, _, observations = gistim.timml_elements.initialize_model(timml_spec) - ttim_model, _, observations = gistim.ttim_elements.initialize_model( - ttim_spec, timml_model - ) + timml_model, _ = initialize_timml(data["timml"]) + ttim_model, ttim_elements = initialize_ttim(data["ttim"], timml_model) timml_model.solve() ttim_model.solve() - extent, crs = gistim.gridspec(inpath, cellsize) - refdate = ttim_spec.temporal_settings["reference_date"].iloc[0] - gdfs_obs = gistim.ttim_elements.head_observations(ttim_model, refdate, observations) - - # If no output times are specified, just compute the TimML steady-state - # heads. - if len(ttim_spec.output_times) > 0: - head = gistim.ttim_elements.headgrid( - ttim_model, - extent, - cellsize, - ttim_spec.output_times, - refdate, - ) - write_raster(head, crs, outpath) - write_ugrid(head, crs, outpath) - else: - head = gistim.timml_elements.headgrid(timml_model, extent, cellsize) - - write_raster(head, crs, outpath) - write_ugrid(head, crs, outpath) + write_output( + ttim_model, + ttim_elements, + data, + path, + ) + return - for name, gdf in gdfs_obs.items(): - write_vector(gdf, crs, outpath, layername=name) +def compute(path: str, transient: bool): + path = pathlib.Path(path) + if not path.exists(): + raise FileExistsError(path) + if transient: + compute_transient(path) + else: + compute_steady(path) return diff --git a/gistim/geomet/__init__.py b/gistim/geomet/__init__.py new file mode 100644 index 0000000..9a13841 --- /dev/null +++ b/gistim/geomet/__init__.py @@ -0,0 +1,17 @@ +# These modules are copied from: https://github.com/geomet/geomet +# Which has the following copyright notice: +# +# +# Copyright 2013 Lars Butler & individual contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. diff --git a/gistim/geomet/geopackage.py b/gistim/geomet/geopackage.py new file mode 100644 index 0000000..e46caab --- /dev/null +++ b/gistim/geomet/geopackage.py @@ -0,0 +1,410 @@ +# Copyright 2020 Tom Caruso & individual contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import struct as _struct + +from gistim.geomet import wkb as _wkb +from gistim.geomet.util import as_bin_str as _as_bin_str +from gistim.geomet.util import endian_token as _endian_token +from gistim.geomet.util import take as _take + + +def dump(obj, dest_file, big_endian=True): + """ + Dump GeoJSON-like `dict` to GeoPackage binary + and write it to the `dest_file`. + + :param dict obj: + A GeoJSON-like dictionary. It must at least the keys 'type' and + 'coordinates'. + :param dest_file: + Open and writable file-like object. + :param bool big_endian: + specify endianess of the dumped object. + + :return: + """ + dest_file.write(dumps(obj, big_endian)) + + +def load(source_file): + """ + Load a GeoJSON `dict` object from a ``source_file`` containing + GeoPackage (as a byte string). + + :param source_file: + Open and readable file-like object. + + :return: + A GeoJSON `dict` representing the geometry read from the file. + """ + return loads(source_file.read()) + + +def dumps(obj, big_endian=True): + """ + Dump a GeoJSON-like dict to a GeoPackage bytestring. + + + If the dict contains a top-level 'meta' key like so: + + ``` + 'meta': {'srid': 4326} + ``` + then the srid will be added to the geopackage header, but *not* + to the WKB geometry header. + + + If the dict contains a top-level 'bbox' key like so: + + ``` + 'bbox': [0, 0, 3, 3] + ``` + + Then an envelope will be added to the geopackage header + with this information. + + + If the geometry's coordinates are empty (an empty list) + then the geopackage header's "empty" flag will be set, + denoting that this geometry has no coordinates. + + Please note that while this library can parse geopackages + with a mixed byte-order in the header, it will only produce + blobs with consistent byte order (albeit properly marked as such). + That means you cannot product a geopackage with e.g. little-endian + header and big-endian WKB geometry. + + :param dict obj: + The geojson geometry to dump + :param bool big_endian: + if True, the geopackage binary will use big-endian + byte order, little-endian otherwise. + + :return bytes: + bytestring representing the geometry in geopackage + format. + """ + header = _build_geopackage_header(obj, not big_endian) + result = _wkb._dumps(obj, big_endian, include_meta=False) + return header + result + + +def loads(string): + """ + Construct a GeoJSON `dict` from geopackage (string). + + This function strips the geopackage header from the + string and passes the remaining WKB geometry to the + `geomet.wkb.loads` function. + + The envelope, if present, is added to the GeoJSON as + a key called 'bbox' as per the GeoJSON spec, [1]. + + If an SRID is specified in the geopackage header + AND the wkb header, the SRID in the geopackage header + will take precedence and will replace that SRID + in the returned dict. + + [1] https://tools.ietf.org/html/rfc7946#section-5 + + :param bytes string: + geopackage byte string. + :return dict: + GeoJSON represented the parsed geopackage binary. + """ + string = iter(string) + + header = _as_bin_str(_take(_GeoPackage.HEADER_LEN, string)) + + _check_is_valid(header) + g, p, version, empty, envelope_indicator, is_little_endian, srid = _parse_header( + header + ) + + wkb_offset = _get_wkb_offset(envelope_indicator) + left_to_take = wkb_offset - _GeoPackage.HEADER_LEN + envelope_data = _as_bin_str(_take(left_to_take, string)) + + if envelope_data: + envelope = _parse_envelope(envelope_indicator, envelope_data, is_little_endian) + + result = _wkb.loads(string) + + if srid: + result["meta"] = {"srid": int(srid)} + result["crs"] = { + "type": "name", + "properties": {"name": "EPSG%s" % srid}, + } + + if envelope_data: + result["bbox"] = envelope + + return result + + +class _GeoPackage: + """ + Much more information on geopackage structure + can be found here: http://www.geopackage.org/spec/#gpb_format + """ + + # The ascii letter 'G' + MAGIC1 = 0x47 + # The ascii letter 'P' + MAGIC2 = 0x50 + VERSION1 = 0x00 + HEADER_LEN = 8 + HEADER_PACK_FMT = "BBBBI" + ENVELOPE_2D_LEN = 32 + ENVELOPE_3D_LEN = 48 + ENVELOPE_4D_LEN = 64 + ENVELOPE_MASK = 0b00001111 + EMPTY_GEOM_MASK = 0b00011111 + ENDIANNESS_MASK = 0b00000001 + + +# map the "envelope indicator" integer we get out of the geopackage header +# to the dimensionality of the envelope. +# more info here: http://www.geopackage.org/spec/#gpb_format +# in the "flags" section, bits 3, 2, 1. +_indicator_to_dim = { + 0: 0, + 1: 4, + 2: 6, + 3: 6, + 4: 8, +} + +# Map the dimensionality of our envelope to the indicator +# integer we will use in the geopackage binary header. +# because we have no way to tell between Z and M values, +# if the geometry has 3 dimensions we default to assume Z. +_dim_to_indicator = {0: 0, 4: 1, 6: 2, 8: 4} + + +def is_valid(data): + """ + Check if the data represents a valid geopackage + geometry. Input can be either the full geometry or + just the header. + + :param bytes data: + bytes representing the geopackage binary. + + :return (bool, str): + Is the geopackage valid, if not, string describing why + """ + g, p, version, _, envelope_indicator, _, _ = _parse_header(data[:8]) + if (g != _GeoPackage.MAGIC1) or (p != _GeoPackage.MAGIC2): + return False, "Missing Geopackage header magic bytes" + if version != _GeoPackage.VERSION1: + return False, "Geopackage version must be 0" + if (envelope_indicator < 0) or (envelope_indicator > 4): + return False, "Envelope indicator must be between 0-4" + return True, "" + + +def _header_is_little_endian(header): + """ + Check to see if the header is encoded + as little endian or big endian. + + Either the entire binary blob or + just the header can be passed in. + + :param bytes header: + geopackage header or binary blob + + :return bool: is the header little endian + """ + (flags,) = _struct.unpack("B", header[3:4]) + return flags & _GeoPackage.ENDIANNESS_MASK + + +def _parse_header(header): + """ + Unpack all information from the geopackage + header, including "magic" GP bytes. Returns + all of them so we can confirm that this + geopackage is validly formed. Can also accept + the full binary blob. + + :param header: + the header or the full geometry. + + :return 7-tuple: + all attributes stored in the binary header. + """ + is_little_endian = _header_is_little_endian(header) + fmt = _endian_token(is_little_endian) + _GeoPackage.HEADER_PACK_FMT + + g, p, version, flags, srid = _struct.unpack(fmt, header[: _GeoPackage.HEADER_LEN]) + empty, envelope_indicator, endianness = _parse_flags(flags) + return g, p, version, empty, envelope_indicator, endianness, srid + + +def _parse_flags(flags): + """ + Parse the bits in the "flags" byte + of the geopackage header to retrieve + useful information. We specifically parse + the endianness, the envelope indicator, + and the "empty" flag. + + Much more info can be found in + the documentation [1]. + + [1] http://www.geopackage.org/spec/#gpb_format + :param byte flags: + The "flags" byte of a geopackage header. + + :return tuple: + """ + endianness = flags & _GeoPackage.ENDIANNESS_MASK + envelope_indicator = (flags & _GeoPackage.ENVELOPE_MASK) >> 1 + empty = (flags & _GeoPackage.EMPTY_GEOM_MASK) >> 4 + return empty, envelope_indicator, endianness + + +def _build_flags(empty, envelope_indicator, is_little_endian=1): + """ + Create the "flags" byte which goes into + the geopackage header. Much more info + can be found in the documentation [1]. + + [1] http://www.geopackage.org/spec/#gpb_format + + :param int empty: + 0 or 1 indicating whether the geometry is empty. + True and False also work as expected. + :param int envelope_indicator: + indicates the dimensionality of the envelope. + :param int is_little_endian: + 0 or 1 (or False / True) indicating + whether the header should be + little-endian encoded. + + :return byte: + geopackage header flags + """ + flags = 0b0 + if empty: + flags = (flags | 1) << 3 + if envelope_indicator: + flags = flags | envelope_indicator + + return (flags << 1) | is_little_endian + + +def _build_geopackage_header(obj, is_little_endian): + """ + Create the geopackage header for the input object. + Looks for a 'bbox' key on the geometry to use + for an envelope, and a 'meta' key with an + SRID to encode into the header. + + :param dict obj: + a geojson object + :param bool is_little_endian: + which endianness to use when + encoding the data. + + :return bytes: geopackage header. + """ + # Collect geometry metadata. + empty = 1 if len(obj["coordinates"]) == 0 else 0 + envelope = obj.get("bbox", []) + srid = obj.get("meta", {}).get("srid", 0) + + try: + envelope_indicator = _dim_to_indicator[len(envelope)] + except KeyError: + raise ValueError( + "Bounding box must be of length 2*n where " + "n is the number of dimensions represented " + "in the contained geometries." + ) + + pack_args = [ + _GeoPackage.MAGIC1, + _GeoPackage.MAGIC2, + _GeoPackage.VERSION1, + # This looks funny, but _build_flags wants a 1 or 0 for + # "little endian" because it uses it to `or` with the bits. + # Conveniently, in Python, False == 0 and True == 1, so + # we can pass the boolean right in and it works as expected. + _build_flags(empty, envelope_indicator, is_little_endian), + srid, + ] + + pack_fmt = _endian_token(is_little_endian) + _GeoPackage.HEADER_PACK_FMT + + # This has no effect if we have a 0 envelope indicator. + pack_fmt += "d" * _indicator_to_dim[envelope_indicator] + pack_args.extend(envelope) + + return _struct.pack(pack_fmt, *pack_args) + + +def _check_is_valid(data): + """ + Raise if the header is not valid geopackage. + + :param bytes data: Geopackage data or header. + + :return None: + """ + valid, reason = is_valid(data) + if not valid: + raise ValueError( + "Could not read Geopackage geometry " "because of errors: " + reason + ) + + +def _get_wkb_offset(envelope_indicator): + """ + Get the full byte offset at which the WKB geometry lies + in the geopackage geometry. + + :param int envelope_indicator: + indicates the dimensionality of the envelope. + + :return int: + number of bytes until the beginning of the + WKB geometry. + + """ + base_len = _GeoPackage.HEADER_LEN + return (base_len * _indicator_to_dim[envelope_indicator]) + base_len + + +def _parse_envelope(envelope_indicator, envelope, is_little_endian): + """ + Parse a geopackage envelope bytestring into an n-tuple + of floats. + + :param int envelope_indicator: + indicates the dimensionality of the envelope. + :param bytes envelope: + Bytestring of the envelope values. + :param bool is_little_endian: + how to pack the bytes in the envelope. + + :return tuple[float]: Geometry envelope. + """ + pack_fmt = _endian_token(is_little_endian) + pack_fmt += "d" * _indicator_to_dim[envelope_indicator] + return _struct.unpack(pack_fmt, envelope) diff --git a/gistim/geomet/util.py b/gistim/geomet/util.py new file mode 100644 index 0000000..804b112 --- /dev/null +++ b/gistim/geomet/util.py @@ -0,0 +1,127 @@ +# Copyright 2013 Lars Butler & individual contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import collections.abc as collections +import itertools + + +def block_splitter(data, block_size): + """ + Creates a generator by slicing ``data`` into chunks of ``block_size``. + + >>> data = range(10) + >>> list(block_splitter(data, 2)) + [[0, 1], [2, 3], [4, 5], [6, 7], [8, 9]] + + If ``data`` cannot be evenly divided by ``block_size``, the last block will + simply be the remainder of the data. Example: + + >>> data = range(10) + >>> list(block_splitter(data, 3)) + [[0, 1, 2], [3, 4, 5], [6, 7, 8], [9]] + + If the ``block_size`` is greater than the total length of ``data``, a + single block will be generated: + + >>> data = range(3) + >>> list(block_splitter(data, 4)) + [[0, 1, 2]] + + :param data: + Any iterable. If ``data`` is a generator, it will be exhausted, + obviously. + :param int block_site: + Desired (maximum) block size. + """ + buf = [] + for i, datum in enumerate(data): + buf.append(datum) + if len(buf) == block_size: + yield buf + buf = [] + + # If there's anything leftover (a partial block), + # yield it as well. + if buf: + yield buf + + +def take(n, iterable): + """ + Return first n items of the iterable as a list + + Copied shamelessly from + http://docs.python.org/2/library/itertools.html#recipes. + """ + return list(itertools.islice(iterable, n)) + + +def as_bin_str(a_list): + return bytes(a_list) + + +def round_geom(geom, precision=None): + """Round coordinates of a geometric object to given precision.""" + if geom["type"] == "Point": + x, y = geom["coordinates"] + xp, yp = [x], [y] + if precision is not None: + xp = [round(v, precision) for v in xp] + yp = [round(v, precision) for v in yp] + new_coords = tuple(zip(xp, yp))[0] + if geom["type"] in ["LineString", "MultiPoint"]: + xp, yp = zip(*geom["coordinates"]) + if precision is not None: + xp = [round(v, precision) for v in xp] + yp = [round(v, precision) for v in yp] + new_coords = tuple(zip(xp, yp)) + elif geom["type"] in ["Polygon", "MultiLineString"]: + new_coords = [] + for piece in geom["coordinates"]: + xp, yp = zip(*piece) + if precision is not None: + xp = [round(v, precision) for v in xp] + yp = [round(v, precision) for v in yp] + new_coords.append(tuple(zip(xp, yp))) + elif geom["type"] == "MultiPolygon": + parts = geom["coordinates"] + new_coords = [] + for part in parts: + inner_coords = [] + for ring in part: + xp, yp = zip(*ring) + if precision is not None: + xp = [round(v, precision) for v in xp] + yp = [round(v, precision) for v in yp] + inner_coords.append(tuple(zip(xp, yp))) + new_coords.append(inner_coords) + return {"type": geom["type"], "coordinates": new_coords} + + +def flatten_multi_dim(sequence): + """Flatten a multi-dimensional array-like to a single dimensional sequence + (as a generator). + """ + for x in sequence: + if isinstance(x, collections.Iterable) and not isinstance(x, str): + for y in flatten_multi_dim(x): + yield y + else: + yield x + + +def endian_token(is_little_endian): + if is_little_endian: + return "<" + else: + return ">" diff --git a/gistim/geomet/wkb.py b/gistim/geomet/wkb.py new file mode 100644 index 0000000..bb848f9 --- /dev/null +++ b/gistim/geomet/wkb.py @@ -0,0 +1,936 @@ +# Copyright 2013 Lars Butler & individual contributors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import binascii +import struct +from itertools import chain + +from gistim.geomet.util import as_bin_str, block_splitter, flatten_multi_dim, take + +#: '\x00': The first byte of any WKB string. Indicates big endian byte +#: ordering for the data. +BIG_ENDIAN = b"\x00" +#: '\x01': The first byte of any WKB string. Indicates little endian byte +#: ordering for the data. +LITTLE_ENDIAN = b"\x01" +#: High byte in a 4-byte geometry type field to indicate that a 4-byte SRID +#: field follows. +SRID_FLAG = b"\x20" + +#: Mapping of GeoJSON geometry types to the "2D" 4-byte binary string +#: representation for WKB. "2D" indicates that the geometry is 2-dimensional, +#: X and Y components. +#: NOTE: Byte ordering is big endian. +WKB_2D = { + "Point": b"\x00\x00\x00\x01", + "LineString": b"\x00\x00\x00\x02", + "Polygon": b"\x00\x00\x00\x03", + "MultiPoint": b"\x00\x00\x00\x04", + "MultiLineString": b"\x00\x00\x00\x05", + "MultiPolygon": b"\x00\x00\x00\x06", + "GeometryCollection": b"\x00\x00\x00\x07", +} + +#: Mapping of GeoJSON geometry types to the "Z" 4-byte binary string +#: representation for WKB. "Z" indicates that the geometry is 3-dimensional, +#: with X, Y, and Z components. +#: NOTE: Byte ordering is big endian. +WKB_Z = { + "Point": b"\x00\x00\x03\xe9", + "LineString": b"\x00\x00\x03\xea", + "Polygon": b"\x00\x00\x03\xeb", + "MultiPoint": b"\x00\x00\x03\xec", + "MultiLineString": b"\x00\x00\x03\xed", + "MultiPolygon": b"\x00\x00\x03\xee", + "GeometryCollection": b"\x00\x00\x03\xef", +} + +#: Mapping of GeoJSON geometry types to the "M" 4-byte binary string +#: representation for WKB. "M" indicates that the geometry is 2-dimensional, +#: with X, Y, and M ("Measure") components. +#: NOTE: Byte ordering is big endian. +WKB_M = { + "Point": b"\x00\x00\x07\xd1", + "LineString": b"\x00\x00\x07\xd2", + "Polygon": b"\x00\x00\x07\xd3", + "MultiPoint": b"\x00\x00\x07\xd4", + "MultiLineString": b"\x00\x00\x07\xd5", + "MultiPolygon": b"\x00\x00\x07\xd6", + "GeometryCollection": b"\x00\x00\x07\xd7", +} + +#: Mapping of GeoJSON geometry types to the "ZM" 4-byte binary string +#: representation for WKB. "ZM" indicates that the geometry is 4-dimensional, +#: with X, Y, Z, and M ("Measure") components. +#: NOTE: Byte ordering is big endian. +WKB_ZM = { + "Point": b"\x00\x00\x0b\xb9", + "LineString": b"\x00\x00\x0b\xba", + "Polygon": b"\x00\x00\x0b\xbb", + "MultiPoint": b"\x00\x00\x0b\xbc", + "MultiLineString": b"\x00\x00\x0b\xbd", + "MultiPolygon": b"\x00\x00\x0b\xbe", + "GeometryCollection": b"\x00\x00\x0b\xbf", +} + +#: Mapping of dimension types to maps of GeoJSON geometry type -> 4-byte binary +#: string representation for WKB. +_WKB = { + "2D": WKB_2D, + "Z": WKB_Z, + "M": WKB_M, + "ZM": WKB_ZM, +} + +#: Mapping from binary geometry type (as a 4-byte binary string) to GeoJSON +#: geometry type. +#: NOTE: Byte ordering is big endian. +_BINARY_TO_GEOM_TYPE = dict( + chain(*((reversed(x) for x in wkb_map.items()) for wkb_map in _WKB.values())) +) + +_INT_TO_DIM_LABEL = {2: "2D", 3: "Z", 4: "ZM"} + + +def _get_geom_type(type_bytes): + """Get the GeoJSON geometry type label from a WKB type byte string. + + :param type_bytes: + 4 byte string in big endian byte order containing a WKB type number. + It may also contain a "has SRID" flag in the high byte (the first type, + since this is big endian byte order), indicated as 0x20. If the SRID + flag is not set, the high byte will always be null (0x00). + :returns: + 3-tuple ofGeoJSON geometry type label, the bytes representing the + geometry type, and a separate "has SRID" flag. If the input + `type_bytes` contains an SRID flag, it will be removed. + + >>> # Z Point, with SRID flag + >>> _get_geom_type(b'\\x20\\x00\\x03\\xe9') == ( + ... 'Point', b'\\x00\\x00\\x03\\xe9', True) + True + + >>> # 2D MultiLineString, without SRID flag + >>> _get_geom_type(b'\\x00\\x00\\x00\\x05') == ( + ... 'MultiLineString', b'\\x00\\x00\\x00\\x05', False) + True + + """ + # slice off the high byte, which may contain the SRID flag + high_byte = bytes([type_bytes[0]]) + has_srid = high_byte == b"\x20" + if has_srid: + # replace the high byte with a null byte + type_bytes = as_bin_str(b"\x00" + type_bytes[1:]) + else: + type_bytes = as_bin_str(type_bytes) + + # look up the geometry type + geom_type = _BINARY_TO_GEOM_TYPE.get(type_bytes) + return geom_type, type_bytes, has_srid + + +def dump(obj, dest_file): + """ + Dump GeoJSON-like `dict` to WKB and write it to the `dest_file`. + + :param dict obj: + A GeoJSON-like dictionary. It must at least the keys 'type' and + 'coordinates'. + :param dest_file: + Open and writable file-like object. + """ + dest_file.write(dumps(obj)) + + +def load(source_file): + """ + Load a GeoJSON `dict` object from a ``source_file`` containing WKB (as a + byte string). + + :param source_file: + Open and readable file-like object. + + :returns: + A GeoJSON `dict` representing the geometry read from the file. + """ + return loads(source_file.read()) + + +def dumps(obj, big_endian=True): + """ + Dump a GeoJSON-like `dict` to a WKB string. + + .. note:: + The dimensions of the generated WKB will be inferred from the first + vertex in the GeoJSON `coordinates`. It will be assumed that all + vertices are uniform. There are 4 types: + + - 2D (X, Y): 2-dimensional geometry + - Z (X, Y, Z): 3-dimensional geometry + - M (X, Y, M): 2-dimensional geometry with a "Measure" + - ZM (X, Y, Z, M): 3-dimensional geometry with a "Measure" + + If the first vertex contains 2 values, we assume a 2D geometry. + If the first vertex contains 3 values, this is slightly ambiguous and + so the most common case is chosen: Z. + If the first vertex contains 4 values, we assume a ZM geometry. + + The WKT/WKB standards provide a way of differentiating normal (2D), Z, + M, and ZM geometries (http://en.wikipedia.org/wiki/Well-known_text), + but the GeoJSON spec does not. Therefore, for the sake of interface + simplicity, we assume that geometry that looks 3D contains XYZ + components, instead of XYM. + + If the coordinates list has no coordinate values (this includes nested + lists, for example, `[[[[],[]], []]]`, the geometry is considered to be + empty. Geometries, with the exception of points, have a reasonable + "empty" representation in WKB; however, without knowing the number of + coordinate values per vertex, the type is ambigious, and thus we don't + know if the geometry type is 2D, Z, M, or ZM. Therefore in this case + we expect a `ValueError` to be raised. + + :param dict obj: + GeoJson-like `dict` object. + :param bool big_endian: + Defaults to `True`. If `True`, data values in the generated WKB will + be represented using big endian byte order. Else, little endian. + + :returns: + A WKB binary string representing of the ``obj``. + """ + return _dumps(obj, big_endian) + + +def _dumps(obj, big_endian=True, include_meta=True): + """ + Basically perform the action of dumps, but with some extra flags for + behavior specifically needed by the geopackage...package. + """ + geom_type = obj["type"] + if include_meta: + meta = obj.get("meta", {}) + else: + meta = {} + + exporter = _dumps_registry.get(geom_type) + if exporter is None: + _unsupported_geom_type(geom_type) + + # Check for empty geometries. GeometryCollections have a slightly different + # JSON/dict structure, but that's handled. + coords_or_geoms = obj.get("coordinates", obj.get("geometries")) + if len(list(flatten_multi_dim(coords_or_geoms))) == 0: + raise ValueError( + "Empty geometries cannot be represented in WKB. Reason: The " + "dimensionality of the WKB would be ambiguous." + ) + + return exporter(obj, big_endian, meta) + + +def loads(string): + """ + Construct a GeoJSON `dict` from WKB (`string`). + + The resulting GeoJSON `dict` will include the SRID as an integer in the + `meta` object. This was an arbitrary decision made by `geomet, the + discussion of which took place here: + https://github.com/geomet/geomet/issues/28. + + In order to be consistent with other libraries [1] and (deprecated) + specifications [2], also include the same information in a `crs` + object. This isn't ideal, but the `crs` member is no longer part of + the GeoJSON standard, according to RFC7946 [3]. However, it's still + useful to include this information in GeoJSON payloads because it + supports conversion to EWKT/EWKB (which are canonical formats used by + PostGIS and the like). + + Example: + + {'type': 'Point', + 'coordinates': [0.0, 1.0], + 'meta': {'srid': 4326}, + 'crs': {'type': 'name', 'properties': {'name': 'EPSG4326'}}} + + NOTE(larsbutler): I'm not sure if it's valid to just prefix EPSG + (European Petroluem Survey Group) to an SRID like this, but we'll + stick with it for now until it becomes a problem. + + NOTE(larsbutler): Ideally, we should use URNs instead of this + notation, according to the new GeoJSON spec [4]. However, in + order to be consistent with [1], we'll stick with this approach + for now. + + References: + + [1] - https://github.com/bryanjos/geo/issues/76 + [2] - http://geojson.org/geojson-spec.html#coordinate-reference-system-objects + [3] - https://tools.ietf.org/html/rfc7946#appendix-B.1 + [4] - https://tools.ietf.org/html/rfc7946#section-4 + """ # noqa + string = iter(string) + # endianness = string[0:1] + endianness = as_bin_str(take(1, string)) + if endianness == BIG_ENDIAN: + big_endian = True + elif endianness == LITTLE_ENDIAN: + big_endian = False + else: + raise ValueError( + "Invalid endian byte: '0x%s'. Expected 0x00 or 0x01" + % binascii.hexlify(endianness.encode()).decode() + ) + + endian_token = ">" if big_endian else "<" + # type_bytes = string[1:5] + type_bytes = as_bin_str(take(4, string)) + if not big_endian: + # To identify the type, order the type bytes in big endian: + type_bytes = type_bytes[::-1] + + geom_type, type_bytes, has_srid = _get_geom_type(type_bytes) + srid = None + if has_srid: + srid_field = as_bin_str(take(4, string)) + [srid] = struct.unpack("%si" % endian_token, srid_field) + + # data_bytes = string[5:] # FIXME: This won't work for GeometryCollections + data_bytes = string + + importer = _loads_registry.get(geom_type) + + if importer is None: + _unsupported_geom_type(geom_type) + + data_bytes = iter(data_bytes) + result = importer(big_endian, type_bytes, data_bytes) + if has_srid: + # As mentioned in the docstring above, include both approaches to + # indicating the SRID. + result["meta"] = {"srid": int(srid)} + result["crs"] = { + "type": "name", + "properties": {"name": "EPSG%s" % srid}, + } + return result + + +def _unsupported_geom_type(geom_type): + raise ValueError("Unsupported geometry type '%s'" % geom_type) + + +# TODO: dont default meta to none +def _header_bytefmt_byteorder(geom_type, num_dims, big_endian, meta=None): + """ + Utility function to get the WKB header (endian byte + type header), byte + format string, and byte order string. + """ + dim = _INT_TO_DIM_LABEL.get(num_dims) + if dim is None: + pass # TODO: raise + + type_byte_str = _WKB[dim][geom_type] + srid = meta.get("srid") + if srid is not None: + # Add the srid flag + type_byte_str = SRID_FLAG + type_byte_str[1:] + + if big_endian: + header = BIG_ENDIAN + byte_fmt = b">" + byte_order = ">" + else: + header = LITTLE_ENDIAN + byte_fmt = b"<" + byte_order = "<" + # reverse the byte ordering for little endian + type_byte_str = type_byte_str[::-1] + + header += type_byte_str + if srid is not None: + srid = int(srid) + + if big_endian: + srid_header = struct.pack(">i", srid) + else: + srid_header = struct.pack(" pd.DataFrame: + # From records? + return pd.DataFrame( + data={ + "table_name": table_names, + "data_type": "features", + "identifier": table_names, + "description": "", + "last_change": pd.Timestamp.now(), + "min_x": [bb.xmin for bb in bounding_boxes], + "min_y": [bb.ymin for bb in bounding_boxes], + "max_x": [bb.xmax for bb in bounding_boxes], + "max_y": [bb.ymax for bb in bounding_boxes], + "srs_id": srs_id, + } + ) + + +def create_gkpg_spatial_ref_sys(crs: CoordinateReferenceSystem) -> pd.DataFrame: + return pd.DataFrame( + data={ + "srs_name": [ + "Undefined Cartesian SRS", + "Undefined geographic SRS", + "WGS 84 geodetic", + crs.description, + ], + "srs_id": [-1, 0, 4326, crs.srs_id], + "organization": ["NONE", "NONE", "EPSG", "EPSG"], + "organization_coordsys_id": [-1, 0, 4326, crs.organization], + "definition": ["undefined", "undefined", WGS84_WKT, crs.wkt], + "description": [ + "undefined Cartesian coordinate reference system", + "undefined geographic coordinate reference system", + "longitude/latitude coordinates in decimal degrees on the WGS 84 spheroid", + "", + ], + } + ) + + +def create_gpkg_geometry_columns( + table_names: List[str], + geometry_type_names: List[str], + srs_id: int, +) -> pd.DataFrame: + return pd.DataFrame( + data={ + "table_name": table_names, + "column_name": "geom", + "geometry_type_name": geometry_type_names, + "srs_id": srs_id, + "z": 0, + "m": 0, + } + ) + + +def collect_bounding_box(features, geometry_type) -> BoundingBox: + if geometry_type == "Point": + x = [] + y = [] + for point in features: + coordinates = point["coordinates"] + x.append(coordinates[0]) + y.append(coordinates[1]) + else: + x, y = zip( + *itertools.chain.from_iterable(line["coordinates"] for line in features) + ) + return BoundingBox(xmin=min(x), ymin=min(y), xmax=max(x), ymax=max(y)) + + +def process_table(dataframe: pd.DataFrame) -> Tuple[pd.DataFrame, BoundingBox, str]: + geometry = dataframe.pop("geometry").to_numpy() + geometry_type = set(f["type"] for f in geometry) + if len(geometry_type) != 1: + raise ValueError( + f"Table should contain exactly one geometry type. Received: {geometry_type}" + ) + # Get first (and only) member of set. + geometry_type = next(iter(geometry_type)) + bounding_box = collect_bounding_box(geometry, geometry_type) + dataframe["geom"] = [geopackage.dumps(f) for f in geometry] + return dataframe, bounding_box, geometry_type + + +def force_sql_datetime(df: pd.DataFrame): + """ + Pandas SQL writes datetimes as SQL Timestamps, which are not accepted by QGIS. + They have to SQL DATETIME instead. + """ + return { + colname: "DATETIME" + for colname, dtype in df.dtypes.to_dict().items() + if np.issubdtype(dtype, np.datetime64) + } + + +def write_geopackage( + tables: Dict[str, pd.DataFrame], crs: CoordinateReferenceSystem, path: Path +) -> None: + """ + Write the content of the tables to a geopackage. Overwrites the geopackage + if it already exists. + + Parameters + ---------- + tables: Dict[str, pd.DataFrame] + Mapping of table name to contents, as a dataframe. + May be an empty dictionary. + crs: CoordinateReferenceSystem + Coordinate reference system of the geometries. + path: Path + Path to the geopackage to write. + + Returns + ------- + None + """ + # We write all tables to a temporary GeoPackage with a dot prefix, and at + # the end move this over the target file. This does not throw a + # PermissionError if the file is open in QGIS. + gpkg_path = path.with_suffix(".output.gpkg") + temp_path = gpkg_path.with_stem(".") + # avoid adding tables to existing geopackage. + temp_path.unlink(missing_ok=True) + + try: + connection = sqlite3.connect(database=temp_path) + connection.execute(f"PRAGMA application_id = {APPLICATION_ID};") + connection.execute(f"PRAGMA user_version = {USER_VERSION};") + + table_names = [] + geometry_types = [] + bounding_boxes = [] + for layername, dataframe in tables.items(): + dataframe, bounding_box, geometry_type = process_table(dataframe) + table_names.append(layername) + geometry_types.append(geometry_type) + bounding_boxes.append(bounding_box) + dataframe.to_sql( + layername, con=connection, dtype=force_sql_datetime(dataframe) + ) + + # Create mandatory geopackage tables. + gpkg_contents = create_gpkg_contents( + table_names=table_names, bounding_boxes=bounding_boxes, srs_id=crs.srs_id + ) + gpkg_geometry_columns = create_gpkg_geometry_columns( + table_names=table_names, + geometry_type_names=geometry_types, + srs_id=crs.srs_id, + ) + gpkg_spatial_ref_sys = create_gkpg_spatial_ref_sys(crs) + # Write to Geopackage database. + gpkg_contents.to_sql(name="gpkg_contents", con=connection) + gpkg_geometry_columns.to_sql(name="gpkg_geometry_columns", con=connection) + gpkg_spatial_ref_sys.to_sql(name="gpkg_spatial_ref_sys", con=connection) + + finally: + connection.commit() + connection.close() + + shutil.move(temp_path, gpkg_path) + return diff --git a/gistim/netcdf.py b/gistim/netcdf.py new file mode 100644 index 0000000..fe9d9e5 --- /dev/null +++ b/gistim/netcdf.py @@ -0,0 +1,47 @@ +import pathlib +from typing import Union + +import numpy as np +import xarray as xr + +from gistim.geopackage import CoordinateReferenceSystem +from gistim.ugrid import to_ugrid2d + + +def write_raster( + head: xr.DataArray, + crs: CoordinateReferenceSystem, + outpath: Union[pathlib.Path, str], +) -> None: + # Write GDAL required metadata. + head["spatial_ref"] = xr.DataArray( + np.int32(0), attrs={"crs_wkt": crs.wkt, "spatial_ref": crs.wkt} + ) + head.attrs["grid_mapping"] = "spatial_ref" + head["x"].attrs = { + "axis": "X", + "long_name": "x coordinate of projection", + "standard_name": "projection_x_coordinate", + } + head["y"].attrs = { + "axis": "Y", + "long_name": "y coordinate of projection", + "standard_name": "projection_y_coordinate", + } + head.to_netcdf(outpath.with_suffix(".nc"), format="NETCDF3_CLASSIC") + return + + +def write_ugrid( + head: xr.DataArray, + crs: CoordinateReferenceSystem, + outpath: Union[pathlib.Path, str], +) -> None: + ugrid_head = to_ugrid2d(head) + # Write MDAL required metadata. + ugrid_head["projected_coordinate_system"] = xr.DataArray( + data=np.int32(0), + attrs={"epsg": np.int32(crs.srs_id)}, + ) + ugrid_head.to_netcdf(outpath.with_suffix(".ugrid.nc"), format="NETCDF3_CLASSIC") + return diff --git a/gistim/timml_elements.py b/gistim/timml_elements.py deleted file mode 100644 index 09c963a..0000000 --- a/gistim/timml_elements.py +++ /dev/null @@ -1,473 +0,0 @@ -""" -This model converts the geodataframe as read from the geopackage into keyword -arguments for TimML. - -These keyword arguments are used to initialize a model, or used to generate a -Python script. -""" -from collections import defaultdict -from typing import Any, Dict, List, Tuple - -import black -import geopandas as gpd -import numpy as np -import timml -import tqdm -import xarray as xr - -from gistim.common import ( - ElementSpecification, - TimmlModelSpecification, - aquifer_data, - dict_to_kwargs_code, - filter_nan, - headgrid_code, - linestring_coordinates, - point_coordinates, - polygon_coordinates, - sanitized, -) - - -def constant(spec: ElementSpecification) -> List[Dict[str, Any]]: - firstrow = filter_nan(spec.dataframe.iloc[0]) - x, y = point_coordinates(firstrow) - return [ - { - "xr": x, - "yr": y, - "hr": firstrow["head"], - "layer": firstrow["layer"], - "label": firstrow["label"], - } - ] - - -def uflow(spec: ElementSpecification) -> List[Dict[str, Any]]: - row = filter_nan(spec.dataframe.iloc[0]) - return [ - { - "slope": row["slope"], - "angle": row["angle"], - "label": row["label"], - } - ] - - -def observation(spec: ElementSpecification) -> List[Dict[str, Any]]: - dataframe = spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - kwargslist = [] - for (row, x, y) in zip(dataframe.to_dict("records"), X, Y): - row = filter_nan(row) - kwargslist.append( - { - "x": x, - "y": y, - "label": row["label"], - } - ) - return kwargslist - - -def well(spec: ElementSpecification) -> List[Dict[str, Any]]: - dataframe = spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - for (row, x, y) in zip(dataframe.to_dict("records"), X, Y): - row = filter_nan(row) - kwargslist.append( - { - "xw": x, - "yw": y, - "Qw": row["discharge"], - "rw": row["radius"], - "res": row["resistance"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def headwell(spec: ElementSpecification) -> List[Dict[str, Any]]: - dataframe = spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - for (row, x, y) in zip(dataframe.to_dict("records"), X, Y): - row = filter_nan(row) - kwargslist.append( - { - "xw": x, - "yw": y, - "hw": row["head"], - "rw": row["radius"], - "res": row["resistance"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def polygoninhom(spec: ElementSpecification) -> List[Dict[str, Any]]: - """ - Parameters - ---------- - dataframe: tuple of geopandas.GeoDataFrame - - Returns - ------- - None - """ - geometry = spec.dataframe - properties = spec.associated_dataframe.set_index("inhomogeneity_id") - # Iterate through the row containing the geometry - # and iterate through the associated table containing k properties. - kwarglist = [] - for row in geometry.to_dict("records"): - row = filter_nan(row) - dataframe = properties.loc[[row["inhomogeneity_id"]]] - kwargs = aquifer_data(dataframe) - kwargs["xy"] = polygon_coordinates(row) - kwargs["order"] = row["order"] - kwargs["ndeg"] = row["ndegrees"] - kwarglist.append(kwargs) - return kwarglist - - -def polygonareasink(spec: ElementSpecification) -> List[Dict[str, Any]]: - geometry = spec.dataframe - properties = spec.associated_dataframe.sort_values(by="layer").set_index("layer") - # Ensure there's no semiconfined top - dataframe = properties.copy() - dataframe.loc[0, "semiconf_head"] = np.nan - kwarglist = [] - for row in geometry.to_dict("records"): - row = filter_nan(row) - kwargs = aquifer_data(dataframe.reset_index()) - kwargs["N"] = row["rate"] - kwargs["xy"] = polygon_coordinates(row) - kwargs["order"] = row["order"] - kwargs["ndeg"] = row["ndegrees"] - kwarglist.append(kwargs) - return kwarglist - - -def polygontop(spec: ElementSpecification) -> List[Dict[str, Any]]: - geometry = spec.dataframe - properties = spec.associated_dataframe.sort_values(by="layer").set_index("layer") - kwarglist = [] - for row in geometry.to_dict("records"): - row = filter_nan(row) - dataframe = properties.copy() - dataframe.loc[0, "aquitard_c"] = row["aquitard_c"] - dataframe.loc[0, "semiconf_top"] = row["semiconf_top"] - dataframe.loc[0, "semiconf_head"] = row["semiconf_head"] - kwargs = aquifer_data(dataframe.reset_index()) - kwargs["xy"] = polygon_coordinates(row) - kwargs["order"] = row["order"] - kwargs["ndeg"] = row["ndegrees"] - kwarglist.append(kwargs) - return kwarglist - - -def buildingpit(spec: ElementSpecification) -> List[Dict[str, Any]]: - geometry = spec.dataframe - properties = spec.associated_dataframe.set_index("inhomogeneity_id") - # Iterate through the row containing the geometry - # and iterate through the associated table containing k properties. - kwarglist = [] - for row in geometry.to_dict("records"): - row = filter_nan(row) - dataframe = properties.loc[row["inhomogeneity_id"]] - kwargs = aquifer_data(dataframe) - kwargs["xy"] = polygon_coordinates(row) - kwargs["order"] = row["order"] - kwargs["ndeg"] = row["ndegrees"] - kwargs["layers"] = np.atleast_1d(row["layer"]) - kwarglist.append(kwargs) - return kwarglist - - -def headlinesink(spec: ElementSpecification) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "hls": row["head"], - "res": row["resistance"], - "wh": row["width"], - "order": row["order"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def linesinkditch(spec: ElementSpecification) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "Qls": row["discharge"], - "res": row["resistance"], - "wh": row["width"], - "order": row["order"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def leakylinedoublet(spec: ElementSpecification) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "res": row["resistance"], - "layers": row["layer"], - "order": row["order"], - "label": row["label"], - } - ) - return kwargslist - - -def implinedoublet(spec: ElementSpecification) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "layers": row["layer"], - "order": row["order"], - "label": row["label"], - } - ) - return kwargslist - - -def circareasink(spec: ElementSpecification) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - xc, yc = np.array(row["geometry"].centroid.coords)[0] - coords = np.array(row["geometry"].exterior.coords) - x, y = coords.T - # Use squared radii - radii2 = (x - xc) ** 2 + (y - yc) ** 2 - radius2 = radii2[0] - # Check whether geometry is close enough to a circle. - # Accept 1% deviation. - tolerance = 0.01 * radius2 - if not np.allclose(radii2, radius2, atol=tolerance): - raise ValueError("Circular Area Sink geometry is not circular") - - radius = np.sqrt(radius2) - kwargslist.append( - { - "xc": xc, - "yc": yc, - "R": radius, - "N": row["rate"], - } - ) - return kwargslist - - -def validate(spec: TimmlModelSpecification) -> None: - if spec.aquifer is None: - raise ValueError("Aquifer entry is missing") - # TODO: more checks - - -def head_observations(model: timml.Model, observations: Dict) -> gpd.GeoDataFrame: - geodataframes = {} - for name, kwargs_collection in observations.items(): - heads = [] - xx = [] - yy = [] - labels = [] - for kwargs in kwargs_collection: - x = kwargs["x"] - y = kwargs["y"] - heads.append(model.head(x=x, y=y)) - xx.append(x) - yy.append(y) - labels.append(kwargs["label"]) - - d = {"geometry": gpd.points_from_xy(xx, yy), "label": labels} - for i, layerhead in enumerate(np.vstack(heads).T): - d[f"head_layer{i}"] = layerhead - - geodataframes[name] = gpd.GeoDataFrame(d) - - return geodataframes - - -def headgrid(model: timml.Model, extent: Tuple[float], cellsize: float) -> xr.DataArray: - """ - Compute the headgrid of the TimML model, and store the results - in an xarray DataArray with the appropriate dimensions. - - Parameters - ---------- - model: timml.Model - Solved model to get heads from - extent: Tuple[float] - xmin, xmax, ymin, ymax - cellsize: float - Desired cell size of the output head grids - - Returns - ------- - head: xr.DataArray - DataArray with dimensions ``("layer", "y", "x")``. - """ - xmin, xmax, ymin, ymax = extent - x = np.arange(xmin, xmax, cellsize) + 0.5 * cellsize - # In geospatial rasters, y is DECREASING with row number - y = np.arange(ymax, ymin, -cellsize) - 0.5 * cellsize - nlayer = model.aq.find_aquifer_data(x[0], y[0]).naq - layer = [i for i in range(nlayer)] - head = np.empty((nlayer, y.size, x.size), dtype=np.float64) - for i in tqdm.tqdm(range(y.size)): - for j in range(x.size): - head[:, i, j] = model.head(x[j], y[i], layer) - - return xr.DataArray( - data=head, - name="head", - coords={"layer": layer, "y": y, "x": x}, - dims=("layer", "y", "x"), - ) - - -# Map the names of the elements to their constructors -MAPPING = { - "Constant": (constant, timml.Constant), - "Uniform Flow": (uflow, timml.Uflow), - "Circular Area Sink": (circareasink, timml.CircAreaSink), - "Well": (well, timml.Well), - "Head Well": (headwell, timml.HeadWell), - "Polygon Inhomogeneity": (polygoninhom, timml.PolygonInhomMaq), - "Polygon Area Sink": (polygonareasink, timml.PolygonInhomMaq), - "Polygon Semi-Confined Top": (polygontop, timml.PolygonInhomMaq), - "Head Line Sink": (headlinesink, timml.HeadLineSinkString), - "Line Sink Ditch": (linesinkditch, timml.LineSinkDitchString), - "Leaky Line Doublet": (leakylinedoublet, timml.LeakyLineDoubletString), - "Impermeable Line Doublet": (implinedoublet, timml.ImpLineDoubletString), - "Building Pit": (buildingpit, timml.BuildingPit), - "Observation": (observation, None), -} - - -def initialize_model(spec: TimmlModelSpecification) -> timml.Model: - """ - Initialize a TimML analytic model based on the data in a geopackage. - - Parameters - ---------- - spec: ModelSpecification - Named tuple with the layer name of the aquifer and a dictionary of - other element names to its construction function. - - Returns - ------- - timml.Model - elements: Dict of number name to TimML element - observations: Dict of list - - Examples - -------- - - >>> import gistim - >>> path = "my-model.gpkg" - >>> spec = gistim.model_specification(path) - >>> model = gistim.initialize_model(spec) - >>> model.solve() - - """ - validate(spec) - model = timml.ModelMaq(**aquifer_data(spec.aquifer)) - elements = {} - observations = defaultdict(list) - for name, element_spec in spec.elements.items(): - if (not element_spec.active) or (len(element_spec.dataframe.index) == 0): - continue - - elementtype = element_spec.elementtype - if elementtype not in MAPPING: - msg = ( - f'Invalid element specification "{elementtype}". ' - f'Available types are: {", ".join(MAPPING.keys())}.' - ) - raise KeyError(msg) - - f_to_kwargs, element = MAPPING[elementtype] - for i, kwargs in enumerate(f_to_kwargs(element_spec)): - if elementtype == "Observation": - observations[name].append(kwargs) - else: - kwargs["model"] = model - elements[f"{name}_{i}"] = element(**kwargs) - - return model, elements, observations - - -def convert_to_script(spec: TimmlModelSpecification) -> str: - """ - Convert model specification to an equivalent Python script. - """ - modelkwargs = dict_to_kwargs_code(aquifer_data(spec.aquifer)) - - observations = {} - strings = [ - "import numpy as np", - "import timml", - f"model = timml.ModelMaq({modelkwargs})", - ] - for name, element_spec in spec.elements.items(): - elementtype = element_spec.elementtype - - if elementtype not in MAPPING: - msg = ( - f'Invalid element specification "{elementtype}". ' - f'Available types are: {", ".join(MAPPING.keys())}.' - ) - raise KeyError(msg) - - f_to_kwargs, element = MAPPING[elementtype] - for i, kwargs in enumerate(f_to_kwargs(element_spec)): - if elementtype == "Observation": - kwargs.pop("label") - kwargs = dict_to_kwargs_code(kwargs) - observations[f"observation_{sanitized(name)}_{i}"] = kwargs - else: - kwargs["model"] = "model" - kwargs = dict_to_kwargs_code(kwargs) - strings.append( - f"{sanitized(name)}_{i} = timml.{element.__name__}({kwargs})" - ) - - strings.append("model.solve()") - - xg, yg = headgrid_code(spec.domain) - strings.append(f"head = model.headgrid(xg={xg}, yg={yg})") - - # Add all the individual observation points - for name, kwargs in observations.items(): - strings.append(f"{name} = model.head({kwargs})") - - return black.format_str("\n".join(strings), mode=black.FileMode()) diff --git a/gistim/ttim_elements.py b/gistim/ttim_elements.py deleted file mode 100644 index 1885752..0000000 --- a/gistim/ttim_elements.py +++ /dev/null @@ -1,442 +0,0 @@ -""" -This model converts the geodataframe as read from the geopackage into keyword -arguments for TTim. - -These keyword arguments are used to initialize a model, or used to generate a -Python script. - -""" -from collections import defaultdict -from typing import Any, Dict, List, Tuple - -import black -import geopandas as gpd -import numpy as np -import pandas as pd -import timml -import tqdm -import xarray as xr - -try: - import ttim - from ttim.model import TimModel -except ImportError: - ttim = None - TimModel = None - -from gistim.common import ( - ElementSpecification, - FloatArray, - TransientElementSpecification, - TtimModelSpecification, - aquifer_data, - dict_to_kwargs_code, - filter_nan, - headgrid_code, - linestring_coordinates, - point_coordinates, - sanitized, -) - - -# Dataframe to ttim element -# -------------------------- -def transient_dataframe(spec): - if spec.dataframe is not None: - return spec.dataframe.set_index("timeseries_id") - else: - return None - - -def transient_input(timeseries_df, row, field, time_start) -> List: - timeseries_id = row["timeseries_id"] - row_start = row.get("time_start") - row_end = row.get("time_end") - transient_value = row.get(f"{field}_transient") - if row_start is not None and row_end is not None and transient_value is not None: - steady_value = row[field] - return [(row_start, transient_value - steady_value), (row_end, 0.0)] - elif timeseries_df is None or timeseries_id not in timeseries_df.index: - return [(time_start, 0.0)] - else: - timeseries = timeseries_df.loc[[timeseries_id]] - timeseries[field] -= row[field] - timeseries = timeseries.sort_values(by="time_start") - return list( - timeseries[["time_start", field]].itertuples(index=False, name=None) - ) - - -def ttim_model( - dataframe: gpd.GeoDataFrame, - temporal_settings: pd.DataFrame, - timml_model: timml.Model, -) -> Dict[str, Any]: - """ - Parameters - ---------- - dataframe: geopandas.GeoDataFrame - - Returns - ------- - ttim.TimModel - """ - data = aquifer_data(dataframe, transient=True) - row = temporal_settings.iloc[0].to_dict() - data["tstart"] = row["time_start"] - data["tmin"] = row["time_min"] - data["tmax"] = row["time_max"] - data["M"] = row["stehfest_M"] - data["timmlmodel"] = timml_model - return data - - -def observation(spec: TransientElementSpecification, _): - df = transient_dataframe(spec) - if len(df) == 0: - return {} - - dataframe = spec.steady_spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - for i, (row, x, y) in enumerate(zip(dataframe.to_dict("records"), X, Y)): - row = filter_nan(row) - timeseries_id = row["timeseries_id"] - kwargslist.append( - { - "x": x, - "y": y, - "t": df.loc[timeseries_id, "time"].values, - "label": row["label"], - "observation_id": i, - } - ) - return kwargslist - - -def well( - spec: TransientElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - df = transient_dataframe(spec) - dataframe = spec.steady_spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - for (row, x, y) in zip(dataframe.to_dict("records"), X, Y): - row = filter_nan(row) - kwargslist.append( - { - "xw": x, - "yw": y, - "tsandQ": transient_input(df, row, "discharge", time_start), - "rw": row["radius"], - "res": row["resistance"], - "layers": row["layer"], - "label": row["label"], - "rc": row["caisson_radius"], - "wbstype": "slug" if row["slug"] else "pumping", - } - ) - return kwargslist - - -def headwell( - spec: TransientElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - df = transient_dataframe(spec) - dataframe = spec.steady_spec.dataframe - X, Y = point_coordinates(dataframe) - kwargslist = [] - for (row, x, y) in zip(dataframe.to_dict("records"), X, Y): - row = filter_nan(row) - kwargslist.append( - { - "xw": x, - "yw": y, - "tsandh": transient_input(df, row, "head", time_start), - "rw": row["radius"], - "res": row["resistance"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def headlinesink( - spec: TransientElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - df = transient_dataframe(spec) - kwargslist = [] - for row in spec.steady_spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "tsandh": transient_input(df, row, "head", time_start), - "res": row["resistance"], - "wh": row["width"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def linesinkditch( - spec: ElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - df = transient_dataframe(spec) - kwargslist = [] - for row in spec.steady_spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "tsandQ": transient_input(df, row, "discharge", time_start), - "res": row["resistance"], - "wh": row["width"], - "layers": row["layer"], - "label": row["label"], - } - ) - return kwargslist - - -def leakylinedoublet( - spec: ElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "res": row["resistance"], - "layers": row["layer"], - "order": row["order"], - "label": row["label"], - } - ) - return kwargslist - - -def implinedoublet( - spec: ElementSpecification, time_start: float -) -> List[Dict[str, Any]]: - kwargslist = [] - for row in spec.dataframe.to_dict("records"): - row = filter_nan(row) - kwargslist.append( - { - "xy": linestring_coordinates(row), - "res": "imp", - "layers": row["layer"], - "order": row["order"], - "label": row["label"], - } - ) - return kwargslist - - -def circareasink(spec: TransientElementSpecification, time_start: float) -> None: - """ - Parameters - ---------- - dataframe: geopandas.GeoDataFrame - - Returns - ------- - None - """ - df = transient_dataframe(spec) - kwargslist = [] - for row in spec.steady_spec.dataframe.to_dict("records"): - row = filter_nan(row) - x, y = np.array(row["geometry"].centroid.coords)[0] - coords = np.array(row["geometry"].exterior.coords) - x0, y0 = coords[0] - radius = np.sqrt((x0 - x) ** 2 + (y0 - y) ** 2) - kwargslist.append( - { - "xc": x, - "yc": y, - "R": radius, - "tsandN": transient_input(df, row, "rate", time_start), - } - ) - return kwargslist - - -def head_observations( - model: TimModel, - reference_date: pd.Timestamp, - observations: Dict, -) -> gpd.GeoDataFrame: - # We'll duplicate all values in time, except head which is unique per time. - refdate = pd.to_datetime(reference_date) - geodataframes = {} - for name, kwargs_collections in observations.items(): - xx = [] - yy = [] - labels = [] - heads = [] - starts = [] - ends = [] - observation_ids = [] - for kwargs in kwargs_collections: - x = kwargs["x"] - y = kwargs["y"] - t = kwargs["t"] - end = refdate + pd.to_timedelta(t, "D") - start = np.insert(end[:-1], 0, refdate) - - starts.append(start) - ends.append(end) - heads.append(model.head(x=x, y=y, t=t)) - xx.append(np.repeat(x, len(t))) - yy.append(np.repeat(y, len(t))) - labels.append(np.repeat(kwargs["label"], len(t))) - observation_ids.append(np.repeat(kwargs["observation_id"], len(t))) - - d = { - "geometry": gpd.points_from_xy(np.concatenate(xx), np.concatenate(yy)), - "datetime_start": np.concatenate(starts), - "datetime_end": np.concatenate(ends), - "label": np.concatenate(labels), - "observation_id": np.concatenate(observation_ids), - } - for i, layerhead in enumerate(np.hstack(heads)): - d[f"head_layer{i}"] = layerhead - - geodataframes[name] = gpd.GeoDataFrame(d) - - return geodataframes - - -def headgrid( - model: TimModel, - extent: Tuple[float], - cellsize: float, - times: FloatArray, - reference_date: pd.Timestamp, -) -> xr.DataArray: - xmin, xmax, ymin, ymax = extent - x = np.arange(xmin, xmax, cellsize) + 0.5 * cellsize - # In geospatial rasters, y is DECREASING with row number - y = np.arange(ymax, ymin, -cellsize) - 0.5 * cellsize - nlayer = model.aq.find_aquifer_data(x[0], y[0]).naq - layer = [i for i in range(nlayer)] - head = np.empty((nlayer, times.size, y.size, x.size), dtype=np.float64) - for i in tqdm.tqdm(range(y.size)): - for j in range(x.size): - head[:, :, i, j] = model.head(x[j], y[i], times, layer) - - time = pd.to_datetime(reference_date) + pd.to_timedelta(times, "D") - return xr.DataArray( - data=head, - name="head", - coords={"layer": layer, "time": time, "y": y, "x": x}, - dims=("layer", "time", "y", "x"), - ) - - -# Map the names of the elements to their constructors -if ttim is not None: - MAPPING = { - "Circular Area Sink": (circareasink, ttim.CircAreaSink), - "Well": (well, ttim.Well), - "Head Well": (headwell, ttim.HeadWell), - "Head Line Sink": (headlinesink, ttim.HeadLineSinkString), - "Line Sink Ditch": (linesinkditch, ttim.LineSinkDitchString), - "Leaky Line Doublet": (leakylinedoublet, ttim.LeakyLineDoubletString), - "Impermeable Line Doublet": (implinedoublet, ttim.LeakyLineDoubletString), - "Observation": (observation, None), - } - - -def initialize_model(spec: TtimModelSpecification, timml_model) -> TimModel: - """ - Initialize a Ttim analytic model based on the data in a geopackage. - - Parameters - ---------- - spec: ModelSpecification - Named tuple with the layer name of the aquifer and a dictionary of - other element names to its construction function. - - Returns - ------- - timml.Model - - Examples - -------- - - >>> import gistim - >>> path = "my-model.gpkg" - >>> spec = gistim.model_specification(path) - >>> model = gistim.initialize_model(spec) - >>> model.solve() - - """ - model = ttim.ModelMaq( - **ttim_model(spec.aquifer, spec.temporal_settings, timml_model) - ) - elements = {} - observations = defaultdict(list) - for name, element_spec in spec.elements.items(): - elementtype = element_spec.elementtype - if (not element_spec.active) or (elementtype not in MAPPING): - continue - - # print(f"adding {name} as {elementtype}") - f_to_kwargs, element = MAPPING[elementtype] - for i, kwargs in enumerate(f_to_kwargs(element_spec, model.tstart)): - if elementtype == "Observation": - observations[name].append(kwargs) - else: - kwargs["model"] = model - elements[f"{name}_{i}"] = element(**kwargs) - - return model, elements, observations - - -def convert_to_script(spec: TtimModelSpecification) -> str: - """ - Convert model specification to an equivalent Python script. - """ - modelkwargs = ttim_model(spec.aquifer, spec.temporal_settings, "model") - strkwargs = dict_to_kwargs_code(modelkwargs) - - observations = {} - strings = ["import ttim", f"ttim_model = ttim.ModelMaq({strkwargs})"] - for name, element_spec in spec.elements.items(): - elementtype = element_spec.elementtype - if elementtype not in MAPPING: - continue - - f_to_kwargs, element = MAPPING[elementtype] - for i, kwargs in enumerate(f_to_kwargs(element_spec, modelkwargs["tstart"])): - if elementtype == "Observation": - kwargs.pop("label") - kwargs.pop("observation_id") - kwargs = dict_to_kwargs_code(kwargs) - observations[f"ttim_observation_{sanitized(name)}_{i}"] = kwargs - else: - kwargs["model"] = "ttim_model" - kwargs = dict_to_kwargs_code(kwargs) - strings.append( - f"ttim_{sanitized(name)}_{i} = ttim.{element.__name__}({kwargs})" - ) - - strings.append("ttim_model.solve()") - - xg, yg = headgrid_code(spec.domain) - times = spec.output_times.tolist() - if len(times) > 0: - strings.append(f"ttim_head = ttim_model.headgrid(xg={xg}, yg={yg}, t={times})") - - # Add all the individual observation points - for name, kwargs in observations.items(): - strings.append(f"{name} = ttim_model.head({kwargs})") - - return black.format_str("\n".join(strings), mode=black.FileMode()) diff --git a/plugin/qgistim/core/dummy_ugrid.py b/plugin/qgistim/core/dummy_ugrid.py deleted file mode 100644 index e17c2c9..0000000 --- a/plugin/qgistim/core/dummy_ugrid.py +++ /dev/null @@ -1,41 +0,0 @@ -""" -This module contains a workaround for the following issue: -https://github.com/qgis/QGIS/issues/44474 - -QGIS does not release the file properly. A mesh layer with new results will not -be visualized properly. This can be worked around by forcing QGIS to load -another (dummy) file with UGRID mesh data. This "passes on" the file lock to -the newly loaded dummy file, and frees the older one. - -This requires always carrying around a dummy netCDF file, or generating it -with Python, then writing it. A netCDF is mostly easily generated using -xarray, at the cost of xarray and netCDF4 as dependencies. - -Instead, this module contains the binary representation of a netCDF4 file -that has been generated with the following lines of code: - -```python -import base64 -import numpy as np -import xarray as xr - -ds = xr.Dataset() -ds["dummy"] = xr.DataArray( - [[np.nan]], {"y": [0.], "x": [0.]}, ("y", "x") -) -ds.to_netcdf("dummy.nc") -with open("dummy.nc", "rb"): - binary_content = f.read() -b64_encoded = base64.b64encode(binary_content) -``` - -""" -import base64 -from pathlib import Path - -b64_encoded = "iUhERg0KGgoAAAAAAAgIAAQAEAAAAAAAAAAAAAAAAAD//////////7UkAAAAAAAA//////////8AAAAAAAAAAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAT0hEUgIM4AIiAAAAAAADBAAAAAAAAAD///////////////////////////////8KAgABAAAAABAQAAAAAMENAAAAAAAAQAAAAAAAAAAVHAAEAAAAAwEA////////////////////////////////DEoAAAAAAQAOAAgACABfTkNQcm9wZXJ0aWVzAAAAEwAAACIAAAABAAAAAAAAAHZlcnNpb249MixuZXRjZGY9NC43LjQsaGRmNT0xLjEwLjYGFAAAAAABBAIAAAAAAAAAAXnfAgAAAAAAAAAIAAAAAAAAAAAAAAAAnsdezk9IRFICDYgBARgAAAAAAQEBAAAAAAABAAAAAAAAAAEAAAAAAAAAAxQAAQAAESA/AAgAAAAAAEAANAsANP8DAAAQEAAAAAC3DwAAAAAAANwAAAAAAAAABAwAAQAACAAAAAAAAAAAAPh/CBIAAAAAAwFBBAAAAAAAAAgAAAAAAAAAFRwABAAAAAMFAP///////////////////////////////wwwAAAAAAEABgAIAAgAQ0xBU1MAAAATAAAAEAAAAAEAAAAAAAAARElNRU5TSU9OX1NDQUxFAAwiAAABAAEABQAIAAgATkFNRQAAAAATAAAAAgAAAAEAAAAAAAAAeAAMNAAAAgABAA4ADAAIAF9OZXRjZGY0RGltaWQAAAAQCAAABAAAAAAAIAAAAAAAAQAAAAAAAAAAAAAADFAAAAMAAQALABQAGABfRmlsbFZhbHVlAAAAAAAAESA/AAgAAAAAAEAANAsANP8DAAAAAAAAAQEBAAAAAAABAAAAAAAAAAEAAAAAAAAAAAAAAAAA+H8KLddoT0hEUgINMgEBGAAAAAABAQEAAAAAAAEAAAAAAAAAAQAAAAAAAAADFAABAAARID8ACAAAAAAAQAA0CwA0/wMAABAQAAAAAEEMAAAAAAAAdAAAAAAAAAAEDAABAAAIAAAAAAAAAAAA+H8IEgAAAAADAUkEAAAAAAAACAAAAAAAAAAVHAAEAAAAAwUA////////////////////////////////DDAAAAAAAQAGAAgACABDTEFTUwAAABMAAAAQAAAAAQAAAAAAAABESU1FTlNJT05fU0NBTEUADCIAAAEAAQAFAAgACABOQU1FAAAAABMAAAACAAAAAQAAAAAAAAB5AAw0AAACAAEADgAMAAgAX05ldGNkZjREaW1pZAAAABAIAAAEAAAAAAAgAAAAAAABAAAAAAAAAAEAAACugT5IZAAAABAIAAAEAAAAAAAgAAAAAAABAAAAAAAAAAEAAAAE/QA+AAAAAAAAAAAAAAAAAAAAAAAAAAAAAPh/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABPQ0hLEBAAAAAA2w4AAAAAAADcAAAAAAAAAAxQAAADAAEACwAUABgAX0ZpbGxWYWx1ZQAAAAAAABEgPwAIAAAAAABAADQLADT/AwAAAAAAAAEBAQAAAAAAAQAAAAAAAAABAAAAAAAAAAAAAAAAAPh/DFP14U9IRFICDQABASgAAAAAAQIBAAAAAAABAAAAAAAAAAEAAAAAAAAAAQAAAAAAAAABAAAAAAAAAAMUAAEAABEgPwAIAAAAAABAADQLADT/AwAAEBAAAAAAAQ4AAAAAAADaAAAAAAAAAAQMAAEAAAgAAAAAAAAAAAD4fwgSAAAAAAMBUQQAAAAAAAAIAAAAAAAAABUcAAQAAAADAwD///////////////////////////////8MUAAAAAABABQADAAYAF9OZXRjZGY0Q29vcmRpbmF0ZXMAAAAAABAIAAAEAAAAAAAgAAAAAAABAQEAAAAAAAIAAAAAAAAAAgAAAAAAAAABAAAAAAAAACnzvHpPQ0hLBhQAAAAAAQQAAAAAAAAAAAF4SwEAAAAAAAAGGAAAAAABBAMAAAAAAAAABWR1bW15tQwAAAAAAACBIi5kT0NISwUQAAEAAAICAgEIAAAAAAAAAAAA+H8MUAAAAQABAAsAFAAYAF9GaWxsVmFsdWUAAAAAAAARID8ACAAAAAAAQAA0CwA0/wMAAAAAAAABAQEAAAAAAAEAAAAAAAAAAQAAAAAAAAAAAAAAAAD4fwxgAAACAAEADwAQABgARElNRU5TSU9OX0xJU1QAABkAAAAQAAAAFwAAAAgAAAABAQEAAAAAAAIAAAAAAAAAAgAAAAAAAAABAAAAtRQAAAAAAAACAAAAAQAAALUUAAAAAAAAAwAAAFHJK5JPQ0hLBRAAAQAAAgICAQgAAAAAAAAAAAD4fwy4AAAEAAEADwB0ABgAUkVGRVJFTkNFX0xJU1QAABYCAAAQAAAAZGF0YXNldAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABcAAAAIAAAAZGltZW5zaW9uAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAEAgAAAQAAAAAACAAAAAAAAEBAQAAAAAAAQAAAAAAAAABAAAAAAAAALUMAAAAAAAAAAAAAAAAAACsBlGrT0NISwUQAAEAAAICAgEIAAAAAAAAAAAA+H8MuAAABAABAA8AdAAYAFJFRkVSRU5DRV9MSVNUAAAWAgAAEAAAAGRhdGFzZXQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAXAAAACAAAAGRpbWVuc2lvbgAAAAAAAAAIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABAIAAAEAAAAAAAgAAAAAAABAQEAAAAAAAEAAAAAAAAAAQAAAAAAAAC1DAAAAAAAAAEAAAAAAAAAsUXkDQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAR0NPTAEAAAAAEAAAAAAAAAEAAAAAAAAACAAAAAAAAADfAgAAAAAAAAIAAAAAAAAACAAAAAAAAADfAgAAAAAAAAMAAAAAAAAACAAAAAAAAABLAQAAAAAAAAAAAAAAAAAAqA8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA==" - - -def write_dummy_ugrid(path: Path) -> None: - with open(path, "wb") as f: - f.write(base64.b64decode(b64_encoded)) diff --git a/plugin/qgistim/core/elements/__init__.py b/plugin/qgistim/core/elements/__init__.py new file mode 100644 index 0000000..4985967 --- /dev/null +++ b/plugin/qgistim/core/elements/__init__.py @@ -0,0 +1,100 @@ +import re +from collections import defaultdict +from functools import partial +from typing import List, Tuple + +from qgistim.core import geopackage +from qgistim.core.elements.aquifer import Aquifer +from qgistim.core.elements.building_pit import BuildingPit +from qgistim.core.elements.circular_area_sink import CircularAreaSink +from qgistim.core.elements.constant import Constant +from qgistim.core.elements.domain import Domain +from qgistim.core.elements.element import Element +from qgistim.core.elements.head_line_sink import HeadLineSink +from qgistim.core.elements.headwell import HeadWell, RemoteHeadWell +from qgistim.core.elements.impermeable_line_doublet import ImpermeableLineDoublet +from qgistim.core.elements.leaky_building_pit import LeakyBuildingPit +from qgistim.core.elements.leaky_line_doublet import LeakyLineDoublet +from qgistim.core.elements.line_sink_ditch import LineSinkDitch +from qgistim.core.elements.observation import HeadObservation +from qgistim.core.elements.polygon_area_sink import PolygonAreaSink +from qgistim.core.elements.polygon_inhomogeneity import PolygonInhomogeneity +from qgistim.core.elements.polygon_semi_confined_top import PolygonSemiConfinedTop +from qgistim.core.elements.uniform_flow import UniformFlow +from qgistim.core.elements.well import Well + +ELEMENTS = { + element.element_type: element + for element in ( + Aquifer, + Domain, + Constant, + UniformFlow, + Well, + HeadWell, + RemoteHeadWell, + HeadLineSink, + LineSinkDitch, + CircularAreaSink, + ImpermeableLineDoublet, + LeakyLineDoublet, + PolygonAreaSink, + PolygonSemiConfinedTop, + PolygonInhomogeneity, + BuildingPit, + LeakyBuildingPit, + HeadObservation, + ) +} + + +def parse_name(layername: str) -> Tuple[str, str, str]: + """ + Based on the layer name find out: + + * whether it's a timml or ttim element; + * which element type it is; + * what the user provided name is. + + For example: + parse_name("timml Headwell:drainage") -> ("timml", "Head Well", "drainage") + """ + prefix, name = layername.split(":") + element_type = re.split("timml |ttim ", prefix)[1] + mapping = { + "Computation Times": "Domain", + "Temporal Settings": "Aquifer", + "Polygon Inhomogeneity Properties": "Polygon Inhomogeneity", + "Building Pit Properties": "Building Pit", + "Leaky Building Pit Properties": "Leaky Building Pit", + } + element_type = mapping.get(element_type, element_type) + if "timml" in prefix: + if "Properties" in prefix: + tim_type = "timml_assoc" + else: + tim_type = "timml" + elif "ttim" in prefix: + tim_type = "ttim" + else: + raise ValueError("Neither timml nor ttim in layername") + return tim_type, element_type, name + + +def load_elements_from_geopackage(path: str) -> List[Element]: + # List the names in the geopackage + gpkg_names = geopackage.layers(path) + + # Group them on the basis of name + dd = defaultdict + grouped_names = dd(partial(dd, partial(dd, list))) + for layername in gpkg_names: + tim_type, element_type, name = parse_name(layername) + grouped_names[element_type][name][tim_type] = layername + + elements = [] + for element_type, group in grouped_names.items(): + for name in group: + elements.append(ELEMENTS[element_type](path, name)) + + return elements diff --git a/plugin/qgistim/core/elements/aquifer.py b/plugin/qgistim/core/elements/aquifer.py new file mode 100644 index 0000000..35fd052 --- /dev/null +++ b/plugin/qgistim/core/elements/aquifer.py @@ -0,0 +1,129 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField +from qgistim.core import geopackage +from qgistim.core.elements.element import ElementExtraction, TransientElement +from qgistim.core.elements.schemata import SingleRowSchema, TableSchema +from qgistim.core.schemata import ( + AllGreaterEqual, + AllRequired, + OffsetAllRequired, + OptionalFirstOnly, + Positive, + Range, + Required, + SemiConfined, + StrictlyDecreasing, +) + + +class AquiferSchema(TableSchema): + timml_schemata = { + "layer": AllRequired(Range()), + "aquifer_top": AllRequired(StrictlyDecreasing()), + "aquifer_bottom": AllRequired(StrictlyDecreasing()), + "aquitard_c": OffsetAllRequired(Positive()), + "aquifer_k": AllRequired(Positive()), + "semiconf_top": OptionalFirstOnly(), + "semiconf_head": OptionalFirstOnly(), + } + timml_consistency_schemata = ( + SemiConfined(), + AllGreaterEqual("aquifer_top", "aquifer_bottom"), + ) + ttim_schemata = { + "aquitard_s": OffsetAllRequired(Positive()), + "aquifer_s": AllRequired(Positive()), + } + + +class TemporalSettingsSchema(SingleRowSchema): + ttim_schemata = { + "time_min": Required(Positive()), + "laplace_inversion_M": Required(Positive()), + "reference_date": Required(), + } + + +class Aquifer(TransientElement): + element_type = "Aquifer" + geometry_type = "No Geometry" + timml_attributes = [ + QgsField("layer", QVariant.Int), + QgsField("aquifer_top", QVariant.Double), + QgsField("aquifer_bottom", QVariant.Double), + QgsField("aquitard_c", QVariant.Double), + QgsField("aquifer_k", QVariant.Double), + QgsField("semiconf_top", QVariant.Double), + QgsField("semiconf_head", QVariant.Double), + QgsField("aquitard_s", QVariant.Double), + QgsField("aquifer_s", QVariant.Double), + QgsField("aquitard_npor", QVariant.Double), + QgsField("aquifer_npor", QVariant.Double), + ] + ttim_attributes = ( + QgsField("time_min", QVariant.Double), + QgsField("laplace_inversion_M", QVariant.Int), + QgsField("reference_date", QVariant.DateTime), + ) + ttim_defaults = { + "time_min": QgsDefaultValue("0.01"), + "laplace_inversion_M": QgsDefaultValue("10"), + } + transient_columns = ( + "aquitard_s", + "aquifer_s", + "aquitard_npor", + "aquifer_npor", + ) + schema = AquiferSchema() + assoc_schema = TemporalSettingsSchema() + + def __init__(self, path: str, name: str): + self._initialize_default(path, name) + self.timml_name = f"timml {self.element_type}:Aquifer" + self.ttim_name = "ttim Temporal Settings:Aquifer" + + def write(self): + self.timml_layer = geopackage.write_layer( + self.path, self.timml_layer, self.timml_name, newfile=True + ) + self.ttim_layer = geopackage.write_layer( + self.path, self.ttim_layer, self.ttim_name + ) + self.set_defaults() + + def remove_from_geopackage(self): + """This element may not be removed.""" + return + + def to_timml(self) -> ElementExtraction: + missing = self.check_timml_columns() + if missing: + return ElementExtraction(errors=missing) + + data = self.table_to_dict(layer=self.timml_layer) + errors = self.schema.validate_timml(name=self.timml_layer.name(), data=data) + return ElementExtraction(errors=errors, data=data) + + def to_ttim(self) -> ElementExtraction: + missing = self.check_ttim_columns() + if missing: + return ElementExtraction(errors=missing) + + data = self.table_to_dict(layer=self.timml_layer) + time_data = self.table_to_records(layer=self.ttim_layer) + errors = { + **self.schema.validate_ttim(name=self.timml_layer.name(), data=data), + **self.assoc_schema.validate_ttim( + name=self.ttim_layer.name(), data=time_data + ), + } + if errors: + return ElementExtraction(errors=errors) + return ElementExtraction(data={**data, **time_data[0]}) + + def extract_data(self, transient: bool) -> ElementExtraction: + if transient: + return self.to_ttim() + else: + return self.to_timml() diff --git a/plugin/qgistim/core/elements/building_pit.py b/plugin/qgistim/core/elements/building_pit.py new file mode 100644 index 0000000..4b399e4 --- /dev/null +++ b/plugin/qgistim/core/elements/building_pit.py @@ -0,0 +1,100 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import RED, TRANSPARENT_RED +from qgistim.core.elements.element import AssociatedElement +from qgistim.core.elements.schemata import RowWiseSchema, TableSchema +from qgistim.core.schemata import ( + AllGreaterEqual, + AllRequired, + AtleastOneTrue, + Membership, + OffsetAllRequired, + OptionalFirstOnly, + Positive, + Range, + Required, + SemiConfined, + StrictlyDecreasing, +) + + +class BuildingPitSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "inhomogeneity_id": Required(Membership("properties inhomogeneity_id")), + "order": Required(Positive()), + "ndegrees": Required(Positive()), + } + + +class AssociatedBuildingPitSchema(TableSchema): + timml_schemata = { + "inhomogeneity_id": AllRequired(), + "layer": AllRequired(Range()), + "aquifer_top": AllRequired(StrictlyDecreasing()), + "aquifer_bottom": AllRequired(StrictlyDecreasing()), + "aquitard_c": OffsetAllRequired(Positive()), + "aquifer_k": AllRequired(Positive()), + "semiconf_top": OptionalFirstOnly(), + "semiconf_head": OptionalFirstOnly(), + "wall_in_layer": AllRequired(AtleastOneTrue()), + } + timml_consistency_schemata = ( + SemiConfined(), + AllGreaterEqual("aquifer_top", "aquifer_bottom"), + ) + + +class BuildingPit(AssociatedElement): + element_type = "Building Pit" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("order", QVariant.Int), + QgsField("ndegrees", QVariant.Int), + ) + assoc_attributes = [ + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("aquifer_top", QVariant.Double), + QgsField("aquifer_bottom", QVariant.Double), + QgsField("aquitard_c", QVariant.Double), + QgsField("aquifer_k", QVariant.Double), + QgsField("semiconf_top", QVariant.Double), + QgsField("semiconf_head", QVariant.Double), + QgsField("wall_in_layer", QVariant.Bool), + QgsField("aquitard_npor", QVariant.Double), + QgsField("aquifer_npor", QVariant.Double), + ] + timml_defaults = { + "order": QgsDefaultValue("4"), + "ndegrees": QgsDefaultValue("6"), + "inhomogeneity_id": QgsDefaultValue("1"), + } + assoc_defaults = { + "inhomogeneity_id": QgsDefaultValue("1"), + "wall_in_layer": QgsDefaultValue("false"), + } + schema = BuildingPitSchema() + assoc_schema = AssociatedBuildingPitSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.polygon_renderer( + color=TRANSPARENT_RED, color_border=RED, width_border="0.75" + ) + + def process_timml_row(self, row: Dict[str, Any], grouped: Dict[int, Any]): + inhom_id = row["inhomogeneity_id"] + raw_data = grouped[inhom_id] + layers = [i for i, active in enumerate(raw_data["wall_in_layer"]) if active] + aquifer_data = self.aquifer_data(raw_data, transient=False) + return { + "xy": self.polygon_xy(row), + "order": row["order"], + "ndeg": row["ndegrees"], + "layers": layers, + **aquifer_data, + } diff --git a/plugin/qgistim/core/elements/circular_area_sink.py b/plugin/qgistim/core/elements/circular_area_sink.py new file mode 100644 index 0000000..205cb2c --- /dev/null +++ b/plugin/qgistim/core/elements/circular_area_sink.py @@ -0,0 +1,100 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from qgis.core import QgsField +from qgistim.core.elements.colors import GREEN, TRANSPARENT_GREEN +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import ( + AllOrNone, + AllRequired, + CircularGeometry, + Membership, + NotBoth, + Optional, + Positive, + Required, +) + + +class CircularAreaSinkSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(CircularGeometry()), + "rate": Required(), + "layer": Required(Membership("aquifer layers")), + } + ttim_schemata = { + "time_start": Optional(Positive()), + "time_end": Optional(Positive()), + "timeseries_id": Optional(Membership("ttim timeseries IDs")), + } + ttim_consistency_schemata = ( + AllOrNone("time_start", "time_end", "rate_transient"), + NotBoth("time_start", "timeseries_id"), + ) + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time_start": AllRequired(Positive()), + "rate": AllRequired(), + } + + +class CircularAreaSink(TransientElement): + element_type = "Circular Area Sink" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("rate", QVariant.Double), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + QgsField("time_start", QVariant.Double), + QgsField("time_end", QVariant.Double), + QgsField("rate_transient", QVariant.Double), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time_start", QVariant.Double), + QgsField("rate", QVariant.Double), + ) + transient_columns = ( + "time_start", + "time_end", + "rate_transient", + "timeseries_id", + ) + schema = CircularAreaSinkSchema() + + @classmethod + def renderer(cls): + return cls.polygon_renderer( + color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75" + ) + + def _centroid_and_radius(self, row): + # Take the first vertex. + x, y = self.point_xy(row) + # Compare with the centroid to derive radius. + xc, yc = row["centroid"] + radius = ((x - xc) ** 2 + (y - yc) ** 2) ** 0.5 + return xc, yc, radius + + def process_timml_row(self, row, other=None) -> Dict[str, Any]: + xc, yc, radius = self._centroid_and_radius(row) + return { + "xc": xc, + "yc": yc, + "R": radius, + "N": row["rate"], + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + xc, yc, radius = self._centroid_and_radius(row) + tsandN, times = self.transient_input(row, grouped, "rate") + return { + "xc": xc, + "yc": yc, + "R": radius, + "tsandN": tsandN, + "label": row["label"], + }, times diff --git a/plugin/qgistim/core/elements/colors.py b/plugin/qgistim/core/elements/colors.py new file mode 100644 index 0000000..302e338 --- /dev/null +++ b/plugin/qgistim/core/elements/colors.py @@ -0,0 +1,13 @@ +RED = "215,48,39,255" +GREEN = "51,160,44,255" +BLUE = "31,120,180,255" +GREY = "135,135,135,255" +BLACK = "0,0,0,255" +LIGHT_BLUE = "166,206,227,255" + +OPACITY = 51 # 20% +TRANSPARENT = "255,255,255,0" +TRANSPARENT_RED = f"215,48,39,{OPACITY}" +TRANSPARENT_GREEN = f"51,160,44,{OPACITY}" +TRANSPARENT_BLUE = f"31,120,180,{OPACITY}" +TRANSPARENT_GREY = f"135,135,135,{OPACITY}" diff --git a/plugin/qgistim/core/elements/constant.py b/plugin/qgistim/core/elements/constant.py new file mode 100644 index 0000000..ddd3465 --- /dev/null +++ b/plugin/qgistim/core/elements/constant.py @@ -0,0 +1,39 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import RED +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import SingleRowSchema +from qgistim.core.schemata import Membership, Required + + +class ConstantSchema(SingleRowSchema): + timml_schemata = { + "geometry": Required(), + "head": Required(), + "layer": Required(Membership("aquifer layers")), + } + + +class Constant(Element): + element_type = "Constant" + geometry_type = "Point" + timml_attributes = ( + QgsField("head", QVariant.Double), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + ) + schema = ConstantSchema() + + @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) + return { + "xr": x, + "yr": y, + "hr": row["head"], + "layer": row["layer"], + "label": row["label"], + } diff --git a/plugin/qgistim/core/elements/domain.py b/plugin/qgistim/core/elements/domain.py new file mode 100644 index 0000000..9149c60 --- /dev/null +++ b/plugin/qgistim/core/elements/domain.py @@ -0,0 +1,102 @@ +from typing import Any, Tuple + +from PyQt5.QtCore import QVariant +from qgis.core import ( + QgsFeature, + QgsField, + QgsGeometry, + QgsPointXY, + QgsSingleSymbolRenderer, +) +from qgistim.core.elements.colors import BLACK +from qgistim.core.elements.element import ElementExtraction, TransientElement +from qgistim.core.elements.schemata import SingleRowSchema +from qgistim.core.schemata import AllRequired, Positive, Required, StrictlyIncreasing + + +class DomainSchema(SingleRowSchema): + timml_schemata = {"geometry": Required()} + timeseries_schemata = { + "time": AllRequired(Positive(), StrictlyIncreasing()), + } + + +class Domain(TransientElement): + element_type = "Domain" + geometry_type = "Polygon" + ttim_attributes = (QgsField("time", QVariant.Double),) + schema = DomainSchema() + + def __init__(self, path: str, name: str): + self._initialize_default(path, name) + self.timml_name = f"timml {self.element_type}:Domain" + self.ttim_name = "ttim Computation Times:Domain" + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + """ + Results in transparent fill, with a medium thick black border line. + """ + return cls.polygon_renderer( + color="255,0,0,0", color_border=BLACK, width_border="0.75" + ) + + def remove_from_geopackage(self): + pass + + def update_extent(self, iface: Any) -> Tuple[float, float]: + provider = self.timml_layer.dataProvider() + provider.truncate() # removes all features + canvas = iface.mapCanvas() + extent = canvas.extent() + xmin = extent.xMinimum() + ymin = extent.yMinimum() + xmax = extent.xMaximum() + ymax = extent.yMaximum() + points = [ + QgsPointXY(xmin, ymax), + QgsPointXY(xmax, ymax), + QgsPointXY(xmax, ymin), + QgsPointXY(xmin, ymin), + ] + feature = QgsFeature() + feature.setGeometry(QgsGeometry.fromPolygonXY([points])) + provider.addFeatures([feature]) + canvas.refresh() + return ymax, ymin + + def to_timml(self, other) -> ElementExtraction: + data = self.table_to_records(layer=self.timml_layer) + errors = self.schema.validate_timml( + name=self.timml_layer.name(), data=data, other=other + ) + if errors: + return ElementExtraction(errors=errors) + else: + x = [point[0] for point in data[0]["geometry"]] + y = [point[1] for point in data[0]["geometry"]] + return ElementExtraction( + data={ + "xmin": min(x), + "xmax": max(x), + "ymin": min(y), + "ymax": max(y), + } + ) + + def to_ttim(self, other) -> ElementExtraction: + timml_extraction = self.to_timml(other) + data = timml_extraction.data + + timeseries = self.table_to_dict(layer=self.ttim_layer) + errors = self.schema.validate_timeseries( + name=self.ttim_layer.name(), data=timeseries + ) + if errors: + return ElementExtraction(errors=errors) + if timeseries["time"]: + data["time"] = timeseries["time"] + times = set(timeseries["time"]) + else: + times = set() + return ElementExtraction(data=data, times=times) diff --git a/plugin/qgistim/core/elements/element.py b/plugin/qgistim/core/elements/element.py new file mode 100644 index 0000000..63227ff --- /dev/null +++ b/plugin/qgistim/core/elements/element.py @@ -0,0 +1,552 @@ +""" +This module contains the classes to represent the TimML and TTim input layers. + +The classes specify: + +* The (unabbreviated) name +* The type of geometry (No geometry, point, linestring, polygon) +* The required attributes of the attribute table + +They contain logic for setting up: + +* Simple input, requiring a single table, e.g. Uniform Flow or Constant +* Transient input, requiring two tables, one with geometry and steady-state + properties, and one containing optional time series input e.g. Well, Head Line + Sink. +* Associated input, requiring two tables, one with geometry and one with input + for layers, e.g. Polygon Inhomogeneity or Building Pit. + +Each element is (optionally) represented in multiple places: + +* It always lives in a GeoPackage. +* While a geopackage is active within plugin, it is always represented in a + Dataset Tree: the Dataset Tree provides a direct look at the state of the + GeoPackage. In this tree, steady and transient input are on the same row. + Associated input is, to potentially enable transient associated data later + on (like a building pit with changing head top boundary). +* It can be added to the Layers Panel in QGIS. This enables a user to visualize + and edit its data. + +Some elements require specific rendering in QGIS (e.g. no fill polygons), which +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 +associated with each other is by comparing names. Names must therefore be +unique within a group of the same type of elements. + +Rendering: + +* Fixed discharge (area sink, well, ditch): green +* Fixed head (head well, semi-confined top, head line): blue +* Inhomogeneity: grey +* Constant: star +* Observations: triangle +* Line Doublets: red +* Polygons: Line and Fill same color, Fill color 15% opacity. + +Nota bene: the order of the columns is important for hiding and showing the +transient columns. QGIS has a bug which causes it to show the wrong column if +hidden columns appear before shown ones. This only affects the attribute table +when it has no features. +""" + +import abc +from collections import defaultdict +from copy import deepcopy +from typing import Any, Dict, List, NamedTuple, Optional, Set, Tuple, Union + +from PyQt5.QtWidgets import ( + QDialog, + QHBoxLayout, + QLabel, + QLineEdit, + QPushButton, + QVBoxLayout, +) +from qgis.core import ( + QgsFillSymbol, + QgsLineSymbol, + QgsMarkerSymbol, + QgsSingleSymbolRenderer, + QgsVectorLayer, +) +from qgistim.core import geopackage +from qgistim.core.extractor import ExtractorMixin + + +class ElementExtraction(NamedTuple): + errors: Optional[Dict[str, Any]] = None + data: Optional[Union[Dict[str, Any], List[Dict[str, Any]]]] = None + times: Optional[Set[float]] = None + + +class NameDialog(QDialog): + def __init__(self, parent=None): + super().__init__(parent) + self.name_line_edit = QLineEdit() + self.ok_button = QPushButton("OK") + self.cancel_button = QPushButton("Cancel") + self.ok_button.clicked.connect(self.accept) + self.cancel_button.clicked.connect(self.reject) + first_row = QHBoxLayout() + first_row.addWidget(QLabel("Layer name")) + first_row.addWidget(self.name_line_edit) + second_row = QHBoxLayout() + second_row.addStretch() + second_row.addWidget(self.ok_button) + second_row.addWidget(self.cancel_button) + layout = QVBoxLayout() + layout.addLayout(first_row) + layout.addLayout(second_row) + self.setLayout(layout) + + +class Element(ExtractorMixin, abc.ABC): + """ + Abstract base class for "ordinary" timml elements. + """ + + element_type = None + geometry_type = None + timml_attributes = () + ttim_attributes = () + assoc_attributes = () + transient_columns = () + timml_defaults = {} + ttim_defaults = {} + assoc_defaults = {} + + def _initialize_default(self, path, name): + self.name = name + self.path = path + self.timml_name = None + self.ttim_name = None + self.assoc_name = None + self.timml_layer = None + self.ttim_layer = None + self.assoc_layer = None + self.item = None + + def __init__(self, path: str, name: str): + self._initialize_default(path, name) + self.timml_name = f"timml {self.element_type}:{name}" + + @classmethod + def dialog(cls, path: str, crs: Any, iface: Any, names: List[str]): + dialog = NameDialog() + dialog.show() + ok = dialog.exec_() + if not ok: + return + + name = dialog.name_line_edit.text() + if name in names: + raise ValueError(f"Name already exists in geopackage: {name}") + + instance = cls(path, name) + instance.create_layers(crs) + return instance + + def create_layer( + self, crs: Any, geometry_type: str, name: str, attributes: List + ) -> QgsVectorLayer: + layer = QgsVectorLayer(geometry_type, name, "memory") + provider = layer.dataProvider() + provider.addAttributes(attributes) + layer.updateFields() + layer.setCrs(crs) + return layer + + def create_timml_layer(self, crs: Any): + self.timml_layer = self.create_layer( + crs=crs, + geometry_type=self.geometry_type, + name=self.timml_name, + attributes=self.timml_attributes, + ) + + def create_ttim_layer(self, crs: Any): + pass + + def create_assoc_layer(self, crs: Any): + pass + + def create_layers(self, crs: Any): + self.create_timml_layer(crs) + self.create_ttim_layer(crs) + self.create_assoc_layer(crs) + + def set_defaults(self): + for layer, defaults in zip( + (self.timml_layer, self.ttim_layer, self.assoc_layer), + (self.timml_defaults, self.ttim_defaults, self.assoc_defaults), + ): + if layer is None: + continue + fields = layer.fields() + for name, definition in defaults.items(): + index = fields.indexFromName(name) + layer.setDefaultValueDefinition(index, definition) + return + + @staticmethod + def marker_renderer(**kwargs): + symbol = QgsMarkerSymbol.createSimple(kwargs) + return QgsSingleSymbolRenderer(symbol) + + @staticmethod + def line_renderer(**kwargs): + symbol = QgsLineSymbol.createSimple(kwargs) + return QgsSingleSymbolRenderer(symbol) + + @staticmethod + def polygon_renderer(**kwargs): + symbol = QgsFillSymbol.createSimple(kwargs) + return QgsSingleSymbolRenderer(symbol) + + @classmethod + def renderer(cls): + return None + + def timml_layer_from_geopackage(self) -> QgsVectorLayer: + self.timml_layer = QgsVectorLayer( + f"{self.path}|layername={self.timml_name}", self.timml_name + ) + + def ttim_layer_from_geopackage(self): + return + + def assoc_layer_from_geopackage(self): + return + + def load_layers_from_geopackage(self) -> None: + self.timml_layer_from_geopackage() + self.ttim_layer_from_geopackage() + self.assoc_layer_from_geopackage() + self.set_defaults() + return + + def write(self): + self.timml_layer = geopackage.write_layer( + self.path, self.timml_layer, self.timml_name + ) + self.set_defaults() + + def remove_from_geopackage(self): + geopackage.remove_layer(self.path, self.timml_name) + + def on_transient_changed(self, _): + return + + @staticmethod + def groupby(data: Dict[str, Any], key: str): + data = deepcopy(data) + grouper = data.pop(key) + grouped = {k: defaultdict(list) for k in grouper} + for key, values in data.items(): + for groupkey, value in zip(grouper, values): + grouped[groupkey][key].append(value) + return grouped + + @staticmethod + def _check_table_columns(attributes, layer) -> Dict[str, List]: + """ + Check if any columns are missing from the table. + + In that case, abort and present an error message. + """ + fields = set(field.name() for field in layer.fields()) + missing = set(attr.name() for attr in attributes) - fields + if missing: + columns = ",".join(missing) + msg = ( + f"Table is missing columns: {columns}. " + "Remove and recreate the layer." + ) + return {"Table:": [msg]} + return {} + + def check_timml_columns(self): + return self._check_table_columns( + attributes=self.timml_attributes, layer=self.timml_layer + ) + + def to_timml(self, other=None) -> ElementExtraction: + missing = self.check_timml_columns() + if missing: + return ElementExtraction(errors=missing) + + data = self.table_to_records(layer=self.timml_layer) + errors = self.schema.validate_timml( + name=self.timml_layer.name(), data=data, other=other + ) + + if errors: + return ElementExtraction(errors=errors) + else: + elements = [self.process_timml_row(row=row, other=other) for row in data] + return ElementExtraction(data=elements) + + def to_ttim(self, other=None) -> ElementExtraction: + return self.to_timml(other) + + def extract_data(self, transient: bool, other=None) -> ElementExtraction: + if transient: + return self.to_ttim(other) + else: + return self.to_timml(other) + + @staticmethod + def aquifer_data(data, transient: bool): + """Used by the Aquifer element and the Inhomogeneities.""" + + def interleave(a, b): + return [value for pair in zip(a, b) for value in pair] + + data = deepcopy(data) # Avoid side-effects + hstar = data["semiconf_head"][0] + semiconfined = ( + data["aquitard_c"][0] is not None + and data["semiconf_top"][0] is not None + and hstar is not None + ) + + kaq = data["aquifer_k"] + c = data["aquitard_c"] + z = [data["semiconf_top"][0]] + interleave( + data["aquifer_top"], data["aquifer_bottom"] + ) + porosity = interleave(data["aquitard_npor"], data["aquifer_npor"]) + s_aquifer = data["aquifer_s"] + s_aquitard = data["aquitard_s"] + + if semiconfined: + topboundary = "semi" + else: + topboundary = "conf" + z.pop(0) + c.pop(0) + porosity.pop(0) + if transient: + s_aquitard.pop(0) + + aquifer = { + "kaq": kaq, + "z": z, + "c": c, + "topboundary": topboundary, + } + + if transient: + aquifer["Saq"] = s_aquifer + aquifer["Sll"] = s_aquitard + aquifer["phreatictop"] = True + aquifer["tmin"] = data["time_min"] + aquifer["tstart"] = 0.0 + aquifer["M"] = data["laplace_inversion_M"] + else: + aquifer["npor"] = porosity + aquifer["hstar"] = hstar + + if "resistance" in data: + aquifer["res"] = data["resistance"][0] + + return aquifer + + +class TransientElement(Element, abc.ABC): + """ + Abstract base class for transient (ttim) elements. + """ + + def __init__(self, path: str, name: str): + self._initialize_default(path, name) + self.timml_name = f"timml {self.element_type}:{name}" + self.ttim_name = f"ttim {self.element_type}:{name}" + + def create_ttim_layer(self, crs: Any): + self.ttim_layer = self.create_layer( + crs=crs, + geometry_type="No Geometry", + name=self.ttim_name, + attributes=self.ttim_attributes, + ) + + def ttim_layer_from_geopackage(self): + self.ttim_layer = QgsVectorLayer( + f"{self.path}|layername={self.ttim_name}", + self.ttim_name, + ) + + def write(self): + self.timml_layer = geopackage.write_layer( + self.path, self.timml_layer, self.timml_name + ) + self.ttim_layer = geopackage.write_layer( + self.path, self.ttim_layer, self.ttim_name + ) + self.set_defaults() + + def remove_from_geopackage(self): + geopackage.remove_layer(self.path, self.timml_name) + geopackage.remove_layer(self.path, self.ttim_name) + + def on_transient_changed(self, transient: bool): + if len(self.transient_columns) == 0: + return + + config = self.timml_layer.attributeTableConfig() + columns = config.columns() + + for i, column in enumerate(columns): + if column.name in self.transient_columns: + config.setColumnHidden(i, not transient) + + self.timml_layer.setAttributeTableConfig(config) + return + + @staticmethod + def transient_input( + row, all_timeseries: Dict[str, Any], variable: str + ) -> Tuple[List[Any], Set[float]]: + timeseries_id = row["timeseries_id"] + row_start = row["time_start"] + row_end = row["time_end"] + steady_value = row[variable] + transient_value = row[f"{variable}_transient"] + + start_and_stop = ( + row_start is not None + and row_end is not None + and transient_value is not None + ) + + if timeseries_id is not None: + timeseries = all_timeseries[timeseries_id] + return [ + (time, value - steady_value) + for time, value in zip(timeseries["time_start"], timeseries[variable]) + ], set(timeseries["time_start"]) + elif start_and_stop: + return [ + (row_start, transient_value - steady_value), + (row_end, 0.0), + ], {row_end} + else: + return [(0.0, 0.0)], {0.0} + + def check_ttim_columns(self): + return self._check_table_columns( + attributes=self.ttim_attributes, layer=self.ttim_layer + ) + + def to_ttim(self, other): + missing = self.check_ttim_columns() + if missing: + return ElementExtraction(errors=missing) + + other = other.copy() # avoid side-effects + timeseries = self.table_to_dict(self.ttim_layer) + if timeseries: + errors = self.schema.validate_timeseries( + name=self.ttim_layer.name(), data=timeseries + ) + if errors: + return ElementExtraction(errors=errors) + + grouped = self.groupby(timeseries, "timeseries_id") + # Store timeseries in other dict for validation. + other["ttim timeseries IDs"] = set(timeseries["timeseries_id"]) + + else: + grouped = {} + other["ttim timeseries IDs"] = {None} + + data = self.table_to_records(self.timml_layer) + errors = self.schema.validate_ttim( + name=self.timml_layer.name(), data=data, other=other + ) + if errors: + return ElementExtraction(errors=errors) + + elements = [] + times = set() + for row in data: + row_data, row_times = self.process_ttim_row(row, grouped) + elements.append(row_data) + times.update(row_times) + + return ElementExtraction(data=elements, times=times) + + +class AssociatedElement(Element, abc.ABC): + """ + Abstract class for elements that require associated tables such as + Inhomogenities. + """ + + def __init__(self, path: str, name: str): + self._initialize_default(path, name) + self.timml_name = f"timml {self.element_type}:{name}" + self.assoc_name = f"timml {self.element_type} Properties:{name}" + + def create_assoc_layer(self, crs: Any): + self.assoc_layer = self.create_layer( + crs=crs, + geometry_type="No Geometry", + name=self.assoc_name, + attributes=self.assoc_attributes, + ) + + def assoc_layer_from_geopackage(self): + self.assoc_layer = QgsVectorLayer( + f"{self.path}|layername={self.assoc_name}", + self.assoc_name, + ) + + def write(self): + self.timml_layer = geopackage.write_layer( + self.path, self.timml_layer, self.timml_name + ) + self.assoc_layer = geopackage.write_layer( + self.path, self.assoc_layer, self.assoc_name + ) + self.set_defaults() + + def remove_from_geopackage(self): + geopackage.remove_layer(self.path, self.timml_name) + geopackage.remove_layer(self.path, self.assoc_name) + + def to_timml(self, other) -> ElementExtraction: + missing = self._check_table_columns( + attributes=self.timml_attributes, layer=self.timml_layer + ) + if missing: + return ElementExtraction(errors=missing) + + properties = self.table_to_dict(self.assoc_layer) + errors = self.assoc_schema.validate_timml( + name=self.assoc_layer.name(), + data=properties, + ) + if errors: + return ElementExtraction(errors=errors) + + other = other.copy() # Avoid side-effects + other["properties inhomogeneity_id"] = list(set(properties["inhomogeneity_id"])) + data = self.table_to_records(self.timml_layer) + errors = self.schema.validate_timml( + name=self.timml_layer.name(), + data=data, + other=other, + ) + if errors: + return ElementExtraction(errors=errors) + + grouped = self.groupby(properties, "inhomogeneity_id") + elements = [self.process_timml_row(row=row, grouped=grouped) for row in data] + return ElementExtraction(data=elements) + + def to_ttim(self, _): + raise NotImplementedError(f"{type(self).__name__} is not supported in TTim.") diff --git a/plugin/qgistim/core/elements/head_line_sink.py b/plugin/qgistim/core/elements/head_line_sink.py new file mode 100644 index 0000000..1bdf709 --- /dev/null +++ b/plugin/qgistim/core/elements/head_line_sink.py @@ -0,0 +1,97 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import BLUE +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import ( + AllOrNone, + AllRequired, + Membership, + NotBoth, + Optional, + Positive, + Required, +) + + +class HeadLineSinkSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "head": Required(), + "resistance": Required(Positive()), + "width": Required(Positive()), + "order": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + ttim_consistency_schemata = ( + AllOrNone("time_start", "time_end", "head_transient"), + NotBoth("time_start", "timeseries_id"), + ) + ttim_schemata = { + "time_start": Optional(Positive()), + "time_end": Optional(Positive()), + "timeseries_id": Optional(Membership("ttim timeseries IDs")), + } + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time_start": AllRequired(Positive()), + "head": AllRequired(), + } + + +class HeadLineSink(TransientElement): + element_type = "Head Line Sink" + geometry_type = "Linestring" + timml_attributes = ( + QgsField("head", QVariant.Double), + QgsField("resistance", QVariant.Double), + QgsField("width", QVariant.Double), + QgsField("order", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + QgsField("time_start", QVariant.Double), + QgsField("time_end", QVariant.Double), + QgsField("head_transient", QVariant.Double), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time_start", QVariant.Double), + QgsField("head", QVariant.Double), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + } + transient_columns = ( + "time_start", + "time_end", + "head_transient", + "timeseries_id", + ) + schema = HeadLineSinkSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.line_renderer(color=BLUE, width="0.75") + + def process_timml_row(self, row, other=None): + return { + "xy": self.linestring_xy(row), + "hls": row["head"], + "res": row["resistance"], + "wh": row["width"], + "order": row["order"], + "layers": row["layer"], + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + tsandh, times = self.transient_input(row, grouped, "head") + return { + "xy": self.linestring_xy(row), + "tsandh": tsandh, + "res": row["resistance"], + "wh": row["width"], + "layers": row["layer"], + "label": row["label"], + }, times diff --git a/plugin/qgistim/core/elements/headwell.py b/plugin/qgistim/core/elements/headwell.py new file mode 100644 index 0000000..22f0e82 --- /dev/null +++ b/plugin/qgistim/core/elements/headwell.py @@ -0,0 +1,139 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from PyQt5.QtGui import QColor +from qgis.core import ( + QgsArrowSymbolLayer, + QgsDefaultValue, + QgsField, + QgsLineSymbol, + QgsSingleSymbolRenderer, +) +from qgistim.core.elements.colors import BLUE +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import ( + AllOrNone, + AllRequired, + Membership, + NotBoth, + Optional, + Positive, + Required, +) + + +class HeadWellSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "head": Required(), + "radius": Required(Positive()), + "resistance": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + ttim_schemata = { + "time_start": Optional(Positive()), + "time_end": Optional(Positive()), + "timeseries_id": Optional(Membership("ttim timeseries IDs")), + } + ttim_consistency_schemata = ( + AllOrNone(("time_start", "time_end", "head_transient")), + NotBoth("time_start", "timeseries_id"), + ) + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time_start": AllRequired(Positive()), + "head": AllRequired(), + } + + +class HeadWell(TransientElement): + element_type = "Head Well" + geometry_type = "Point" + timml_attributes = ( + QgsField("head", QVariant.Double), + QgsField("radius", QVariant.Double), + QgsField("resistance", QVariant.Double), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + QgsField("time_start", QVariant.Double), + QgsField("time_end", QVariant.Double), + QgsField("head_transient", QVariant.Double), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time_start", QVariant.Double), + QgsField("head", QVariant.Double), + ) + timml_defaults = { + "radius": QgsDefaultValue("0.1"), + "resistance": QgsDefaultValue("0.0"), + } + transient_columns = ( + "time_start", + "time_end", + "head_transient", + "timeseries_id", + ) + schema = HeadWellSchema() + + @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) + return { + "xw": x, + "yw": y, + "hw": row["head"], + "rw": row["radius"], + "res": row["resistance"], + "layers": row["layer"], + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + x, y = self.point_xy(row) + tsandh, times = self.transient_input(row, grouped, "head") + return { + "xw": x, + "yw": y, + "tsandh": tsandh, + "rw": row["radius"], + "res": row["resistance"], + "layers": row["layer"], + "label": row["label"], + }, times + + +class RemoteHeadWell(HeadWell): + element_type = "Remote Head Well" + geometry_type = "Linestring" + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + arrow = QgsArrowSymbolLayer() + red, green, blue, _ = [int(v) for v in BLUE.split(",")] + arrow.setColor(QColor(red, green, blue)) + arrow.setHeadLength(2.5) + symbol = QgsLineSymbol.createSimple({}) + symbol.changeSymbolLayer(0, arrow) + return QgsSingleSymbolRenderer(symbol) + + def process_timml_row(self, row, other=None) -> Dict[str, Any]: + xy = self.linestring_xy(row) + xw, yw = xy[-1] + xc, yc = xy[0] + return { + "xw": xw, + "yw": yw, + "hw": row["head"], + "rw": row["radius"], + "res": row["resistance"], + "layers": row["layer"], + "label": row["label"], + "xc": xc, + "yc": yc, + } diff --git a/plugin/qgistim/core/elements/impermeable_line_doublet.py b/plugin/qgistim/core/elements/impermeable_line_doublet.py new file mode 100644 index 0000000..a61b9d8 --- /dev/null +++ b/plugin/qgistim/core/elements/impermeable_line_doublet.py @@ -0,0 +1,56 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import RED +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import Membership, Positive, Required + + +class ImpermeableLineDoubletSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "order": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + + +class ImpermeableLineDoublet(Element): + element_type = "Impermeable Line Doublet" + geometry_type = "Linestring" + timml_attributes = ( + QgsField("order", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + } + schema = ImpermeableLineDoubletSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.line_renderer(color=RED, width="0.75") + + def process_timml_row(self, row, other=None): + return { + "xy": self.linestring_xy(row), + "layers": row["layer"], + "order": row["order"], + "label": row["label"], + } + + def to_ttim(self, other): + # TTim doesn't have an ImpermeableLineDoublet, we need to add "imp" as + # the resistance entry. + _, data = self.to_timml(other) + out = [] + for row in data: + out.append( + { + "xy": row["xy"], + "res": "imp", + "order": row["order"], + "label": row["label"], + } + ) + return None, out diff --git a/plugin/qgistim/core/elements/leaky_building_pit.py b/plugin/qgistim/core/elements/leaky_building_pit.py new file mode 100644 index 0000000..a4be992 --- /dev/null +++ b/plugin/qgistim/core/elements/leaky_building_pit.py @@ -0,0 +1,107 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import RED, TRANSPARENT_RED +from qgistim.core.elements.element import AssociatedElement +from qgistim.core.elements.schemata import RowWiseSchema, TableSchema +from qgistim.core.schemata import ( + AllGreaterEqual, + AllRequired, + AtleastOneTrue, + Membership, + OffsetAllRequired, + OptionalFirstOnly, + Positive, + Range, + Required, + RequiredFirstOnly, + SemiConfined, + StrictlyDecreasing, +) + + +class LeakyBuildingPitSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "inhomogeneity_id": Required(Membership("properties inhomogeneity_id")), + "order": Required(Positive()), + "ndegrees": Required(Positive()), + } + + +class AssociatedLeakyBuildingPitchema(TableSchema): + timml_schemata = { + "inhomogeneity_id": AllRequired(), + "layer": AllRequired(Range()), + "aquifer_top": AllRequired(StrictlyDecreasing()), + "aquifer_bottom": AllRequired(StrictlyDecreasing()), + "aquitard_c": OffsetAllRequired(Positive()), + "aquifer_k": AllRequired(Positive()), + "semiconf_top": OptionalFirstOnly(), + "semiconf_head": OptionalFirstOnly(), + "resistance": RequiredFirstOnly(), + "wall_in_layer": AllRequired(AtleastOneTrue()), + } + timml_consistency_schemata = ( + SemiConfined(), + AllGreaterEqual("aquifer_top", "aquifer_bottom"), + ) + + +class LeakyBuildingPit(AssociatedElement): + element_type = "Leaky Building Pit" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("order", QVariant.Int), + QgsField("ndegrees", QVariant.Int), + ) + assoc_attributes = [ + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("aquifer_top", QVariant.Double), + QgsField("aquifer_bottom", QVariant.Double), + QgsField("aquitard_c", QVariant.Double), + QgsField("aquifer_k", QVariant.Double), + QgsField("semiconf_top", QVariant.Double), + QgsField("semiconf_head", QVariant.Double), + QgsField("resistance", QVariant.Double), + QgsField("wall_in_layer", QVariant.Bool), + QgsField("aquitard_npor", QVariant.Double), + QgsField("aquifer_npor", QVariant.Double), + ] + timml_defaults = { + "order": QgsDefaultValue("4"), + "ndegrees": QgsDefaultValue("6"), + "inhomogeneity_id": QgsDefaultValue("1"), + } + assoc_defaults = { + "inhomogeneity_id": QgsDefaultValue("1"), + "wall_in_layer": QgsDefaultValue("false"), + } + schema = LeakyBuildingPitSchema() + assoc_schema = AssociatedLeakyBuildingPitchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.polygon_renderer( + color=TRANSPARENT_RED, + color_border=RED, + width_border="0.75", + outline_style="dash", + ) + + def process_timml_row(self, row: Dict[str, Any], grouped: Dict[int, Any]): + inhom_id = row["inhomogeneity_id"] + raw_data = grouped[inhom_id] + aquifer_data = self.aquifer_data(raw_data, transient=False) + layers = [i for i, active in enumerate(raw_data["wall_in_layer"]) if active] + return { + "xy": self.polygon_xy(row), + "order": row["order"], + "ndeg": row["ndegrees"], + "layers": layers, + "res": raw_data["resistance"][0], + **aquifer_data, + } diff --git a/plugin/qgistim/core/elements/leaky_line_doublet.py b/plugin/qgistim/core/elements/leaky_line_doublet.py new file mode 100644 index 0000000..936c8f4 --- /dev/null +++ b/plugin/qgistim/core/elements/leaky_line_doublet.py @@ -0,0 +1,43 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import RED +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import Membership, Positive, Required + + +class LeakyLineDoubletSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "resistance": Required(Positive()), + "order": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + + +class LeakyLineDoublet(Element): + element_type = "Leaky Line Doublet" + geometry_type = "Linestring" + timml_attributes = ( + QgsField("resistance", QVariant.Double), + QgsField("order", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + } + schema = LeakyLineDoubletSchema() + + @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 { + "xy": self.linestring_xy(row), + "res": row["resistance"], + "layers": row["layer"], + "order": row["order"], + "label": row["label"], + } diff --git a/plugin/qgistim/core/elements/line_sink_ditch.py b/plugin/qgistim/core/elements/line_sink_ditch.py new file mode 100644 index 0000000..83d8daa --- /dev/null +++ b/plugin/qgistim/core/elements/line_sink_ditch.py @@ -0,0 +1,97 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import GREEN +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import ( + AllOrNone, + AllRequired, + Membership, + NotBoth, + Optional, + Positive, + Required, +) + + +class LineSinkDitchSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "discharge": Required(), + "resistance": Required(Positive()), + "width": Required(Positive()), + "order": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + ttim_consistency_schemata = ( + AllOrNone("time_start", "time_end", "discharge_transient"), + NotBoth("time_start", "timeseries_id"), + ) + ttim_schemata = { + "time_start": Optional(Positive()), + "time_end": Optional(Positive()), + "timeseries_id": Optional(Membership("ttim timeseries IDs")), + } + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time_start": AllRequired(Positive()), + "discharge": AllRequired(), + } + + +class LineSinkDitch(TransientElement): + element_type = "Line Sink Ditch" + geometry_type = "Linestring" + timml_attributes = ( + QgsField("discharge", QVariant.Double), + QgsField("resistance", QVariant.Double), + QgsField("width", QVariant.Double), + QgsField("order", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + QgsField("time_start", QVariant.Double), + QgsField("time_end", QVariant.Double), + QgsField("discharge_transient", QVariant.Double), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time_start", QVariant.Double), + QgsField("discharge", QVariant.Double), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + } + transient_columns = ( + "time_start", + "time_end", + "discharge_transient", + "timeseries_id", + ) + schema = LineSinkDitchSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.line_renderer(color=GREEN, width="0.75") + + def process_timml_row(self, row, other=None): + return { + "xy": self.linestring_xy(row), + "Qls": row["discharge"], + "res": row["resistance"], + "wh": row["width"], + "order": row["order"], + "layers": row["layer"], + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + tsandQ, times = self.transient_input(row, grouped, "discharge") + return { + "xy": self.linestring_xy(row), + "tsandQ": tsandQ, + "res": row["resistance"], + "wh": row["width"], + "layers": row["layer"], + "label": row["label"], + }, times diff --git a/plugin/qgistim/core/elements/observation.py b/plugin/qgistim/core/elements/observation.py new file mode 100644 index 0000000..7fefaef --- /dev/null +++ b/plugin/qgistim/core/elements/observation.py @@ -0,0 +1,59 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import LIGHT_BLUE +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import AllRequired, Positive, Required + + +class HeadObservationSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + } + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time": AllRequired(Positive()), + } + + +class HeadObservation(TransientElement): + element_type = "Head Observation" + geometry_type = "Point" + timml_attributes = ( + QgsField("label", QVariant.String), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time", QVariant.Double), + ) + timml_defaults = { + "timeseries_id": QgsDefaultValue("1"), + } + ttim_defaults = { + "timeseries_id": QgsDefaultValue("1"), + } + transient_columns = ("timeseries_id",) + schema = HeadObservationSchema() + + @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) + return { + "x": x, + "y": y, + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + x, y = self.point_xy(row) + times = grouped[row["timeseries_id"]]["time"] + return { + "x": x, + "y": y, + "t": times, + "label": row["label"], + }, times diff --git a/plugin/qgistim/core/elements/polygon_area_sink.py b/plugin/qgistim/core/elements/polygon_area_sink.py new file mode 100644 index 0000000..47e016f --- /dev/null +++ b/plugin/qgistim/core/elements/polygon_area_sink.py @@ -0,0 +1,52 @@ +from copy import deepcopy + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import GREEN, TRANSPARENT_GREEN +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import Positive, Required + + +class PolygonAreaSinkSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "rate": Required(), + "order": Required(Positive()), + "ndegrees": Required(Positive()), + } + + +class PolygonAreaSink(Element): + element_type = "Polygon Area Sink" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("rate", QVariant.Double), + QgsField("order", QVariant.Int), + QgsField("ndegrees", QVariant.Int), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + "ndegrees": QgsDefaultValue("6"), + } + schema = PolygonAreaSinkSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.polygon_renderer( + color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75" + ) + + def process_timml_row(self, row, other): + raw_data = deepcopy(other["global_aquifer"]) + raw_data["aquitard_c"][0] = None + raw_data["semiconf_top"][0] = None + raw_data["semiconf_head"][0] = None + aquifer_data = self.aquifer_data(raw_data, transient=False) + return { + "xy": self.polygon_xy(row), + "order": row["order"], + "ndeg": row["ndegrees"], + "N": row["rate"], + **aquifer_data, + } diff --git a/plugin/qgistim/core/elements/polygon_inhomogeneity.py b/plugin/qgistim/core/elements/polygon_inhomogeneity.py new file mode 100644 index 0000000..4b7bddd --- /dev/null +++ b/plugin/qgistim/core/elements/polygon_inhomogeneity.py @@ -0,0 +1,96 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import GREY, TRANSPARENT_GREY +from qgistim.core.elements.element import AssociatedElement +from qgistim.core.elements.schemata import RowWiseSchema, TableSchema +from qgistim.core.schemata import ( + AllGreaterEqual, + AllRequired, + Membership, + OffsetAllRequired, + OptionalFirstOnly, + Positive, + Range, + Required, + SemiConfined, + StrictlyDecreasing, +) + + +class PolygonInhomogeneitySchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "inhomogeneity_id": Required(Membership("properties inhomogeneity_id")), + "order": Required(Positive()), + "ndegrees": Required(Positive()), + } + + +class AssociatedPolygonInhomogeneitySchema(TableSchema): + timml_schemata = { + "inhomogeneity_id": AllRequired(), + "layer": AllRequired(Range()), + "aquifer_top": AllRequired(StrictlyDecreasing()), + "aquifer_bottom": AllRequired(StrictlyDecreasing()), + "aquitard_c": OffsetAllRequired(Positive()), + "aquifer_k": AllRequired(Positive()), + "semiconf_top": OptionalFirstOnly(), + "semiconf_head": OptionalFirstOnly(), + "rate": OptionalFirstOnly(), + } + timml_consistency_schemata = ( + SemiConfined(), + AllGreaterEqual("aquifer_top", "aquifer_bottom"), + ) + + +class PolygonInhomogeneity(AssociatedElement): + element_type = "Polygon Inhomogeneity" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("order", QVariant.Int), + QgsField("ndegrees", QVariant.Int), + ) + assoc_attributes = [ + QgsField("inhomogeneity_id", QVariant.Int), + QgsField("layer", QVariant.Int), + QgsField("aquifer_top", QVariant.Double), + QgsField("aquifer_bottom", QVariant.Double), + QgsField("aquitard_c", QVariant.Double), + QgsField("aquifer_k", QVariant.Double), + QgsField("semiconf_top", QVariant.Double), + QgsField("semiconf_head", QVariant.Double), + QgsField("rate", QVariant.Double), + QgsField("aquitard_npor", QVariant.Double), + QgsField("aquifer_npor", QVariant.Double), + ] + timml_defaults = { + "order": QgsDefaultValue("4"), + "ndegrees": QgsDefaultValue("6"), + "inhomogeneity_id": QgsDefaultValue("1"), + } + assoc_defaults = { + "inhomogeneity_id": QgsDefaultValue("1"), + } + schema = PolygonInhomogeneitySchema() + assoc_schema = AssociatedPolygonInhomogeneitySchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.polygon_renderer( + color=TRANSPARENT_GREY, color_border=GREY, width_border="0.75" + ) + + def process_timml_row(self, row: Dict[str, Any], grouped: Dict[int, Any]): + inhom_id = row["inhomogeneity_id"] + raw_data = grouped[inhom_id] + aquifer_data = self.aquifer_data(raw_data, transient=False) + return { + "xy": self.polygon_xy(row), + "order": row["order"], + "ndeg": row["ndegrees"], + **aquifer_data, + } diff --git a/plugin/qgistim/core/elements/polygon_semi_confined_top.py b/plugin/qgistim/core/elements/polygon_semi_confined_top.py new file mode 100644 index 0000000..f5aa018 --- /dev/null +++ b/plugin/qgistim/core/elements/polygon_semi_confined_top.py @@ -0,0 +1,55 @@ +from copy import deepcopy + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import BLUE, TRANSPARENT_BLUE +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import Positive, Required + + +class PolygonSemiConfinedTopSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "aquitard_c": Required(Positive()), + "semiconf_top": Required(), + "semiconf_head": Required(), + "order": Required(Positive()), + "ndegrees": Required(Positive()), + } + + +class PolygonSemiConfinedTop(Element): + element_type = "Polygon Semi-Confined Top" + geometry_type = "Polygon" + timml_attributes = ( + QgsField("aquitard_c", QVariant.Double), + QgsField("semiconf_top", QVariant.Double), + QgsField("semiconf_head", QVariant.Double), + QgsField("order", QVariant.Int), + QgsField("ndegrees", QVariant.Int), + ) + timml_defaults = { + "order": QgsDefaultValue("4"), + "ndegrees": QgsDefaultValue("6"), + } + schema = PolygonSemiConfinedTopSchema() + + @classmethod + def renderer(cls) -> QgsSingleSymbolRenderer: + return cls.polygon_renderer( + color=TRANSPARENT_BLUE, color_border=BLUE, width_border="0.75" + ) + + def process_timml_row(self, row, other): + raw_data = deepcopy(other["global_aquifer"]) + raw_data["aquitard_c"][0] = row["aquitard_c"] + raw_data["semiconf_top"][0] = row["semiconf_top"] + raw_data["semiconf_head"][0] = row["semiconf_head"] + aquifer_data = self.aquifer_data(raw_data, transient=False) + return { + "xy": self.polygon_xy(row), + "order": row["order"], + "ndeg": row["ndegrees"], + **aquifer_data, + } diff --git a/plugin/qgistim/core/elements/schemata.py b/plugin/qgistim/core/elements/schemata.py new file mode 100644 index 0000000..b46aba0 --- /dev/null +++ b/plugin/qgistim/core/elements/schemata.py @@ -0,0 +1,134 @@ +import abc +from collections import defaultdict +from typing import Any, Dict, List, NamedTuple, Optional, Tuple, Union + +from qgistim.core.schemata import ( + ConsistencySchema, + IterableSchemaContainer, + SchemaContainer, +) + + +class ValidationData(NamedTuple): + schemata: Dict[str, SchemaContainer] + consistency_schemata: Tuple[ConsistencySchema] + name: str + data: Dict[str, Any] + other: Optional[Dict[str, Any]] = None + + +class SchemaBase(abc.ABC): + # TODO: check for presence of columns + timml_schemata: Dict[str, Union[SchemaContainer, IterableSchemaContainer]] = {} + timml_consistency_schemata: Tuple[ConsistencySchema] = () + ttim_schemata: Dict[str, Union[SchemaContainer, IterableSchemaContainer]] = {} + ttim_consistency_schemata: Tuple[ConsistencySchema] = () + timeseries_schemata: Dict[str, Union[SchemaContainer, IterableSchemaContainer]] = {} + + @staticmethod + def _validate_table(vd: ValidationData) -> Dict[str, List]: + errors = defaultdict(list) + for variable, schema in vd.schemata.items(): + _errors = schema.validate(vd.data[variable], vd.other) + if _errors: + errors[variable].extend(_errors) + + # The consistency schema rely on the row input being valid. + # Hence, they are tested second. + if not errors: + for schema in vd.consistency_schemata: + _error = schema.validate(vd.data, vd.other) + if _error: + errors["Table:"].append(_error) + + return errors + + @classmethod + def validate_timeseries( + cls, name: str, data: Dict[str, Any], other=None + ) -> Dict[str, List]: + vd = ValidationData(cls.timeseries_schemata, (), name, data, other) + return cls._validate_table(vd) + + @classmethod + def validate_timml( + cls, name: str, data: Dict[str, Any], other=None + ) -> Dict[str, List]: + vd = ValidationData( + cls.timml_schemata, cls.timml_consistency_schemata, name, data, other + ) + return cls._validate(vd) + + @classmethod + def validate_ttim( + cls, name: str, data: Dict[str, Any], other=None + ) -> Dict[str, List]: + vd = ValidationData( + cls.ttim_schemata, cls.ttim_consistency_schemata, name, data, other + ) + return cls._validate(vd) + + @abc.abstractclassmethod + def _validate(vd: ValidationData) -> Dict[str, List]: + pass + + +class TableSchema(SchemaBase, abc.ABC): + """ + Schema for Tabular data, such as Aquifer properties. + """ + + @classmethod + def _validate( + cls, + vd: ValidationData, + ) -> Dict[str, List]: + return cls._validate_table(vd) + + +class RowWiseSchema(SchemaBase, abc.ABC): + """ + Schema for entries that should be validated row-by-row, such as Wells. + """ + + @staticmethod + def _validate(vd: ValidationData) -> Dict[str, List]: + errors = defaultdict(list) + + for i, row in enumerate(vd.data): + row_errors = defaultdict(list) + + for variable, schema in vd.schemata.items(): + _errors = schema.validate(row[variable], vd.other) + if _errors: + row_errors[variable].extend(_errors) + + # Skip consistency tests if the individual values are not good. + if not row_errors: + for schema in vd.consistency_schemata: + _error = schema.validate(row, vd.other) + if _error: + row_errors["Row:"].append(_error) + + if row_errors: + errors[f"Row {i + 1}:"] = row_errors + + return errors + + +class SingleRowSchema(RowWiseSchema, abc.ABC): + """ + Schema for entries that should contain only one row, which should be + validated as a row, such as Constant, Domain, Uniform Flow. + """ + + @staticmethod + def _validate(vd: ValidationData) -> Dict[str, List]: + nrow = len(vd.data) + if nrow > 1: + return { + vd.name: [ + "Table must contain a single row. " f"Table contains {nrow} rows." + ] + } + return RowWiseSchema._validate(vd) diff --git a/plugin/qgistim/core/elements/uniform_flow.py b/plugin/qgistim/core/elements/uniform_flow.py new file mode 100644 index 0000000..6993f3a --- /dev/null +++ b/plugin/qgistim/core/elements/uniform_flow.py @@ -0,0 +1,30 @@ +from PyQt5.QtCore import QVariant +from qgis.core import QgsField +from qgistim.core.elements.element import Element +from qgistim.core.elements.schemata import SingleRowSchema +from qgistim.core.schemata import Required + + +class UniformFlowSchema(SingleRowSchema): + timml_schemata = { + "slope": Required(), + "angle": Required(), + } + + +class UniformFlow(Element): + element_type = "Uniform Flow" + geometry_type = "No geometry" + timml_attributes = ( + QgsField("slope", QVariant.Double), + QgsField("angle", QVariant.Double), + QgsField("label", QVariant.String), + ) + schema = UniformFlowSchema() + + def process_timml_row(self, row, other=None): + return { + "slope": row["slope"], + "angle": row["angle"], + "label": row["label"], + } diff --git a/plugin/qgistim/core/elements/well.py b/plugin/qgistim/core/elements/well.py new file mode 100644 index 0000000..46aee89 --- /dev/null +++ b/plugin/qgistim/core/elements/well.py @@ -0,0 +1,111 @@ +from typing import Any, Dict + +from PyQt5.QtCore import QVariant +from qgis.core import QgsDefaultValue, QgsField, QgsSingleSymbolRenderer +from qgistim.core.elements.colors import GREEN +from qgistim.core.elements.element import TransientElement +from qgistim.core.elements.schemata import RowWiseSchema +from qgistim.core.schemata import ( + AllOrNone, + AllRequired, + Membership, + NotBoth, + Optional, + Positive, + Required, +) + + +class WellSchema(RowWiseSchema): + timml_schemata = { + "geometry": Required(), + "discharge": Required(), + "radius": Required(Positive()), + "resistance": Required(Positive()), + "layer": Required(Membership("aquifer layers")), + } + ttim_schemata = { + "caisson_radius": Required(Positive), + "slug": Required(), + "time_start": Optional(Positive()), + "time_end": Optional(Positive()), + "timeseries_id": Optional(Membership("ttim timeseries IDs")), + } + ttim_consistency_schemata = ( + AllOrNone(("time_start", "time_end", "discharge_transient")), + NotBoth("time_start", "timeseries_id"), + ) + timeseries_schemata = { + "timeseries_id": AllRequired(), + "time_start": AllRequired(Positive()), + "discharge": AllRequired(), + } + + +class Well(TransientElement): + element_type = "Well" + geometry_type = "Point" + timml_attributes = ( + QgsField("discharge", QVariant.Double), + QgsField("radius", QVariant.Double), + QgsField("resistance", QVariant.Double), + QgsField("layer", QVariant.Int), + QgsField("label", QVariant.String), + QgsField("time_start", QVariant.Double), + QgsField("time_end", QVariant.Double), + QgsField("discharge_transient", QVariant.Double), + QgsField("caisson_radius", QVariant.Double), + QgsField("slug", QVariant.Bool), + QgsField("timeseries_id", QVariant.Int), + ) + ttim_attributes = ( + QgsField("timeseries_id", QVariant.Int), + QgsField("time_start", QVariant.Double), + QgsField("discharge", QVariant.Double), + ) + timml_defaults = { + "radius": QgsDefaultValue("0.1"), + "resistance": QgsDefaultValue("0.0"), + "caisson_radius": QgsDefaultValue("0.0"), + "slug": QgsDefaultValue("True"), + } + transient_columns = ( + "time_start", + "time_end", + "discharge_transient", + "caisson_radius", + "slug", + "timeseries_id", + ) + schema = WellSchema() + + @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) + return { + "xw": x, + "yw": y, + "Qw": row["discharge"], + "rw": row["radius"], + "res": row["resistance"], + "layers": row["layer"], + "label": row["label"], + } + + def process_ttim_row(self, row, grouped): + x, y = self.point_xy(row) + tsandQ, times = self.transient_input(row, grouped, "discharge") + return { + "xw": x, + "yw": y, + "tsandQ": tsandQ, + "rw": row["radius"], + "res": row["resistance"], + "layers": row["layer"], + "label": row["label"], + "rc": row["caisson_radius"], + "wbstype": "slug" if row["slug"] else "pumping", + }, times diff --git a/plugin/qgistim/core/extractor.py b/plugin/qgistim/core/extractor.py new file mode 100644 index 0000000..482c12b --- /dev/null +++ b/plugin/qgistim/core/extractor.py @@ -0,0 +1,67 @@ +import abc +from collections import defaultdict +from typing import Any, Dict, List, Tuple + +from qgis.core import NULL, QgsVectorLayer + + +def remove_zero_length(geometry): + # This removes repeated vertices + return list(dict.fromkeys(geometry)) + + +class ExtractorMixin(abc.ABC): + """ + Mixin class to extract all data from QgsVectorLayers. + """ + + @staticmethod + def argsort(seq): + return sorted(range(len(seq)), key=seq.__getitem__) + + @staticmethod + def extract_coordinates(feature): + geometry = feature.geometry() + coordinates = [] + for vertex in geometry.vertices(): + coordinates.append((vertex.x(), vertex.y())) + centroid = geometry.centroid().asPoint() + return (centroid.x(), centroid.y()), coordinates + + @classmethod + def table_to_records(cls, layer: QgsVectorLayer) -> List[Dict[str, Any]]: + features = [] + for feature in layer.getFeatures(): + data = feature.attributeMap() + for key, value in data.items(): + if value == NULL: + data[key] = None + geomtype = layer.geometryType() + # Skip if no geometry is present. + if geomtype != geomtype.Null: + data["centroid"], data["geometry"] = cls.extract_coordinates(feature) + features.append(data) + return features + + def table_to_dict(cls, layer: QgsVectorLayer) -> Dict[str, Any]: + features = defaultdict(list) + for feature in layer.getFeatures(): + for key, value in feature.attributeMap().items(): + if value == NULL: + features[key].append(None) + else: + features[key].append(value) + return features + + @staticmethod + def point_xy(row) -> Tuple[List[float], List[float]]: + point = row["geometry"][0] + return point[0], point[1] + + @staticmethod + def linestring_xy(row) -> List: + return remove_zero_length(row["geometry"]) + + @staticmethod + def polygon_xy(row) -> List: + return remove_zero_length(row["geometry"]) diff --git a/plugin/qgistim/core/formatting.py b/plugin/qgistim/core/formatting.py new file mode 100644 index 0000000..a9ebc44 --- /dev/null +++ b/plugin/qgistim/core/formatting.py @@ -0,0 +1,324 @@ +""" +Format the content of a collection of dictionaries into a Python script or a +dictionary to be serialized to JSON. +""" + +import pprint +import re +import textwrap +from typing import Any, Dict, Tuple, Union + +import numpy as np + +TIMML_MAPPING = { + "Constant": "Constant", + "Uniform Flow": "Uflow", + "Circular Area Sink": "CircAreaSink", + "Well": "Well", + "Head Well": "HeadWell", + "Remote Head Well": "HeadWell", + "Polygon Inhomogeneity": "PolygonInhomMaq", + "Polygon Area Sink": "PolygonInhomMaq", + "Polygon Semi-Confined Top": "PolygonInhomMaq", + "Head Line Sink": "HeadLineSinkString", + "Line Sink Ditch": "LineSinkDitchString", + "Leaky Line Doublet": "LeakyLineDoubletString", + "Impermeable Line Doublet": "ImpLineDoubletString", + "Building Pit": "BuildingPit", + "Leaky Building Pit": "LeakyBuildingPit", + "Head Observation": "Head Observation", + "Discharge Observation": "Discharge Observation", +} +# In TTim, a constant or uniform flow may be added, but they have no effect on +# the transient superposed result. +TTIM_MAPPING = { + "Constant": None, + "Uniform Flow": None, + "Circular Area Sink": "CircAreaSink", + "Well": "Well", + "Head Well": "HeadWell", + "Head Line Sink": "HeadLineSinkString", + "Line Sink Ditch": "LineSinkDitchString", + "Leaky Line Doublet": "LeakyLineDoubletString", + "Impermeable Line Doublet": "LeakyLineDoubletString", + "Head Observation": "Head Observation", +} +PREFIX = " " + + +def sanitized(name: str) -> str: + return name.split(":")[-1].replace(" ", "_") + + +def format_kwargs(data: Dict[str, Any]) -> str: + return textwrap.indent( + "\n".join(f"{k}={pprint.pformat(v)}," for k, v in data.items()), prefix=PREFIX + ) + + +def round_spacing(ymin: float, ymax: float) -> float: + """ + Some reasonable defaults for grid spacing. + + We attempt to get around 50 rows in the computed grid, with grid sizes a + multiple of 1.0, 5.0, 50.0, or 500.0. + """ + dy = (ymax - ymin) / 50.0 + if dy > 500.0: + dy = round(dy / 500.0) * 500.0 + elif dy > 50.0: + dy = round(dy / 50.0) * 50.0 + elif dy > 5.0: # round to five + dy = round(dy / 5.0) * 5.0 + elif dy > 1.0: + dy = round(dy) + return dy + + +def round_extent(domain: Dict[str, float], cellsize: float) -> Tuple[float]: + """ + Increases the extent until all sides lie on a coordinate + divisible by cellsize. + + Parameters + ---------- + extent: Tuple[float] + xmin, xmax, ymin, ymax + cellsize: float + Desired cell size of the output head grids + + Returns + ------- + extent: Tuple[float] + xmin, xmax, ymin, ymax + """ + xmin = domain["xmin"] + ymin = domain["ymin"] + xmax = domain["xmax"] + ymax = domain["ymax"] + xmin = np.floor(xmin / cellsize) * cellsize + ymin = np.floor(ymin / cellsize) * cellsize + xmax = np.ceil(xmax / cellsize) * cellsize + ymax = np.ceil(ymax / cellsize) * cellsize + xmin += 0.5 * cellsize + xmax += 0.5 * cellsize + ymax -= 0.5 * cellsize + xmin -= 0.5 * cellsize + return xmin, xmax, ymin, ymax + + +def headgrid_entry(domain: Dict[str, float], cellsize: float) -> Dict[str, float]: + (xmin, xmax, ymin, ymax) = round_extent(domain, cellsize) + return { + "xmin": xmin, + "xmax": xmax, + "ymin": ymin, + "ymax": ymax, + "spacing": cellsize, + "time": domain.get("time"), + } + + +def headgrid_code(domain) -> Tuple[str, str]: + ymin = domain["ymin"] + ymax = domain["ymax"] + dy = round_spacing(ymin, ymax) + (xmin, xmax, ymin, ymax) = round_extent(domain, dy) + xg = textwrap.indent(f"xg=np.arange({xmin}, {xmax}, {dy})", prefix=PREFIX) + yg = textwrap.indent(f"yg=np.arange({ymax}, {ymin}, -{dy})", prefix=PREFIX) + time = domain.get("time") + t = textwrap.indent(f"t={time}", prefix=PREFIX) + return xg, yg, t + + +def elements_and_observations(data, mapping: Dict[str, str], tim: str): + strings = [] + observations = [] + model_string = textwrap.indent(f"model={tim}_model,", prefix=PREFIX) + for layername, element_data in data.items(): + prefix, name = layername.split(":") + plugin_name = re.split("timml |ttim ", prefix)[1] + tim_name = mapping[plugin_name] + if tim_name is None: + continue + + for i, kwargs in enumerate(element_data): + if plugin_name == "Head Observation": + # Should not be added to the model. + # Would result in e.g.: + # observation_piezometer_0 = timml.head( + # x=10.0, + # y=20.0, + # ) + kwargs.pop("label") + observations.append( + f"observation_{sanitized(name)}_{i}={tim}_model.head(\n{format_kwargs(kwargs)}\n)" + ) + else: + # Has to be added to the model. + # Would result in e.g.: + # timml_extraction_0 = timml.Well( + # model=timml_model, + # ... + # ) + kwargs = format_kwargs(kwargs) + strings.append( + f"{tim}_{sanitized(name)}_{i} = {tim}.{tim_name}(\n{model_string}\n{kwargs}\n)" + ) + + return strings, observations + + +def timml_script_content(data: Dict[str, Any]): + data = data.copy() # avoid side-effects + aquifer_data = data.pop("timml Aquifer:Aquifer") + data.pop("timml Domain:Domain") + + strings = [ + "import numpy as np", + "import timml", + "", + f"timml_model = timml.ModelMaq(\n{format_kwargs(aquifer_data)}\n)", + ] + + element_strings, observations = elements_and_observations( + data, TIMML_MAPPING, tim="timml" + ) + strings = strings + element_strings + return strings, observations + + +def timml_script(data: Dict[str, Any]) -> str: + strings, observations = timml_script_content(data) + strings.append("\ntimml_model.solve()\n") + xg, yg, _ = headgrid_code(data["timml Domain:Domain"]) + strings.append(f"head = timml_model.headgrid(\n{xg},\n{yg}\n)") + strings.append("\n") + strings.extend(observations) + return "\n".join(strings) + + +def ttim_script(timml_data: Dict[str, Any], ttim_data: Dict[str, Any]) -> str: + strings, _ = timml_script_content(timml_data) + strings.insert(2, "import ttim") + + data = ttim_data.copy() # avoid side-effects + aquifer_data = data.pop("timml Aquifer:Aquifer") + domain_data = data.pop("timml Domain:Domain") + data.pop("reference_date") + + strings.append( + f"\nttim_model = ttim.ModelMaq(\n{format_kwargs(aquifer_data)}\n{PREFIX}timmlmodel=timml_model,\n)" + ) + + element_strings, observations = elements_and_observations( + data, TTIM_MAPPING, tim="ttim" + ) + strings = strings + element_strings + strings.append("\ntimml_model.solve()\nttim_model.solve()\n") + + if domain_data.get("time"): + xg, yg, t = headgrid_code(domain_data) + strings.append(f"head = ttim_model.headgrid(\n{xg},\n{yg},\n{t}\n)") + strings.append("\n") + + strings.extend(observations) + return "\n".join(strings) + + +def data_to_script( + timml_data: Dict[str, Any], + ttim_data: Union[Dict[str, Any], None], +) -> str: + if ttim_data is None: + return timml_script(timml_data) + else: + return ttim_script(timml_data, ttim_data) + + +def json_elements_and_observations(data, mapping: Dict[str, str]): + aquifer_data = data.pop("timml Aquifer:Aquifer") + + observations = {} + tim_data = {"ModelMaq": aquifer_data} + for layername, element_data in data.items(): + prefix, name = layername.split(":") + plugin_name = re.split("timml |ttim ", prefix)[1] + tim_name = mapping[plugin_name] + if tim_name is None: + continue + + entry = {"type": tim_name, "name": name, "data": element_data} + if tim_name == "Head Observation": + observations[layername] = entry + else: + tim_data[layername] = entry + + return tim_data, observations + + +def timml_json( + timml_data: Dict[str, Any], + cellsize: float, + output_options: Dict[str, bool], +) -> Dict[str, Any]: + """ + Take the data and add: + + * the TimML type + * the layer name + + Parameters + ---------- + data: Dict[str, Any] + + Returns + ------- + json_data: Dict[str, Any] + Data ready to dump to JSON. + """ + # Process TimML elements + data = timml_data.copy() # avoid side-effects + domain_data = data.pop("timml Domain:Domain") + timml_json, observations = json_elements_and_observations( + data, mapping=TIMML_MAPPING + ) + json_data = { + "timml": timml_json, + "output_options": output_options, + "headgrid": headgrid_entry(domain_data, cellsize), + "observations": observations, + } + return json_data + + +def ttim_json( + timml_data: Dict[str, Any], + ttim_data: Dict[str, Any], + cellsize: float, + output_options: Dict[str, bool], +) -> Dict[str, Any]: + json_data = timml_json(timml_data, cellsize, output_options) + + data = ttim_data.copy() + domain_data = data.pop("timml Domain:Domain") + reference_date = data.pop("reference_date") + ttim_json, observations = json_elements_and_observations(data, mapping=TTIM_MAPPING) + + json_data["ttim"] = ttim_json + json_data["headgrid"] = headgrid_entry(domain_data, cellsize) + json_data["reference_date"] = reference_date + json_data["observations"] = observations + return json_data + + +def data_to_json( + timml_data: Dict[str, Any], + ttim_data: Union[Dict[str, Any], None], + cellsize: float, + output_options: Dict[str, bool], +) -> Dict[str, Any]: + if ttim_data is None: + return timml_json(timml_data, cellsize, output_options) + else: + return ttim_json(timml_data, ttim_data, cellsize, output_options) diff --git a/plugin/qgistim/core/processing.py b/plugin/qgistim/core/processing.py index b885973..32e540d 100644 --- a/plugin/qgistim/core/processing.py +++ b/plugin/qgistim/core/processing.py @@ -5,6 +5,7 @@ into line contours. """ import datetime +from pathlib import Path from typing import NamedTuple import numpy as np @@ -22,6 +23,7 @@ QgsVectorLayer, QgsVectorLayerTemporalProperties, ) +from qgistim.core import geopackage def raster_steady_contours( @@ -78,7 +80,6 @@ def steady_contours( stop: float, step: float, ) -> QgsVectorLayer: - contourer = QgsMeshContours(layer) # Collect contours from mesh layer @@ -94,7 +95,7 @@ def steady_contours( feature_data.append(SteadyContourData(geom, value)) # Setup output layer - contour_layer = QgsVectorLayer("Linestring", f"contours-{name}", "memory") + contour_layer = QgsVectorLayer("Linestring", name, "memory") contour_layer.setCrs(layer.crs()) provider = contour_layer.dataProvider() provider.addAttributes( @@ -110,7 +111,7 @@ def steady_contours( f.setGeometry(item.geometry) # Make sure to convert to the appropriate Qt types # e.g. no numpy floats allowed - f.setAttributes([float(item.head)]) + f.setAttributes([round(float(item.head), 3)]) provider.addFeature(f) contour_layer.updateExtents() @@ -131,7 +132,7 @@ def transient_contours( contourer = QgsMeshContours(layer) # Setup output layer - contour_layer = QgsVectorLayer("Linestring", f"contours-{name}", "memory") + contour_layer = QgsVectorLayer("Linestring", name, "memory") contour_layer.setCrs(layer.crs()) provider = contour_layer.dataProvider() provider.addAttributes( @@ -169,7 +170,7 @@ def transient_contours( end_dates = { a: b - datetime.timedelta(minutes=1) for a, b in zip(dates[:-1], dates[1:]) } - end_dates[dates[-1]] = dates[-1] + datetime.timedelta(days=1) + end_dates[dates[-1]] = dates[-1] + datetime.timedelta(minutes=1) # Add items to layer for item in feature_data: @@ -179,7 +180,7 @@ def transient_contours( # e.g. no numpy floats allowed f.setAttributes( [ - float(item.head), + float(round(item.head, 3)), QDateTime(item.datetime), QDateTime(end_dates[item.datetime]), ] @@ -187,19 +188,24 @@ def transient_contours( provider.addFeature(f) contour_layer.updateExtents() - # Set the temporal properties - temporal_properties = contour_layer.temporalProperties() - temporal_properties.setStartField("datetime_start") - temporal_properties.setEndField("datetime_end") - temporal_properties.setMode( - QgsVectorLayerTemporalProperties.ModeFeatureDateTimeStartAndEndFromFields - ) - temporal_properties.setIsActive(True) - return contour_layer +def set_temporal_properties(layer: QgsVectorLayer) -> None: + fields = [field.name() for field in layer.fields()] + if ("datetime_start" in fields) and ("datetime_end" in fields): + temporal_properties = layer.temporalProperties() + temporal_properties.setStartField("datetime_start") + temporal_properties.setEndField("datetime_end") + temporal_properties.setMode( + QgsVectorLayerTemporalProperties.ModeFeatureDateTimeStartAndEndFromFields + ) + temporal_properties.setIsActive(True) + return + + def mesh_contours( + gpkg_path: str, layer: QgsMeshLayer, index: int, name: str, @@ -208,6 +214,16 @@ def mesh_contours( step: float, ) -> QgsVectorLayer: if layer.firstValidTimeStep().isValid(): - return transient_contours(layer, index, name, start, stop, step) + vector_layer = transient_contours(layer, index, name, start, stop, step) else: - return steady_contours(layer, index, name, start, stop, step) + vector_layer = steady_contours(layer, index, name, start, stop, step) + + newfile = not Path(gpkg_path).exists() + written_layer = geopackage.write_layer( + path=gpkg_path, + layer=vector_layer, + layername=name, + newfile=newfile, + ) + set_temporal_properties(written_layer) + return written_layer diff --git a/plugin/qgistim/core/schemata.py b/plugin/qgistim/core/schemata.py new file mode 100644 index 0000000..022add2 --- /dev/null +++ b/plugin/qgistim/core/schemata.py @@ -0,0 +1,337 @@ +import abc +import operator +from typing import List, Sequence, Union + +OPERATORS = { + "<": operator.lt, + "<=": operator.le, + "==": operator.eq, + "!=": operator.ne, + ">=": operator.ge, + ">": operator.gt, +} + + +MaybeError = Union[None, str] +ErrorList = List[str] + + +def format(data) -> str: + return ", ".join(map(str, data)) + + +# Base classes + + +class BaseSchema(abc.ABC): + """Base class for single value.""" + + def __init__(self): + pass + + @abc.abstractmethod + def validate(self, data, other) -> MaybeError: + pass + + def validate_many(self, data, other) -> ErrorList: + errors = [] + for value in data: + error = self.validate(value, other) + if error is not None: + errors.append(error) + + return errors + + +class ConsistencySchema(BaseSchema, abc.ABC): + pass + + +class IterableSchema(abc.ABC): + """Base class for collection of values.""" + + def __init__(self, *schemata): + self.schemata = schemata + + def validate_many(self, data, other) -> ErrorList: + error = self.validate(data, other) + if error: + return [error] + else: + return [] + + +class SchemaContainer(abc.ABC): + def __init__(self, *schemata): + self.schemata = schemata + + @abc.abstractmethod + def validate(self): + pass + + def _validate_schemata(self, data, other=None) -> ErrorList: + errors = [] + for schema in self.schemata: + _error = schema.validate(data, other) + if _error: + errors.append(_error) + return errors + + +class IterableSchemaContainer(abc.ABC): + def __init__(self, *schemata): + self.schemata = schemata + + @abc.abstractmethod + def validate(self): + pass + + def _validate_schemata(self, data, other=None) -> ErrorList: + errors = [] + for schema in self.schemata: + _errors = schema.validate_many(data, other) + if _errors: + errors.extend(_errors) + return errors + + +# Schema containers for a single value. + + +class OptionalFirstOnly(SchemaContainer): + """Exclusively the first value may be provided.""" + + def validate(self, data, other=None) -> ErrorList: + if any(v is not None for v in data[1:]): + return ["Only the first value may be filled in."] + elif data[0] is None: + return [] + else: + return self._validate_schemata(data[0], other) + + +class RequiredFirstOnly(SchemaContainer): + """Exclusively the first value must be provided.""" + + def validate(self, data, other=None) -> ErrorList: + if data[0] is None: + return ["The first value must be filled in."] + elif any(v is not None for v in data[1:]): + return ["Only the first value may be filled in."] + else: + return self._validate_schemata(data[0], other) + + +class Optional(SchemaContainer): + def validate(self, data, other=None) -> ErrorList: + if data is None: + return [] + return self._validate_schemata(data, other) + + +class Required(SchemaContainer): + def validate(self, data, other=None) -> ErrorList: + if data is None: + return ["a value is required."] + return self._validate_schemata(data, other) + + +# SchemataContainer for multiple values. + + +class AllRequired(IterableSchemaContainer): + def validate(self, data, other=None) -> ErrorList: + missing = [i + 1 for i, v in enumerate(data) if v is None] + if missing: + return [f"No values provided at row(s): {format(missing)}"] + return self._validate_schemata(data, other) + + +class OffsetAllRequired(IterableSchemaContainer): + def validate(self, data, other=None) -> ErrorList: + missing = [i + 2 for i, v in enumerate(data[1:]) if v is None] + if missing: + return [f"No values provided at row(s): {format(missing)}"] + if data[0] is None: + return self._validate_schemata(data[1:], other) + else: + return self._validate_schemata(data, other) + + +class AllOptional(IterableSchemaContainer): + def validate(self, data, other=None) -> ErrorList: + missing = [i + 1 for i, v in enumerate(data) if v is None] + if len(missing) == len(data): + return [] + return self._validate_schemata(data, other) + + +# Schemata for a single value. + + +class Positive(BaseSchema): + def validate(self, data, _=None) -> MaybeError: + if data < 0: + return f"Non-positive value: {data}" + return None + + +class AllOrNone(BaseSchema): + def __init__(self, *variables: Sequence[str]): + self.variables = variables + + def validate(self, data, _=None) -> MaybeError: + present = [data.get(v) is not None for v in self.variables] + if any(present) != all(present): + vars = ", ".join(self.variables) + return ( + "Exactly all or none of the following variables must be " + f"provided: {vars}" + ) + return None + + +class NotBoth(BaseSchema): + def __init__(self, x: str, y: str): + self.x = x + self.y = y + + def validate(self, data, _=None) -> MaybeError: + if (data.get(self.x) is not None) and (data.get(self.y) is not None): + return f"Either {self.x} or {self.y} should be provided, not both." + return None + + +class Membership(BaseSchema): + def __init__(self, members_key: str): + self.members_key = members_key + + def validate(self, data, other=None) -> MaybeError: + if data is None: + return None + member_values = other[self.members_key] + if data not in member_values: + return ( + f"Value {data} not found in {self.members_key}: {format(member_values)}" + ) + return None + + +class CircularGeometry(BaseSchema): + def validate(self, data, _=None) -> MaybeError: + coordinates = data + # Compute centroid. + n_vertex = len(data) + x_mean = 0.0 + y_mean = 0.0 + for x, y in coordinates: + x_mean += x / n_vertex + y_mean += y / n_vertex + # Compute distances to centroid. + distances = [(x - x_mean) ** 2 + (y - y_mean) ** 2 for (x, y) in coordinates] + min_distance = min(distances) ** 0.5 + max_distance = max(distances) ** 0.5 + # Accept 1% deviation in squared distance from a circle. + if (max_distance - min_distance) > (0.01 * min_distance): + return "Geometry is not circular." + return None + + +# Schemata for a collection of values. + + +class Range(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + expected = list(range(len(data))) + if not data == expected: + return f"Expected {format(expected)}; received {format(data)}" + return None + + +class Increasing(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + monotonic = all(a <= b for a, b in zip(data, data[1:])) + if not monotonic: + return f"Values are not increasing: {format(data)}" + return None + + +class StrictlyIncreasing(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + monotonic = all(a < b for a, b in zip(data, data[1:])) + if not monotonic: + return f"Values are not strictly increasing (no repeated values): {format(data)}" + return None + + +class Decreasing(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + monotonic = all(a >= b for a, b in zip(data, data[1:])) + if not monotonic: + return f"Values are not decreasing: {format(data)}" + return None + + +class StrictlyDecreasing(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + monotonic = all(a > b for a, b in zip(data, data[1:])) + if not monotonic: + return f"Values are not strictly decreasing (no repeated values): {format(data)}" + return None + + +class AllGreaterEqual(IterableSchema): + def __init__(self, x, y): + self.x = x + self.y = y + + def validate(self, data, _=None) -> MaybeError: + x = data[self.x] + y = data[self.y] + wrong = [i + 1 for i, (a, b) in enumerate(zip(x, y)) if a < b] + if wrong: + return f"{self.x} is not greater or equal to {self.y} at row(s): {format(wrong)}" + return None + + +class AtleastOneTrue(IterableSchema): + def validate(self, data, _=None) -> MaybeError: + if not any(value for value in data): + return "Atleast one row value must be true." + return None + + +# Consistency schemata + + +class SemiConfined(ConsistencySchema): + def validate(self, data, _=None) -> MaybeError: + semiconf_data = { + "aquitard_c": data["aquitard_c"][0], + "semiconf_top": data["semiconf_top"][0], + "semiconf_head": data["semiconf_head"][0], + } + present = [v is not None for v in semiconf_data.values()] + if any(present) != all(present): + variables = format(semiconf_data.keys()) + values = format(semiconf_data.values()) + return ( + "To enable a semi-confined top, the first row must be fully " + f"filled in for {variables}. To disable semi-confined top, none " + f"of the values must be filled in. Found: {values}" + ) + semitop = data["semiconf_top"][0] + if semitop is not None and semitop <= data["aquifer_top"][0]: + return "semiconf_top must be greater than first aquifer_top." + if "rate" in data: + if data["rate"][0] is not None and semitop: + return "A rate cannot be given when a semi-confined is enabled." + return None + + +class SingleRow(ConsistencySchema): + def validate(self, data, _=None) -> MaybeError: + nrow = len(data) + if nrow != 1: + return f"Table must contain one row, found {nrow} rows." + return None diff --git a/plugin/qgistim/core/tim_elements.py b/plugin/qgistim/core/tim_elements.py deleted file mode 100644 index 1020c1a..0000000 --- a/plugin/qgistim/core/tim_elements.py +++ /dev/null @@ -1,925 +0,0 @@ -""" -This module contains the classes to represent the TimML and TTim input layers. - -The classes specify: - -* The (unabbreviated) name -* The type of geometry (No geometry, point, linestring, polygon) -* The required attributes of the attribute table - -They contain logic for setting up: - -* Simple input, requiring a single table, e.g. Uniform Flow or Constant -* Transient input, requiring two tables, one with geometry and steady-state - properties, and one containing optional time series input e.g. Well, Head Line - Sink. -* Associated input, requiring two tables, one with geometry and one with input - for layers, e.g. Polygon Inhomogeneity or Building Pit. - -Each element is (optionally) represented in multiple places: - -* It always lives in a GeoPackage. -* While a geopackage is active within plugin, it is always represented in a - Dataset Tree: the Dataset Tree provides a direct look at the state of the - GeoPackage. In this tree, steady and transient input are on the same row. - Associated input is, to potentially enable transient associated data later - on (like a building pit with changing head top boundary). -* It can be added to the Layers Panel in QGIS. This enables a user to visualize - and edit its data. - -Some elements require specific rendering in QGIS (e.g. no fill polygons), which -are supplied by the `.renderer` property. - -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 -associated with each other is by comparing names. Names must therefore be -unique within a group of the same type of elements. - -Rendering: - -* Fixed discharge (area sink, well, ditch): green -* Fixed head (head well, semi-confined top, head line): blue -* Inhomogeneity: grey -* Constant: star -* Observations: triangle -* Line Doublets: red -* Polygons: Line and Fill same color, Fill color 15% opacity. -""" - -import re -from collections import defaultdict -from functools import partial -from typing import Any, List, Tuple - -from PyQt5.QtCore import QVariant -from PyQt5.QtWidgets import ( - QDialog, - QHBoxLayout, - QLabel, - QLineEdit, - QPushButton, - QVBoxLayout, -) -from qgis.core import ( - QgsDefaultValue, - QgsFeature, - QgsField, - QgsFillSymbol, - QgsGeometry, - QgsLineSymbol, - QgsMarkerSymbol, - QgsPointXY, - QgsSingleSymbolRenderer, - QgsVectorLayer, -) -from qgistim.core import geopackage - -RED = "215,48,39,255" -GREEN = "51,160,44,255" -BLUE = "31,120,180,255" -GREY = "135,135,135,255" -BLACK = "0,0,0,255" -LIGHT_BLUE = "166,206,227,255" - -OPACITY = 51 # 20% -TRANSPARENT = "255,255,255,0" -TRANSPARENT_RED = f"215,48,39,{OPACITY}" -TRANSPARENT_GREEN = f"51,160,44,{OPACITY}" -TRANSPARENT_BLUE = f"31,120,180,{OPACITY}" -TRANSPARENT_GREY = f"135,135,135,{OPACITY}" - - -# These columns are reused by Aquifer and Polygon Inhom, Building pit Aquitards -# are on top of the aquifer, so it comes first Nota bene: the order of these is -# important for hiding and showing the transient columns. QGIS has a bug which -# causes it to show the wrong column if hidden columns appear before shown -# ones. This only affects the attribute table when it has no features. -AQUIFER_ATTRIBUTES = [ - QgsField("layer", QVariant.Int), - QgsField("aquifer_top", QVariant.Double), - QgsField("aquifer_bottom", QVariant.Double), - QgsField("aquitard_c", QVariant.Double), - QgsField("aquifer_k", QVariant.Double), - QgsField("semiconf_top", QVariant.Double), - QgsField("semiconf_head", QVariant.Double), - QgsField("aquitard_s", QVariant.Double), - QgsField("aquifer_s", QVariant.Double), - QgsField("aquitard_npor", QVariant.Double), - QgsField("aquifer_npor", QVariant.Double), -] -INHOM_ATTRIBUTES = [ - QgsField("inhomogeneity_id", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("aquifer_top", QVariant.Double), - QgsField("aquifer_bottom", QVariant.Double), - QgsField("aquitard_c", QVariant.Double), - QgsField("aquifer_k", QVariant.Double), - QgsField("semiconf_top", QVariant.Double), - QgsField("semiconf_head", QVariant.Double), - QgsField("rate", QVariant.Double), - QgsField("aquitard_s", QVariant.Double), - QgsField("aquifer_s", QVariant.Double), - QgsField("aquitard_npor", QVariant.Double), - QgsField("aquifer_npor", QVariant.Double), -] -BUILDING_PIT_ATTRIBUTES = [ - QgsField("inhomogeneity_id", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("aquifer_top", QVariant.Double), - QgsField("aquifer_bottom", QVariant.Double), - QgsField("aquitard_c", QVariant.Double), - QgsField("aquifer_k", QVariant.Double), - QgsField("semiconf_top", QVariant.Double), - QgsField("semiconf_head", QVariant.Double), - QgsField("aquitard_s", QVariant.Double), - QgsField("aquifer_s", QVariant.Double), - QgsField("aquitard_npor", QVariant.Double), - QgsField("aquifer_npor", QVariant.Double), -] - - -class NameDialog(QDialog): - def __init__(self, parent=None): - super().__init__(parent) - self.name_line_edit = QLineEdit() - self.ok_button = QPushButton("OK") - self.cancel_button = QPushButton("Cancel") - self.ok_button.clicked.connect(self.accept) - self.cancel_button.clicked.connect(self.reject) - first_row = QHBoxLayout() - first_row.addWidget(QLabel("Layer name")) - first_row.addWidget(self.name_line_edit) - second_row = QHBoxLayout() - second_row.addStretch() - second_row.addWidget(self.ok_button) - second_row.addWidget(self.cancel_button) - layout = QVBoxLayout() - layout.addLayout(first_row) - layout.addLayout(second_row) - self.setLayout(layout) - - -class Element: - """ - Abstract base class for "ordinary" timml elements. - """ - - element_type = None - geometry_type = None - timml_attributes = () - ttim_attributes = () - assoc_attributes = () - transient_columns = () - timml_defaults = {} - ttim_defaults = {} - assoc_defaults = {} - - def _initialize_default(self, path, name): - self.name = name - self.path = path - self.timml_name = None - self.ttim_name = None - self.assoc_name = None - self.timml_layer = None - self.ttim_layer = None - self.assoc_layer = None - self.item = None - - def __init__(self, path: str, name: str): - self._initialize_default(path, name) - self.timml_name = f"timml {self.element_type}:{name}" - - @classmethod - def dialog(cls, path: str, crs: Any, iface: Any, names: List[str]): - dialog = NameDialog() - dialog.show() - ok = dialog.exec_() - if not ok: - return - - name = dialog.name_line_edit.text() - if name in names: - raise ValueError(f"Name already exists in geopackage: {name}") - - instance = cls(path, name) - instance.create_layers(crs) - return instance - - def create_layer( - self, crs: Any, geometry_type: str, name: str, attributes: List - ) -> QgsVectorLayer: - layer = QgsVectorLayer(geometry_type, name, "memory") - provider = layer.dataProvider() - provider.addAttributes(attributes) - layer.updateFields() - layer.setCrs(crs) - return layer - - def create_timml_layer(self, crs: Any): - self.timml_layer = self.create_layer( - crs=crs, - geometry_type=self.geometry_type, - name=self.timml_name, - attributes=self.timml_attributes, - ) - - def create_ttim_layer(self, crs: Any): - pass - - def create_assoc_layer(self, crs: Any): - pass - - def create_layers(self, crs: Any): - self.create_timml_layer(crs) - self.create_ttim_layer(crs) - self.create_assoc_layer(crs) - - def set_defaults(self): - for layer, defaults in zip( - (self.timml_layer, self.ttim_layer, self.assoc_layer), - (self.timml_defaults, self.ttim_defaults, self.assoc_defaults), - ): - if layer is None: - continue - fields = layer.fields() - for name, definition in defaults.items(): - index = fields.indexFromName(name) - layer.setDefaultValueDefinition(index, definition) - return - - @staticmethod - def marker_renderer(**kwargs): - symbol = QgsMarkerSymbol.createSimple(kwargs) - return QgsSingleSymbolRenderer(symbol) - - @staticmethod - def line_renderer(**kwargs): - symbol = QgsLineSymbol.createSimple(kwargs) - return QgsSingleSymbolRenderer(symbol) - - @staticmethod - def polygon_renderer(**kwargs): - symbol = QgsFillSymbol.createSimple(kwargs) - return QgsSingleSymbolRenderer(symbol) - - @property - def renderer(self): - return None - - def timml_layer_from_geopackage(self) -> QgsVectorLayer: - self.timml_layer = QgsVectorLayer( - f"{self.path}|layername={self.timml_name}", self.timml_name - ) - - def ttim_layer_from_geopackage(self): - return - - def assoc_layer_from_geopackage(self): - return - - def load_layers_from_geopackage(self) -> None: - self.timml_layer_from_geopackage() - self.ttim_layer_from_geopackage() - self.assoc_layer_from_geopackage() - self.set_defaults() - return - - def write(self): - self.timml_layer = geopackage.write_layer( - self.path, self.timml_layer, self.timml_name - ) - self.set_defaults() - - def remove_from_geopackage(self): - geopackage.remove_layer(self.path, self.timml_name) - - def on_transient_changed(self, transient: bool): - if len(self.transient_columns) == 0: - return - - config = self.timml_layer.attributeTableConfig() - columns = config.columns() - - for i, column in enumerate(columns): - if column.name in self.transient_columns: - config.setColumnHidden(i, not transient) - - self.timml_layer.setAttributeTableConfig(config) - return - - -class TransientElement(Element): - """ - Abstract base class for transient (ttim) elements. - """ - - def __init__(self, path: str, name: str): - self._initialize_default(path, name) - self.timml_name = f"timml {self.element_type}:{name}" - self.ttim_name = f"ttim {self.element_type}:{name}" - - def create_ttim_layer(self, crs: Any): - self.ttim_layer = self.create_layer( - crs=crs, - geometry_type="No Geometry", - name=self.ttim_name, - attributes=self.ttim_attributes, - ) - - def ttim_layer_from_geopackage(self): - self.ttim_layer = QgsVectorLayer( - f"{self.path}|layername={self.ttim_name}", - self.ttim_name, - ) - - def write(self): - self.timml_layer = geopackage.write_layer( - self.path, self.timml_layer, self.timml_name - ) - self.ttim_layer = geopackage.write_layer( - self.path, self.ttim_layer, self.ttim_name - ) - self.set_defaults() - - def remove_from_geopackage(self): - geopackage.remove_layer(self.path, self.timml_name) - geopackage.remove_layer(self.path, self.ttim_name) - - -class AssociatedElement(Element): - """ - Abstract class for elements that require associated tables such as - Inhomogenities. - """ - - def __init__(self, path: str, name: str): - self._initialize_default(path, name) - self.timml_name = f"timml {self.element_type}:{name}" - self.assoc_name = f"timml {self.element_type} Properties:{name}" - - def create_assoc_layer(self, crs: Any): - self.assoc_layer = self.create_layer( - crs=crs, - geometry_type="No Geometry", - name=self.assoc_name, - attributes=self.assoc_attributes, - ) - - def assoc_layer_from_geopackage(self): - self.assoc_layer = QgsVectorLayer( - f"{self.path}|layername={self.assoc_name}", - self.assoc_name, - ) - - def write(self): - self.timml_layer = geopackage.write_layer( - self.path, self.timml_layer, self.timml_name - ) - self.assoc_layer = geopackage.write_layer( - self.path, self.assoc_layer, self.assoc_name - ) - self.set_defaults() - - def remove_from_geopackage(self): - geopackage.remove_layer(self.path, self.timml_name) - geopackage.remove_layer(self.path, self.assoc_name) - - -class Domain(TransientElement): - element_type = "Domain" - geometry_type = "Polygon" - ttim_attributes = (QgsField("time", QVariant.Double),) - - def __init__(self, path: str, name: str): - self._initialize_default(path, name) - self.timml_name = f"timml {self.element_type}:Domain" - self.ttim_name = "ttim Computation Times:Domain" - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - """ - Results in transparent fill, with a medium thick black border line. - """ - return self.polygon_renderer( - color="255,0,0,0", color_border=BLACK, width_border="0.75" - ) - - def remove_from_geopackage(self): - pass - - def update_extent(self, iface: Any) -> Tuple[float, float]: - provider = self.timml_layer.dataProvider() - provider.truncate() # removes all features - canvas = iface.mapCanvas() - extent = canvas.extent() - xmin = extent.xMinimum() - ymin = extent.yMinimum() - xmax = extent.xMaximum() - ymax = extent.yMaximum() - points = [ - QgsPointXY(xmin, ymax), - QgsPointXY(xmax, ymax), - QgsPointXY(xmax, ymin), - QgsPointXY(xmin, ymin), - ] - feature = QgsFeature() - feature.setGeometry(QgsGeometry.fromPolygonXY([points])) - provider.addFeatures([feature]) - canvas.refresh() - return ymax, ymin - - -class Aquifer(TransientElement): - element_type = "Aquifer" - geometry_type = "No Geometry" - timml_attributes = AQUIFER_ATTRIBUTES.copy() - ttim_attributes = ( - QgsField("time_min", QVariant.Double), - QgsField("time_max", QVariant.Double), - QgsField("time_start", QVariant.Double), - QgsField("stehfest_M", QVariant.Int), - QgsField("reference_date", QVariant.DateTime), - ) - ttim_defaults = { - "time_min": QgsDefaultValue("0.01"), - "time_max": QgsDefaultValue("10.0"), - "time_start": QgsDefaultValue("0.0"), - "stehfest_M": QgsDefaultValue("10"), - } - transient_columns = ( - "aquitard_s", - "aquifer_s", - "aquitard_npor", - "aquifer_npor", - ) - - def __init__(self, path: str, name: str): - self._initialize_default(path, name) - self.timml_name = f"timml {self.element_type}:Aquifer" - self.ttim_name = "ttim Temporal Settings:Aquifer" - - def write(self): - self.timml_layer = geopackage.write_layer( - self.path, self.timml_layer, self.timml_name, newfile=True - ) - self.ttim_layer = geopackage.write_layer( - self.path, self.ttim_layer, self.ttim_name - ) - self.set_defaults() - - def remove_from_geopackage(self): - """This element may not be removed.""" - return - - -class UniformFlow(Element): - element_type = "Uniform Flow" - geometry_type = "No geometry" - timml_attributes = ( - QgsField("slope", QVariant.Double), - QgsField("angle", QVariant.Double), - QgsField("label", QVariant.String), - ) - - -class Constant(Element): - element_type = "Constant" - geometry_type = "Point" - timml_attributes = ( - QgsField("head", QVariant.Double), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.marker_renderer(color=RED, name="star", size="5") - - -class Observation(TransientElement): - element_type = "Observation" - geometry_type = "Point" - timml_attributes = ( - QgsField("label", QVariant.String), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time", QVariant.Double), - ) - timml_defaults = { - "timeseries_id": QgsDefaultValue("1"), - } - ttim_defaults = { - "timeseries_id": QgsDefaultValue("1"), - } - transient_columns = ("timeseries_id",) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.marker_renderer(color=LIGHT_BLUE, name="triangle", size="3") - - -class Well(TransientElement): - element_type = "Well" - geometry_type = "Point" - timml_attributes = ( - QgsField("discharge", QVariant.Double), - QgsField("radius", QVariant.Double), - QgsField("resistance", QVariant.Double), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - QgsField("time_start", QVariant.Double), - QgsField("time_end", QVariant.Double), - QgsField("discharge_transient", QVariant.Double), - QgsField("caisson_radius", QVariant.Double), - QgsField("slug", QVariant.Bool), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time_start", QVariant.Double), - QgsField("discharge", QVariant.Double), - ) - transient_columns = ( - "time_start", - "time_end", - "discharge_transient", - "caisson_radius", - "slug", - "timeseries_id", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.marker_renderer(color=GREEN, size="3") - - -class HeadWell(TransientElement): - element_type = "Head Well" - geometry_type = "Point" - timml_attributes = ( - QgsField("head", QVariant.Double), - QgsField("radius", QVariant.Double), - QgsField("resistance", QVariant.Double), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - QgsField("time_start", QVariant.Double), - QgsField("time_end", QVariant.Double), - QgsField("head_transient", QVariant.Double), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time_start", QVariant.Double), - QgsField("head", QVariant.Double), - ) - transient_columns = ( - "time_start", - "time_end", - "head_transient", - "timeseries_id", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.marker_renderer(color=BLUE, size="3") - - -class HeadLineSink(TransientElement): - element_type = "Head Line Sink" - geometry_type = "Linestring" - timml_attributes = ( - QgsField("head", QVariant.Double), - QgsField("resistance", QVariant.Double), - QgsField("width", QVariant.Double), - QgsField("order", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - QgsField("time_start", QVariant.Double), - QgsField("time_end", QVariant.Double), - QgsField("head_transient", QVariant.Double), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time_start", QVariant.Double), - QgsField("head", QVariant.Double), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - } - transient_columns = ( - "time_start", - "time_end", - "head_transient", - "timeseries_id", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.line_renderer(color=BLUE, width="0.75") - - -class LineSinkDitch(TransientElement): - element_type = "Line Sink Ditch" - geometry_type = "Linestring" - timml_attributes = ( - QgsField("discharge", QVariant.Double), - QgsField("resistance", QVariant.Double), - QgsField("width", QVariant.Double), - QgsField("order", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - QgsField("time_start", QVariant.Double), - QgsField("time_end", QVariant.Double), - QgsField("discharge_transient", QVariant.Double), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time_start", QVariant.Double), - QgsField("discharge", QVariant.Double), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - } - transient_columns = ( - "time_start", - "time_end", - "discharge_transient", - "timeseries_id", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.line_renderer(color=GREEN, width="0.75") - - -class ImpermeableLineDoublet(Element): - element_type = "Impermeable Line Doublet" - geometry_type = "Linestring" - timml_attributes = ( - QgsField("order", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - } - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.line_renderer(color=RED, width="0.75") - - -class LeakyLineDoublet(Element): - element_type = "Leaky Line Doublet" - geometry_type = "Linestring" - timml_attributes = ( - QgsField("resistance", QVariant.Double), - QgsField("order", QVariant.Int), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - } - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.line_renderer(color=RED, width="0.75", outline_style="dash") - - -class CircularAreaSink(TransientElement): - element_type = "Circular Area Sink" - geometry_type = "Polygon" - timml_attributes = ( - QgsField("rate", QVariant.Double), - QgsField("layer", QVariant.Int), - QgsField("label", QVariant.String), - QgsField("time_start", QVariant.Double), - QgsField("time_end", QVariant.Double), - QgsField("rate_transient", QVariant.Double), - QgsField("timeseries_id", QVariant.Int), - ) - ttim_attributes = ( - QgsField("timeseries_id", QVariant.Int), - QgsField("time_start", QVariant.Double), - QgsField("rate", QVariant.Double), - ) - transient_columns = ( - "time_start", - "time_end", - "rate_transient", - "timeseries_id", - ) - - @property - def renderer(self): - return self.polygon_renderer( - color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75" - ) - - -class PolygonAreaSink(Element): - element_type = "Polygon Area Sink" - geometry_type = "Polygon" - timml_attributes = ( - QgsField("rate", QVariant.Double), - QgsField("order", QVariant.Int), - QgsField("ndegrees", QVariant.Int), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - "ndegrees": QgsDefaultValue("6"), - } - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.polygon_renderer( - color=TRANSPARENT_GREEN, color_border=GREEN, width_border="0.75" - ) - - -class PolygonSemiConfinedTop(Element): - element_type = "Polygon Semi-Confined Top" - geometry_type = "Polygon" - timml_attributes = ( - QgsField("aquitard_c", QVariant.Double), - QgsField("semiconf_top", QVariant.Double), - QgsField("semiconf_head", QVariant.Double), - QgsField("order", QVariant.Int), - QgsField("ndegrees", QVariant.Int), - ) - timml_defaults = { - "order": QgsDefaultValue("4"), - "ndegrees": QgsDefaultValue("6"), - } - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.polygon_renderer( - color=TRANSPARENT_BLUE, color_border=BLUE, width_border="0.75" - ) - - -class PolygonInhomogeneity(AssociatedElement): - element_type = "Polygon Inhomogeneity" - geometry_type = "Polygon" - timml_attributes = ( - QgsField("inhomogeneity_id", QVariant.Int), - QgsField("order", QVariant.Int), - QgsField("ndegrees", QVariant.Int), - ) - assoc_attributes = INHOM_ATTRIBUTES.copy() - timml_defaults = { - "order": QgsDefaultValue("4"), - "ndegrees": QgsDefaultValue("6"), - "inhomogeneity_id": QgsDefaultValue("1"), - } - assoc_defaults = { - "inhomogeneity_id": QgsDefaultValue("1"), - } - transient_columns = ( - "aquitard_s", - "aquifer_s", - "aquitard_npor", - "aquifer_npor", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.polygon_renderer( - color=TRANSPARENT_GREY, color_border=GREY, width_border="0.75" - ) - - def on_transient_changed(self, transient: bool): - config = self.assoc_layer.attributeTableConfig() - columns = config.columns() - - for i, column in enumerate(columns): - if column.name in self.transient_columns: - config.setColumnHidden(i, not transient) - - self.assoc_layer.setAttributeTableConfig(config) - return - - -class BuildingPit(AssociatedElement): - element_type = "Building Pit" - geometry_type = "Polygon" - timml_attributes = ( - QgsField("inhomogeneity_id", QVariant.Int), - QgsField("order", QVariant.Int), - QgsField("ndegrees", QVariant.Int), - QgsField("layer", QVariant.Int), - ) - assoc_attributes = BUILDING_PIT_ATTRIBUTES.copy() - timml_defaults = { - "inhomogeneity_id": QgsDefaultValue("1"), - "order": QgsDefaultValue("4"), - "ndegrees": QgsDefaultValue("6"), - } - assoc_defaults = { - "inhomogeneity_id": QgsDefaultValue("1"), - } - transient_columns = ( - "aquitard_s", - "aquifer_s", - "aquitard_npor", - "aquifer_npor", - ) - - @property - def renderer(self) -> QgsSingleSymbolRenderer: - return self.polygon_renderer( - color=TRANSPARENT_RED, color_border=RED, width_border="0.75" - ) - - def on_transient_changed(self, transient: bool): - config = self.assoc_layer.attributeTableConfig() - columns = config.columns() - - for i, column in enumerate(columns): - if column.name in self.transient_columns: - config.setColumnHidden(i, not transient) - - self.assoc_layer.setAttributeTableConfig(config) - return - - -ELEMENTS = { - element.element_type: element - for element in ( - Aquifer, - Domain, - Constant, - UniformFlow, - Well, - HeadWell, - HeadLineSink, - LineSinkDitch, - CircularAreaSink, - ImpermeableLineDoublet, - LeakyLineDoublet, - PolygonAreaSink, - PolygonSemiConfinedTop, - PolygonInhomogeneity, - BuildingPit, - Observation, - ) -} - - -def parse_name(layername: str) -> Tuple[str, str, str]: - """ - Based on the layer name find out: - - * whether it's a timml or ttim element; - * which element type it is; - * what the user provided name is. - - For example: - parse_name("timml Headwell: drainage") -> ("timml", "Head Well", "drainage") - - This function can also be found in gistim.common - """ - prefix, name = layername.split(":") - element_type = re.split("timml |ttim ", prefix)[1] - mapping = { - "Computation Times": "Domain", - "Temporal Settings": "Aquifer", - "Polygon Inhomogeneity Properties": "Polygon Inhomogeneity", - "Building Pit Properties": "Building Pit", - } - element_type = mapping.get(element_type, element_type) - if "timml" in prefix: - if "Properties" in prefix: - tim_type = "timml_assoc" - else: - tim_type = "timml" - elif "ttim" in prefix: - tim_type = "ttim" - else: - raise ValueError("Neither timml nor ttim in layername") - return tim_type, element_type, name - - -def load_elements_from_geopackage(path: str) -> List[Element]: - # List the names in the geopackage - gpkg_names = geopackage.layers(path) - - # Group them on the basis of name - dd = defaultdict - grouped_names = dd(partial(dd, partial(dd, list))) - for layername in gpkg_names: - tim_type, element_type, name = parse_name(layername) - grouped_names[element_type][name][tim_type] = layername - - elements = [] - for element_type, group in grouped_names.items(): - for name in group: - elements.append(ELEMENTS[element_type](path, name)) - - return elements diff --git a/plugin/qgistim/widgets/compute_widget.py b/plugin/qgistim/widgets/compute_widget.py index 1e45a8f..9e187f4 100644 --- a/plugin/qgistim/widgets/compute_widget.py +++ b/plugin/qgistim/widgets/compute_widget.py @@ -6,7 +6,6 @@ from PyQt5.QtCore import Qt from PyQt5.QtWidgets import ( QCheckBox, - QComboBox, QDoubleSpinBox, QFileDialog, QGroupBox, @@ -20,20 +19,17 @@ from qgis.core import ( QgsApplication, QgsMapLayerProxyModel, - QgsMarkerSymbol, QgsMeshDatasetIndex, QgsMeshLayer, QgsProject, QgsRasterLayer, - QgsSingleSymbolRenderer, QgsTask, QgsVectorLayer, - QgsVectorLayerTemporalProperties, ) from qgis.gui import QgsMapLayerComboBox from qgistim.core import geopackage, layer_styling -from qgistim.core.dummy_ugrid import write_dummy_ugrid -from qgistim.core.processing import mesh_contours +from qgistim.core.elements import ELEMENTS, parse_name +from qgistim.core.processing import mesh_contours, set_temporal_properties from qgistim.core.task import BaseServerTask @@ -59,9 +55,14 @@ def finished(self, result): self.parent.set_interpreter_interaction(True) if result: self.push_success_message() - self.parent.load_mesh_result(self.data["outpath"]) - self.parent.load_raster_result(self.data["outpath"]) - self.parent.load_vector_result(self.data["outpath"]) + self.parent.clear_outdated_output(self.data["path"]) + if self.data["head_observations"] or self.data["discharge"]: + self.parent.load_vector_result(self.data["path"]) + if self.data["mesh"]: + self.parent.load_mesh_result(self.data["path"], self.data["contours"]) + if self.data["raster"]: + self.parent.load_raster_result(self.data["path"]) + else: self.push_failure_message() return @@ -74,39 +75,53 @@ def __init__(self, parent=None): self.start_task = None self.parent = parent - self.dummy_ugrid_path = Path(tempfile.mkdtemp()) / "qgistim-dummy-ugrid.nc" - write_dummy_ugrid(self.dummy_ugrid_path) - self.domain_button = QPushButton("Set to current extent") - self.transient_combo_box = QComboBox() - self.transient_combo_box.addItems(["Steady-state", "Transient"]) - self.transient_combo_box.currentTextChanged.connect(self.on_transient_changed) self.compute_button = QPushButton("Compute") self.compute_button.clicked.connect(self.compute) + + self.mesh_checkbox = QCheckBox("Mesh") + self.raster_checkbox = QCheckBox("Raster") + self.contours_checkbox = QCheckBox("Contours") + self.head_observations_checkbox = QCheckBox("Head Observations") + self.discharge_checkbox = QCheckBox("Discharge") + self.discharge_observations_checkbox = QCheckBox("Discharge Observations") + self.cellsize_spin_box = QDoubleSpinBox() self.cellsize_spin_box.setMinimum(0.0) self.cellsize_spin_box.setMaximum(10_000.0) self.cellsize_spin_box.setSingleStep(1.0) self.cellsize_spin_box.setValue(25.0) self.domain_button.clicked.connect(self.domain) + # By default: all output + self.mesh_checkbox.setChecked(True) + self.raster_checkbox.setChecked(False) + self.contours_checkbox.setChecked(False) + self.head_observations_checkbox.setChecked(True) + self.discharge_checkbox.setChecked(False) + + self.mesh_checkbox.toggled.connect(self.contours_checkbox.setEnabled) + self.mesh_checkbox.toggled.connect( + lambda checked: not checked and self.contours_checkbox.setChecked(False) + ) + # self.mesh_checkbox = QCheckBox("Trimesh") self.output_line_edit = QLineEdit() self.output_button = QPushButton("Set path as ...") self.output_button.clicked.connect(self.set_output_path) - self.contour_checkbox = QCheckBox("Auto-generate contours") - self.contour_button = QPushButton("Export contours") - self.contour_button.clicked.connect(self.export_contours) + self.contour_button = QPushButton("Redraw contours") + self.contour_button.clicked.connect(self.redraw_contours) self.contour_layer = QgsMapLayerComboBox() self.contour_layer.setFilters(QgsMapLayerProxyModel.MeshLayer) self.contour_min_box = QDoubleSpinBox() self.contour_max_box = QDoubleSpinBox() self.contour_step_box = QDoubleSpinBox() + self.contour_max_box.setMaximum(1000.0) + self.contour_max_box.setValue(5.0) + # Ensure the maximum cannot dip below the min box value. + self.contour_min_box.valueChanged.connect(self.set_minimum_contour_stop) self.contour_min_box.setMinimum(-1000.0) self.contour_min_box.setMaximum(1000.0) self.contour_min_box.setValue(-5.0) - self.contour_max_box.setMinimum(-1000.0) - self.contour_max_box.setMaximum(1000.0) - self.contour_max_box.setValue(5.0) self.contour_step_box.setSingleStep(0.1) self.contour_step_box.setValue(0.5) @@ -128,14 +143,21 @@ def __init__(self, parent=None): domain_layout.addWidget(self.domain_button) domain_layout.addLayout(domain_row) - result_row1 = QHBoxLayout() - result_row1.addWidget(self.transient_combo_box) - result_row1.addWidget(self.compute_button) - result_row2 = QHBoxLayout() - result_row2.addWidget(self.output_line_edit) - result_row2.addWidget(self.output_button) - result_layout.addLayout(result_row1) - result_layout.addLayout(result_row2) + output_row = QHBoxLayout() + output_row.addWidget(self.output_line_edit) + output_row.addWidget(self.output_button) + + button_row = QHBoxLayout() + button_row.addWidget(self.compute_button) + result_layout.addLayout(output_row) + + result_layout.addWidget(self.mesh_checkbox) + result_layout.addWidget(self.raster_checkbox) + result_layout.addWidget(self.contours_checkbox) + result_layout.addWidget(self.head_observations_checkbox) + result_layout.addWidget(self.discharge_checkbox) + + result_layout.addLayout(button_row) contour_row1 = QHBoxLayout() to_label = QLabel("to") @@ -150,7 +172,6 @@ def __init__(self, parent=None): contour_row2 = QHBoxLayout() contour_row2.addWidget(self.contour_layer) contour_row2.addWidget(self.contour_button) - contour_layout.addWidget(self.contour_checkbox) contour_layout.addLayout(contour_row1) contour_layout.addLayout(contour_row2) @@ -162,27 +183,27 @@ def __init__(self, parent=None): def reset(self): self.cellsize_spin_box.setValue(25.0) - self.contour_checkbox.setCheckState(False) self.transient_combo_box.setCurrentIndex(0) self.output_line_edit.setText("") + self.raster_checkbox.setCheckState(False) + self.mesh_checkbox.setCheckState(True) + self.contours_checkbox.setCheckState(False) + self.head_observations_checkbox.setCheckState(True) + self.discharge_checkbox.setCheckState(False) self.contour_min_box.setValue(-5.0) self.contour_max_box.setValue(5.0) self.contour_step_box.setValue(0.5) return + def set_minimum_contour_stop(self) -> None: + self.contour_max_box.setMinimum(self.contour_min_box.value() + 0.05) + def set_interpreter_interaction(self, value: bool): self.parent.set_interpreter_interaction(value) def shutdown_server(self): self.parent.shutdown_server() - @property - def transient(self) -> bool: - return self.transient_combo_box.currentText() == "Transient" - - def on_transient_changed(self) -> None: - self.parent.on_transient_changed() - def contour_range(self) -> Tuple[float, float, float]: return ( float(self.contour_min_box.value()), @@ -198,32 +219,78 @@ def add_contour_layer(self, layer) -> None: self.parent.output_group.add_layer( layer, "vector", renderer=renderer, on_top=True, labels=labels ) + return + + @property + def output_path(self) -> str: + return self.output_line_edit.text() - def export_contours(self) -> None: + def clear_outdated_output(self, path: str) -> None: + path = Path(path) + gpkg_path = path.with_suffix(".output.gpkg") + netcdf_paths = (path.with_suffix(".nc"), path.with_suffix(".ugrid.nc")) + for layer in QgsProject.instance().mapLayers().values(): + source = layer.source() + if ( + Path(source) in netcdf_paths + or Path(source.partition("|")[0]) == gpkg_path + ): + QgsProject.instance().removeMapLayer(layer.id()) + return + + def redraw_contours(self) -> None: + path = Path(self.output_path) layer = self.contour_layer.currentLayer() + if layer is None: + return + + start, stop, step = self.contour_range() + if (start == stop) or (step == 0.0): + return + renderer = layer.rendererSettings() index = renderer.activeScalarDatasetGroup() qgs_index = QgsMeshDatasetIndex(group=index, dataset=0) name = layer.datasetGroupMetadata(qgs_index).name() - start, stop, step = self.contour_range() + contours_name = f"{path.stem}-contours-{name}" + gpkg_path = str(path.with_suffix(".output.gpkg")) + layer = mesh_contours( + gpkg_path=gpkg_path, layer=layer, index=index, - name=name, + name=contours_name, start=start, stop=stop, step=step, ) - self.add_contour_layer(layer) + + # Re-use layer if it already exists. Otherwise add a new layer. + project_layers = { + layer.name(): layer for layer in QgsProject.instance().mapLayers().values() + } + project_layer = project_layers.get(contours_name) + if ( + (project_layer is not None) + and (project_layer.name() == layer.name()) + and (project_layer.source() == layer.source()) + ): + project_layer.reload() + else: + self.add_contour_layer(layer) return def set_output_path(self) -> None: - current = self.output_line_edit.text() - path, _ = QFileDialog.getSaveFileName(self, "Save output as...", current, "*") + current = self.output_path + path, _ = QFileDialog.getSaveFileName( + self, "Save output as...", current, "*.gpkg" + ) + if path != "": # Empty string in case of cancel button press - self.output_line_edit.setText(path) + self.output_line_edit.setText(str(Path(path).with_suffix(""))) # Note: Qt does pretty good validity checking of the Path in the # Dialog, there is no real need to validate path here. + return def set_default_path(self, text: str) -> None: """ @@ -233,28 +300,37 @@ def set_default_path(self, text: str) -> None: return path = Path(text) self.output_line_edit.setText(str(path.parent / path.stem)) + return def compute(self) -> None: """ Run a TimML computation with the current state of the currently active GeoPackage dataset. """ - active_elements = self.parent.active_elements() cellsize = self.cellsize_spin_box.value() - inpath = Path(self.parent.path).absolute() - outpath = Path(self.output_line_edit.text()).absolute() - mode = self.transient_combo_box.currentText().lower() + transient = self.parent.dataset_widget.transient + output_options = { + "raster": self.raster_checkbox.isChecked(), + "mesh": self.mesh_checkbox.isChecked(), + "contours": self.contours_checkbox.isChecked(), + "head_observations": self.head_observations_checkbox.isChecked(), + "discharge": self.discharge_checkbox.isChecked(), + } + + path = Path(self.output_path).absolute().with_suffix(".json") + invalid_input = self.parent.dataset_widget.convert_to_json( + path, cellsize=cellsize, transient=transient, output_options=output_options + ) + # Early return in case some problems are found. + if invalid_input: + return + data = { "operation": "compute", - "inpath": str(inpath), - "outpath": str(outpath), - "cellsize": cellsize, - "mode": mode, - "active_elements": active_elements, + "path": str(path), + "transient": transient, + **output_options, } - # import json - # print(json.dumps(data)) - # # https://gis.stackexchange.com/questions/296175/issues-with-qgstask-and-task-manager # It seems the task goes awry when not associated with a Python object! # -- we just assign it to the widget here. @@ -264,7 +340,7 @@ def compute(self) -> None: # task.finished(result) # Remove the output layers from QGIS, otherwise they cannot be overwritten. - gpkg_path = str(outpath) + gpkg_path = str(path) for layer in QgsProject.instance().mapLayers().values(): if Path(gpkg_path) == Path(layer.source()): QgsProject.instance().removeMapLayer(layer.id()) @@ -277,6 +353,7 @@ def compute(self) -> None: ) self.set_interpreter_interaction(False) QgsApplication.taskManager().addTask(self.compute_task) + return def domain(self) -> None: """ @@ -286,6 +363,7 @@ def domain(self) -> None: ymax, ymin = item.element.update_extent(self.parent.iface) self.set_cellsize_from_domain(ymax, ymin) self.parent.iface.mapCanvas().refreshAllLayers() + return def set_cellsize_from_domain(self, ymax: float, ymin: float) -> None: # Guess a reasonable value for the cellsize: about 50 rows @@ -299,30 +377,15 @@ def set_cellsize_from_domain(self, ymax: float, ymin: float) -> None: elif dy > 1.0: dy = round(dy) self.cellsize_spin_box.setValue(dy) + return - def contouring(self) -> Tuple[bool, float, float, float]: - contour = self.contour_checkbox.isChecked() - start, stop, step = self.contour_range() - return contour, start, stop, step - - def load_mesh_result(self, path: Union[Path, str]) -> None: + def load_mesh_result(self, path: Union[Path, str], load_contours: bool) -> None: path = Path(path) # String for QGIS functions netcdf_path = str(path.with_suffix(".ugrid.nc")) - # Loop through layers first. If the path already exists as a layer source, remove it. - # Otherwise QGIS will not the load the new result (this feels like a bug?). - for layer in QgsProject.instance().mapLayers().values(): - if Path(netcdf_path) == Path(layer.source()): - QgsProject.instance().removeMapLayer(layer.id()) - # Ensure the file is properly released by loading a dummy - QgsMeshLayer(str(self.dummy_ugrid_path), "", "mdal") - layer = QgsMeshLayer(netcdf_path, f"{path.stem}", "mdal") indexes = layer.datasetGroupsIndexes() - - contour, start, stop, step = self.contouring() contour_layers = [] - for index in indexes: qgs_index = QgsMeshDatasetIndex(group=index, dataset=0) name = layer.datasetGroupMetadata(qgs_index).name() @@ -339,11 +402,20 @@ def load_mesh_result(self, path: Union[Path, str]) -> None: index_layer.setRendererSettings(renderer) self.parent.output_group.add_layer(index_layer, "mesh") - if contour: + if load_contours: + # Should generally result in 20 contours. + start = scalar_settings.classificationMinimum() + stop = scalar_settings.classificationMaximum() + step = (stop - start) / 21 + # If no head differences are present, no contours can be drawn. + if step == 0.0: + return + contour_layer = mesh_contours( + gpkg_path=str(path.with_suffix(".output.gpkg")), layer=index_layer, index=index, - name=name, + name=f"{path.stem}-contours-{name}", start=start, stop=stop, step=step, @@ -351,7 +423,7 @@ def load_mesh_result(self, path: Union[Path, str]) -> None: contour_layers.append(contour_layer) # Add the contours in the appropriate order: highest (deepest) layer first! - if contour: + if load_contours: for contour_layer in contour_layers[::-1]: self.add_contour_layer(contour_layer) @@ -368,12 +440,6 @@ def steady_or_first(name: str) -> bool: # String for QGIS functions path = Path(path) raster_path = str(path.with_suffix(".nc")) - for layer in QgsProject.instance().mapLayers().values(): - if Path(raster_path) == Path(layer.source()): - QgsProject.instance().removeMapLayer(layer.id()) - - # contour, start, stop, step = self.contouring() - layer = QgsRasterLayer(raster_path, "", "gdal") bands = [i + 1 for i in range(layer.bandCount())] @@ -388,71 +454,39 @@ def steady_or_first(name: str) -> bool: layer.setRenderer(renderer) self.parent.output_group.add_layer(layer, "raster") - # if contour: - # contour_layer = raster_steady_contours( - # layer=layer, - # name=name, - # start=start, - # stop=stop, - # step=step, - # ) - # self.add_contour_layer(contour_layer) - return def load_vector_result(self, path: Union[Path, str]) -> None: path = Path(path) - project_layers = { - layer.name(): layer for layer in QgsProject.instance().mapLayers().values() - } gpkg_path = path.with_suffix(".output.gpkg") if not gpkg_path.exists(): return for layername in geopackage.layers(str(gpkg_path)): - add = False layers_panel_name = f"{path.stem}-{layername}" - project_layer = project_layers.get(layers_panel_name) + + layer = QgsVectorLayer( + f"{gpkg_path}|layername={layername}", layers_panel_name + ) + # Set the temporal properties if it's a temporal layer + set_temporal_properties(layer) + + # Special-case the labelling for observations and discharge. if ( - project_layer is not None - and Path(project_layer.source().partition("|")[0]) == gpkg_path + "timml Head Observation:" in layername + or "ttim Head Observation" in layername ): - # Shares name and source. Just reload the layer. - layer = project_layer - layer.reload() + labels = layer_styling.number_labels("head_layer0") + elif "discharge-" in layername: + labels = layer_styling.number_labels("discharge_layer0") else: - layer = QgsVectorLayer( - f"{gpkg_path}|layername={layername}", layers_panel_name - ) - add = True + labels = None - # Set the temporal properties if it's a temporal layer - temporal_properties = layer.temporalProperties() - fields = [field.name() for field in layer.fields()] - if ("datetime_start" in fields) and ("datetime_end" in fields): - temporal_properties.setStartField("datetime_start") - temporal_properties.setEndField("datetime_end") - temporal_properties.setMode( - QgsVectorLayerTemporalProperties.ModeFeatureDateTimeStartAndEndFromFields - ) - temporal_properties.setIsActive(True) - else: - temporal_properties.setIsActive(False) - - if add: - if "timml Observation:" in layername or "ttim Observation" in layername: - labels = layer_styling.number_labels("head_layer0") - light_blue = "166,206,227,255" - symbol = QgsMarkerSymbol.createSimple( - dict(color=light_blue, name="triangle", size="3") - ) - renderer = QgsSingleSymbolRenderer(symbol) - else: - labels = None - renderer = None - self.parent.output_group.add_layer( - layer, "vector", renderer=renderer, labels=labels - ) + _, element_type, _ = parse_name(layername) + renderer = ELEMENTS[element_type].renderer() + self.parent.output_group.add_layer( + layer, "vector", renderer=renderer, labels=labels + ) return diff --git a/plugin/qgistim/widgets/dataset_widget.py b/plugin/qgistim/widgets/dataset_widget.py index 9e7457f..3eeb7bb 100644 --- a/plugin/qgistim/widgets/dataset_widget.py +++ b/plugin/qgistim/widgets/dataset_widget.py @@ -10,14 +10,16 @@ user chooses the transient simulation mode, a number of elements must be disabled (such as inhomogeneities). """ +import json from pathlib import Path from shutil import copy -from typing import List, Set +from typing import Any, Dict, List, NamedTuple, Set, Tuple from PyQt5.QtCore import Qt from PyQt5.QtWidgets import ( QAbstractItemView, QCheckBox, + QComboBox, QFileDialog, QHBoxLayout, QHeaderView, @@ -30,9 +32,10 @@ QVBoxLayout, QWidget, ) -from qgis.core import QgsApplication, QgsProject, QgsTask -from qgistim.core.task import BaseServerTask -from qgistim.core.tim_elements import Aquifer, Domain, load_elements_from_geopackage +from qgis.core import Qgis, QgsProject +from qgistim.core.elements import Aquifer, Domain, load_elements_from_geopackage +from qgistim.core.formatting import data_to_json, data_to_script +from qgistim.widgets.error_window import ValidationDialog SUPPORTED_TTIM_ELEMENTS = set( [ @@ -47,15 +50,15 @@ "Circular Area Sink", "Impermeable Line Doublet", "Leaky Line Doublet", - "Observation", + "Head Observation", ] ) -class ConvertTask(BaseServerTask): - @property - def task_description(self): - return "converting GeoPackage to Python script" +class Extraction(NamedTuple): + timml: Dict[str, Any] = None + ttim: Dict[str, Any] = None + success: bool = True class DatasetTreeWidget(QTreeWidget): @@ -196,6 +199,66 @@ def remove_geopackage_layers(self) -> None: return + def extract_data(self, transient: bool) -> Tuple[Dict[str, Any], Dict[str, Any]]: + """ + Extract the data of the Geopackage. + + Validates all data while converting, and returns a list of validation + errors if something is amiss. + """ + data = {} + errors = {} + elements = { + item.text(1): item.element + for item in self.items() + if item.timml_checkbox.isChecked() + } + + # First convert the aquifer, since we need its data to validate + # other elements. + name = "timml Aquifer:Aquifer" + aquifer = elements.pop(name) + aquifer_extraction = aquifer.extract_data(transient) + if aquifer_extraction.errors: + errors[name] = aquifer_extraction.errors + return errors, None + + raw_data = aquifer_extraction.data + aquifer_data = aquifer.aquifer_data(raw_data, transient=transient) + data[name] = aquifer_data + if transient: + data["reference_date"] = str(raw_data["reference_date"].toPyDateTime()) + + times = set() + other = {"aquifer layers": raw_data["layer"], "global_aquifer": raw_data} + for name, element in elements.items(): + print(name) + try: + extraction = element.extract_data(transient, other) + if extraction.errors: + errors[name] = extraction.errors + elif extraction.data: # skip empty tables + data[name] = extraction.data + if extraction.times: + times.update(extraction.times) + except RuntimeError as e: + if ( + e.args[0] + == "wrapped C/C++ object of type QgsVectorLayer has been deleted" + ): + # Delay of Qt garbage collection to blame? + pass + else: + raise e + + if transient: + if times: + data["timml Aquifer:Aquifer"]["tmax"] = max(times) + else: + errors["Model"] = {"TTim input:": ["No transient forcing defined."]} + + return errors, data + class DatasetWidget(QWidget): def __init__(self, parent): @@ -203,14 +266,16 @@ def __init__(self, parent): self.parent = parent self.dataset_tree = DatasetTreeWidget() self.start_task = None - self.convert_task = None self.dataset_tree.setSizePolicy(QSizePolicy.Preferred, QSizePolicy.Expanding) self.dataset_line_edit = QLineEdit() self.dataset_line_edit.setEnabled(False) # Just used as a viewing port self.new_geopackage_button = QPushButton("New") self.open_geopackage_button = QPushButton("Open") self.copy_geopackage_button = QPushButton("Copy") - self.remove_button = QPushButton("Remove from Dataset") + self.transient_combo_box = QComboBox() + self.transient_combo_box.addItems(["Steady-state", "Transient"]) + self.transient_combo_box.currentTextChanged.connect(self.on_transient_changed) + self.remove_button = QPushButton("Remove from GeoPackage") self.add_button = QPushButton("Add to QGIS") self.new_geopackage_button.clicked.connect(self.new_geopackage) self.open_geopackage_button.clicked.connect(self.open_geopackage) @@ -219,10 +284,11 @@ def __init__(self, parent): self.suppress_popup_checkbox.stateChanged.connect(self.suppress_popup_changed) self.remove_button.clicked.connect(self.remove_geopackage_layer) self.add_button.clicked.connect(self.add_selection_to_qgis) - self.convert_button = QPushButton("Convert GeoPackage to Python script") - self.convert_button.clicked.connect(self.convert) + self.convert_button = QPushButton("Export to Python script") + self.convert_button.clicked.connect(self.convert_to_python) # Layout dataset_layout = QVBoxLayout() + mode_row = QHBoxLayout() dataset_row = QHBoxLayout() layer_row = QHBoxLayout() dataset_row.addWidget(self.dataset_line_edit) @@ -230,6 +296,8 @@ def __init__(self, parent): dataset_row.addWidget(self.new_geopackage_button) dataset_row.addWidget(self.copy_geopackage_button) dataset_layout.addLayout(dataset_row) + mode_row.addWidget(self.transient_combo_box) + dataset_layout.addLayout(mode_row) dataset_layout.addWidget(self.dataset_tree) dataset_layout.addWidget(self.suppress_popup_checkbox) layer_row.addWidget(self.add_button) @@ -237,6 +305,7 @@ def __init__(self, parent): dataset_layout.addLayout(layer_row) dataset_layout.addWidget(self.convert_button) self.setLayout(dataset_layout) + self.validation_dialog = None @property def path(self) -> str: @@ -246,7 +315,6 @@ def path(self) -> str: def reset(self): # Set state back to defaults self.start_task = None - self.convert_task = None self.dataset_line_edit.setText("") self.dataset_tree.reset() return @@ -258,7 +326,7 @@ def add_item_to_qgis(self, item) -> None: suppress = self.suppress_popup_checkbox.isChecked() # Start adding the layers maplayer = self.parent.input_group.add_layer( - element.timml_layer, "timml", element.renderer, suppress + element.timml_layer, "timml", element.renderer(), suppress ) self.parent.input_group.add_layer(element.ttim_layer, "ttim") self.parent.input_group.add_layer(element.assoc_layer, "timml") @@ -293,15 +361,14 @@ def load_geopackage(self) -> None: for element in elements: self.dataset_tree.add_element(element) + transient = self.transient for item in self.dataset_tree.items(): self.add_item_to_qgis(item) - item.element.on_transient_changed(self.parent.transient) + item.element.on_transient_changed(transient) self.dataset_tree.sortByColumn(0, Qt.SortOrder.AscendingOrder) self.parent.toggle_element_buttons(True) - self.parent.on_transient_changed() - return - + self.on_transient_changed() return def new_geopackage(self) -> None: @@ -367,6 +434,14 @@ def remove_geopackage_layer(self) -> None: self.dataset_tree.remove_geopackage_layers() return + @property + def transient(self) -> bool: + return self.transient_combo_box.currentText() == "Transient" + + def on_transient_changed(self) -> None: + self.dataset_tree.on_transient_changed(self.transient) + return + def suppress_popup_changed(self): suppress = self.suppress_popup_checkbox.isChecked() for item in self.dataset_tree.items(): @@ -401,10 +476,6 @@ def selection_names(self) -> Set[str]: selection.append(item.assoc_item) return set([item.element.name for item in selection]) - def on_transient_changed(self, transient: bool) -> None: - self.dataset_tree.on_transient_changed(transient) - return - def add_element(self, element) -> None: self.dataset_tree.add_element(element) return @@ -413,23 +484,94 @@ def set_interpreter_interaction(self, value: bool): self.parent.set_interpreter_interaction(value) return - def convert(self) -> None: + def _extract_data(self, transient: bool) -> Extraction: + if self.validation_dialog: + self.validation_dialog.close() + self.validation_dialog = None + + errors, timml_data = self.dataset_tree.extract_data(transient=False) + if errors: + self.validation_dialog = ValidationDialog(errors) + return Extraction(success=False) + + ttim_data = None + if transient: + errors, ttim_data = self.dataset_tree.extract_data(transient=True) + if errors: + self.validation_dialog = ValidationDialog(errors) + return Extraction(success=False) + + return Extraction(timml=timml_data, ttim=ttim_data) + + def convert_to_python(self) -> None: + transient = self.transient outpath, _ = QFileDialog.getSaveFileName(self, "Select file", "", "*.py") if outpath == "": # Empty string in case of cancel button press return - data = { - "operation": "convert", - "inpath": self.path, - "outpath": outpath, - } - self.convert_task = ConvertTask(self, data, self.parent.message_bar) - self.start_task = self.parent.start_interpreter_task() - if self.start_task is not None: - self.convert_task.addSubTask( - self.start_task, [], QgsTask.ParentDependsOnSubTask - ) - - self.parent.set_interpreter_interaction(False) - QgsApplication.taskManager().addTask(self.convert_task) + extraction = self._extract_data(transient=transient) + if not extraction.success: + return + + script = data_to_script(extraction.timml, extraction.ttim) + with open(outpath, "w") as f: + f.write(script) + + self.parent.message_bar.pushMessage( + title="Info", + text=f"Converted geopackage to Python script: {outpath}", + level=Qgis.Info, + ) return + + def convert_to_json( + self, + path: str, + cellsize: float, + transient: bool, + output_options: Dict[str, bool], + ) -> bool: + """ + Parameters + ---------- + path: str + Path to JSON file to write. + cellsize: float + Cell size to use to compute the head grid. + transient: bool + Steady-state (False) or transient (True). + + Returns + ------- + invalid_input: bool + Whether validation has succeeded. + """ + extraction = self._extract_data(transient=transient) + if not extraction.success: + return True + + json_data = data_to_json( + extraction.timml, + extraction.ttim, + cellsize=cellsize, + output_options=output_options, + ) + + crs = self.parent.crs + organization, srs_id = crs.authid().split(":") + json_data["crs"] = { + "description": crs.description(), + "organization": organization, + "srs_id": srs_id, + "wkt": crs.toWkt(), + } + + with open(path, "w") as fp: + json.dump(json_data, fp=fp, indent=4) + + self.parent.message_bar.pushMessage( + title="Info", + text=f"Converted geopackage to JSON: {path}", + level=Qgis.Info, + ) + return False diff --git a/plugin/qgistim/widgets/elements_widget.py b/plugin/qgistim/widgets/elements_widget.py index 9ac6454..c60c73a 100644 --- a/plugin/qgistim/widgets/elements_widget.py +++ b/plugin/qgistim/widgets/elements_widget.py @@ -1,7 +1,8 @@ from functools import partial from PyQt5.QtWidgets import QGridLayout, QPushButton, QVBoxLayout, QWidget -from qgistim.core.tim_elements import ELEMENTS +from qgis.core import Qgis +from qgistim.core.elements import ELEMENTS class ElementsWidget(QWidget): @@ -60,7 +61,13 @@ def tim_element(self, element_type: str) -> None: except ValueError: return - element = klass.dialog(self.parent.path, crs, self.parent.iface, names) + try: + element = klass.dialog(self.parent.path, crs, self.parent.iface, names) + except ValueError as e: + msg = str(e) + self.parent.message_bar.pushMessage("Error", msg, level=Qgis.Critical) + return + if element is None: # dialog cancelled return # Write to geopackage diff --git a/plugin/qgistim/widgets/error_window.py b/plugin/qgistim/widgets/error_window.py new file mode 100644 index 0000000..fc7be20 --- /dev/null +++ b/plugin/qgistim/widgets/error_window.py @@ -0,0 +1,63 @@ +from typing import Any, Dict + +from PyQt5.QtWidgets import QDialog, QHBoxLayout, QPushButton, QTextBrowser, QVBoxLayout + + +def format_list(errors: Dict[str, Any]): + """ + Format the a list of errors to HTML lists. + + Since the level of nesting may vary, this function is called recursively. + """ + messages = [] + for variable, var_errors in errors.items(): + if isinstance(var_errors, list): + messages.append(f"
{variable}
{variable}