"""
Model-only base classes for PyVertexvectors (no wx, no OpenGL dependency).
Author: HECE - University of Liege, Pierre Archambeau
Date: 2024-2026
Copyright (c) 2024 University of Liege. All rights reserved.
This script and its content are protected by copyright law. Unauthorized
copying or distribution of this file, via any medium, is strictly prohibited.
Provides
--------
- :class:`TriangulationModel` — triangulation geometry, I/O, transformations.
- :class:`vectorpropertiesModel` — visual property storage, I/O, color helpers.
GUI/OpenGL subclasses live in :mod:`wolfhece.PyVertexvectors`.
"""
from os import path
from typing import Union, Literal
import numpy as np
import struct
import io
import logging
import copy
import warnings
from pathlib import Path
from tqdm import tqdm
import triangle as tri_lib
import pyvista as pv
import pygltflib
from shapely.geometry import Point, Polygon, LineString, MultiLineString, MultiPoint, MultiPolygon
from shapely.ops import nearest_points
from shapely import prepare, is_prepared, destroy_prepared
from matplotlib.colors import Colormap
from matplotlib import cm
from scipy.interpolate import interp1d
from scipy.spatial import cKDTree
import geopandas as gpd
from concurrent.futures import ThreadPoolExecutor, wait
from ..PyTranslate import _
from ..PyVertex import wolfvertex, cloud_vertices
from ..color_constants import getRGBfromI, getIfromRGB, Colors, RGB
from ..PyPalette import wolfpalette
from .polygon_pbr_material import PolygonPBRMaterial
from enum import Enum
[docs]
class VectorOGLRenderer(Enum):
"""Rendering backend for vector objects."""
[docs]
LIST = 0 # Legacy display-list / immediate mode
[docs]
SHADER = 1 # Modern GLSL shader pipeline
# ====================================================================
# TriangulationModel
# ====================================================================
[docs]
class TriangulationModel:
"""Pure-data triangulation: vertices, triangles, I/O, transforms.
Holds arrays of points and triangle indices and provides geometry
operations, file reading/writing, and coordinate transformations
without any OpenGL or wx dependency.
The GUI subclass :class:`~wolfhece.PyVertexvectors.Triangulation`
adds ``Element_To_Draw`` integration, OpenGL display lists, and
matplotlib plotting helpers.
"""
def __init__(self, fn: str = '', pts=None, tri=None, idx: str = '') -> None:
[docs]
self.tri = tri if tri is not None else []
[docs]
self.pts = pts if pts is not None else []
[docs]
self._used_tri = [True] * self.nb_tri
[docs]
self._move_start = None
[docs]
self._rotation_center = None
[docs]
self._rotation_step = None
if fn != '':
self.filename = fn
self.read(fn)
else:
self.validate_format()
def __getstate__(self):
"""Get the state of the object for pickling."""
state = self.__dict__.copy()
state.pop('mapviewer', None)
return state
def __setstate__(self, state):
"""Set the state of the object for unpickling."""
self.__dict__.update(state)
# ----------------------------------------------------------------
# Properties
# ----------------------------------------------------------------
@property
[docs]
def nb_tri(self) -> int:
"""Return the number of triangles."""
return len(self.tri)
@property
[docs]
def nb_pts(self) -> int:
"""Return the number of points."""
return len(self.pts)
# ----------------------------------------------------------------
# Format helpers
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# PyVista conversion
# ----------------------------------------------------------------
[docs]
def as_polydata(self) -> pv.PolyData:
"""Convert the triangulation to a PyVista PolyData object."""
return pv.PolyData(
np.asarray(self.pts),
np.column_stack([[3] * self.nb_tri, self.tri]),
self.nb_tri,
)
[docs]
def from_polydata(self, poly: pv.PolyData):
"""Populate from a PyVista PolyData object."""
self.pts = np.asarray(poly.points.copy())
self.tri = np.asarray(
poly.faces.reshape([int(len(poly.faces) / 4), 4])[:, 1:4]
)
[docs]
def clip_surface(self, other: "TriangulationModel", invert=True, subdivide=0):
"""Clip the triangulation with another one."""
if subdivide == 0:
mypoly = self.as_polydata()
mycrop = other.as_polydata()
else:
mypoly = self.as_polydata().subdivide(subdivide)
mycrop = other.as_polydata().subdivide(subdivide)
res = mypoly.clip_surface(mycrop, invert=invert)
if len(res.faces) > 0:
self.from_polydata(res)
# ----------------------------------------------------------------
# Geometry analysis
# ----------------------------------------------------------------
[docs]
def get_mask(self, eps: float = 1e-10):
"""Return a boolean mask where triangle area ≤ *eps*."""
if self.nb_tri > 0:
v1 = [self.pts[curtri[1]][:2] - self.pts[curtri[0]][:2] for curtri in self.tri]
v2 = [self.pts[curtri[2]][:2] - self.pts[curtri[0]][:2] for curtri in self.tri]
self.areas = np.cross(v2, v1) / 2
return self.areas <= eps
[docs]
def _get_polygons(self) -> list[Polygon]:
"""Get a list of Shapely polygons from each triangle."""
polygons = []
for curtri in self.tri:
if len(curtri) == 3:
poly = Polygon([
self.pts[curtri[0]][:2],
self.pts[curtri[1]][:2],
self.pts[curtri[2]][:2],
])
if poly.is_valid:
polygons.append(poly)
else:
logging.debug('Invalid polygon found in triangulation: {}'.format(poly))
else:
logging.error('Triangle with {} vertices found in triangulation: {}'.format(len(curtri), curtri))
return polygons
[docs]
def unuse_triangles_containing_points(self, points: list[Point]):
"""Mark triangles containing any of *points* as unused."""
polys = self._get_polygons()
for point in points:
for i, poly in enumerate(polys):
if poly.contains(point):
self._used_tri[i] = False
[docs]
def get_triangles_as_listwolfvertices(self, used_only: bool = True) -> list[list[wolfvertex]]:
"""Return triangles as lists of :class:`wolfvertex`."""
triangles = []
for i, curtri in enumerate(self.tri):
if not used_only or (used_only and self._used_tri[i]):
if len(curtri) == 3:
v1 = wolfvertex(self.pts[curtri[0]][0], self.pts[curtri[0]][1], self.pts[curtri[0]][2])
v2 = wolfvertex(self.pts[curtri[1]][0], self.pts[curtri[1]][1], self.pts[curtri[1]][2])
v3 = wolfvertex(self.pts[curtri[2]][0], self.pts[curtri[2]][1], self.pts[curtri[2]][2])
triangles.append([v1, v2, v3])
else:
logging.warning('Triangle with {} vertices found in triangulation: {}'.format(len(curtri), curtri))
return triangles
# ----------------------------------------------------------------
# Bounds
# ----------------------------------------------------------------
[docs]
def find_minmax(self, force):
"""Compute bounding box from point coordinates.
:param force: if ``True``, always recompute.
"""
if force:
if self.nb_pts > 0:
self.xmin = float(np.min(self.pts[:, 0]))
self.ymin = float(np.min(self.pts[:, 1]))
self.xmax = float(np.max(self.pts[:, 0]))
self.ymax = float(np.max(self.pts[:, 1]))
# ----------------------------------------------------------------
# Transformations
# ----------------------------------------------------------------
[docs]
def set_cache(self):
"""Snapshot current points for cached transformations."""
self._cache = self.pts.copy()
[docs]
def clear_cache(self):
"""Discard cached point snapshot."""
self._cache = None
[docs]
def move(self, delta_x: float, delta_y: float, use_cache: bool = True):
"""Translate all points by *(delta_x, delta_y)*.
:param delta_x: delta x [m]
:param delta_y: delta y [m]
"""
if use_cache and self._cache is None:
self.set_cache()
if use_cache:
self.pts[:, 0] = self._cache[:, 0] + delta_x
self.pts[:, 1] = self._cache[:, 1] + delta_y
else:
self.pts[:, 0] += delta_x
self.pts[:, 1] += delta_y
[docs]
def rotate(self, angle: float, center: tuple | wolfvertex, use_cache: bool = True):
"""Rotate all points around *center*.
:param angle: angle in degrees — positive for clockwise rotation.
:param center: centre of rotation.
"""
if isinstance(center, wolfvertex):
center = (center.x, center.y)
angle = np.radians(angle)
c, s = np.cos(angle), np.sin(angle)
R = np.array([[c, -s], [s, c]])
if use_cache and self._cache is None:
self.set_cache()
if use_cache:
locxy = self._cache[:, :2] - np.array(center)
self.pts[:, :2] = np.dot(R, locxy.T).T + np.array(center)
else:
locxy = self.pts[:, :2] - np.array(center)
self.pts[:, :2] = np.dot(R, locxy.T).T + np.array(center)
[docs]
def rotate_xy(self, x: float, y: float, use_cache: bool = True):
"""Rotate around :attr:`_rotation_center` toward *(x, y)*."""
if self._rotation_center is None:
logging.error('No rotation center defined -- set it before rotating by this routine')
return self
angle = np.degrees(np.arctan2(-(y - self._rotation_center[1]), x - self._rotation_center[0]))
if self._rotation_step is not None:
angle = np.round(angle / self._rotation_step) * self._rotation_step
return self.rotate(-angle, center=self._rotation_center, use_cache=use_cache)
# ----------------------------------------------------------------
# Copy
# ----------------------------------------------------------------
[docs]
def copy(self) -> "TriangulationModel":
"""Return a deep copy of the triangulation."""
newtri = TriangulationModel(pts=self.pts.copy(), tri=self.tri.copy(), idx=self.idx + '_copy')
return newtri
# ----------------------------------------------------------------
# File I/O — GLTF
# ----------------------------------------------------------------
[docs]
def import_from_gltf(self, fn: str = ''):
"""Import a GLTF/GLB file.
:param fn: path to the ``.glb`` or ``.gltf`` file (required).
:raises ValueError: if *fn* is empty.
"""
if fn == '':
raise ValueError('A filename must be provided to import_from_gltf')
fn = str(fn)
gltf = pygltflib.GLTF2().load(fn)
mesh = gltf.meshes[gltf.scenes[gltf.scene].nodes[0]]
for primitive in mesh.primitives:
accessor = gltf.accessors[primitive.attributes.POSITION]
bufferView = gltf.bufferViews[accessor.bufferView]
buffer = gltf.buffers[bufferView.buffer]
data = gltf.get_data_from_buffer_uri(buffer.uri)
logging.info(_('Importing GLTF file : {}').format(fn))
logging.info(_('Number of vertices : {}').format(accessor.count))
points = np.zeros([accessor.count, 3], order='F', dtype=np.float64)
for i in range(accessor.count):
index = bufferView.byteOffset + accessor.byteOffset + i * 12
d = data[index:index + 12]
v = struct.unpack("<fff", d)
points[i, :] = np.asarray([v[0], -v[2], v[1]])
if np.mod(i, 100000) == 0:
logging.info(_('Reading vertex {}').format(i))
accessor = gltf.accessors[primitive.indices]
bufferView = gltf.bufferViews[accessor.bufferView]
buffer = gltf.buffers[bufferView.buffer]
data = gltf.get_data_from_buffer_uri(buffer.uri)
triangles = []
if accessor.componentType == pygltflib.UNSIGNED_SHORT:
size = 6
format = '<HHH'
elif accessor.componentType == pygltflib.UNSIGNED_INT:
size = 12
format = '<LLL'
logging.info(_('Number of triangles : {}').format(int(accessor.count / 3)))
for i in range(int(accessor.count / 3)):
index = bufferView.byteOffset + accessor.byteOffset + i * size
d = data[index:index + size]
v = struct.unpack(format, d)
triangles.append(list(v))
if np.mod(i, 100000) == 0:
logging.info(_('Reading triangle {}').format(i))
logging.info(_('Sorting information ...'))
xyz_u, indices = np.unique(np.array(points), return_inverse=True, axis=0)
logging.info(_('Creating triangles ...'))
triangles = [[indices[curtri[0]], indices[curtri[1]], indices[curtri[2]]] for curtri in list(triangles)]
self.pts = xyz_u
self.tri = triangles
[docs]
def export_to_gltf(self, fn: str = '') -> str:
"""Export the triangulation to a GLTF/GLB file.
:param fn: path to the output file (required).
:raises ValueError: if *fn* is empty.
:return: the path written to.
"""
if fn == '':
raise ValueError('A filename must be provided to export_to_gltf')
triangles = np.asarray(self.tri).astype(np.uint32)
points = np.column_stack([self.pts[:, 0], self.pts[:, 2], -self.pts[:, 1]]).astype(np.float32)
triangles_binary_blob = triangles.flatten().tobytes()
points_binary_blob = points.flatten().tobytes()
gltf = pygltflib.GLTF2(
scene=0,
scenes=[pygltflib.Scene(nodes=[0])],
nodes=[pygltflib.Node(mesh=0)],
meshes=[
pygltflib.Mesh(
primitives=[
pygltflib.Primitive(
attributes=pygltflib.Attributes(POSITION=1), indices=0
)
]
)
],
accessors=[
pygltflib.Accessor(
bufferView=0,
componentType=pygltflib.UNSIGNED_INT,
count=self.nb_tri * 3,
type=pygltflib.SCALAR,
max=[int(triangles.max())],
min=[int(triangles.min())],
),
pygltflib.Accessor(
bufferView=1,
componentType=pygltflib.FLOAT,
count=self.nb_pts,
type=pygltflib.VEC3,
max=points.max(axis=0).tolist(),
min=points.min(axis=0).tolist(),
),
],
bufferViews=[
pygltflib.BufferView(
buffer=0,
byteLength=len(triangles_binary_blob),
target=pygltflib.ELEMENT_ARRAY_BUFFER,
),
pygltflib.BufferView(
buffer=0,
byteOffset=len(triangles_binary_blob),
byteLength=len(points_binary_blob),
target=pygltflib.ARRAY_BUFFER,
),
],
buffers=[
pygltflib.Buffer(
byteLength=len(triangles_binary_blob) + len(points_binary_blob)
)
],
)
gltf.set_binary_blob(triangles_binary_blob + points_binary_blob)
gltf.save(fn)
return fn
# ----------------------------------------------------------------
# File I/O — binary .TRI
# ----------------------------------------------------------------
[docs]
def saveas(self, fn: str = ''):
"""Write the triangulation to a binary ``.TRI`` file.
Binary layout (little-endian):
- int32 number_of_points
- int32 64 or 32 (float64 / float32 flag)
- VEC3 points (float64 or float32)
- int32 number_of_triangles
- VEC3 triangle indices (uint32)
"""
if self.filename == '' and fn == '':
return
if fn != '':
self.filename = fn
triangles = np.asarray(self.tri).astype(np.uint32)
points = np.asarray(self.pts)
with open(self.filename, 'wb') as f:
f.write(self.nb_pts.to_bytes(4, 'little'))
if points.dtype == np.float64:
f.write(int(64).to_bytes(4, 'little'))
else:
f.write(int(32).to_bytes(4, 'little'))
points.tofile(f, "")
f.write(self.nb_tri.to_bytes(4, 'little'))
triangles.tofile(f, "")
[docs]
def read(self, fn: str):
"""Read a binary ``.TRI``, ``.dxf``, or ``.gltf/.glb`` file.
After reading, :meth:`validate_format` and :meth:`find_minmax`
are called automatically.
"""
fn = str(fn)
if fn.endswith('.dxf'):
self.import_dxf(fn)
elif fn.endswith('.gltf') or fn.endswith('.glb'):
self.import_from_gltf(fn)
elif fn.endswith('.tri'):
with open(fn, 'rb') as f:
nb_pts = int.from_bytes(f.read(4), 'little')
which_type = int.from_bytes(f.read(4), 'little')
if which_type == 64:
buf = np.frombuffer(f.read(8 * nb_pts * 3), dtype=np.float64)
self.pts = np.array(buf.copy(), dtype=np.float64).reshape([nb_pts, 3])
elif which_type == 32:
buf = np.frombuffer(f.read(4 * nb_pts * 3), dtype=np.float32)
self.pts = np.array(buf.copy(), dtype=np.float32).reshape([nb_pts, 3]).astype(np.float64)
nb_tri = int.from_bytes(f.read(4), 'little')
buf = np.frombuffer(f.read(4 * nb_tri * 3), dtype=np.uint32)
self.tri = np.array(buf.copy(), dtype=np.uint32).reshape([nb_tri, 3]).astype(np.int32)
self.validate_format()
self.find_minmax(True)
# ----------------------------------------------------------------
# File I/O — DXF
# ----------------------------------------------------------------
[docs]
def import_dxf(self, fn: str):
"""Import 3DFACE entities from a DXF file."""
import ezdxf
if not path.exists(fn):
logging.warning(_('File not found !') + ' ' + fn)
return
doc = ezdxf.readfile(fn)
msp = doc.modelspace()
xyz = []
for e in msp:
if e.dxftype() == "3DFACE":
xyz += [e[0], e[1], e[2]]
xyz_u, indices = np.unique(np.array(xyz), return_inverse=True, axis=0)
triangles = indices.reshape([int(len(indices) / 3), 3])
self.tri = triangles.astype(np.int32)
self.pts = xyz_u
self.validate_format()
# ====================================================================
# vectorpropertiesModel
# ====================================================================
[docs]
class vectorpropertiesModel:
"""Pure-data vector properties: colors, widths, legend, persistence.
Stores all visual/styling attributes for a :class:`vector` and
handles file I/O (*save* / *load*) without any wx or OpenGL
dependency.
The GUI subclass :class:`~wolfhece.PyVertexvectors.vectorproperties`
adds ``Wolf_Param`` dialog creation, ``genericImagetexture`` management,
and property-editing callbacks.
"""
def __init__(self, lines: list[str] = None, parent=None) -> None:
"""Initialise vector properties.
Compatibility with VB6/Fortran file formats — modify with care.
:param lines: property lines read from a ``.vecz`` file.
:param parent: owning :class:`vector` (or ``None``).
"""
if lines is None:
lines = []
# Defaults
self.color = 0
self.width = 1
self.style = 1
self.closed = False
self.filled = False
self.legendvisible = False
self.transparent = False
self.alpha = 0
self.flash = False
self.legendtext = ''
self.legendrelpos = 5
self.legendx = 0.
self.legendy = 0.
self.legendbold = False
self.legenditalic = False
self.legendfontname = 'Arial'
self.legendfontsize = 10
self.legendcolor = 0
self.legendunderlined = False
self.used = True
self.extrude = False
[docs]
self.plot_indices = False
[docs]
self.plot_lengths = False
if len(lines) > 0:
try:
line1 = lines[0].split(',')
line2 = lines[1].split(',')
self.color = int(line1[0])
self.width = int(float(line1[1]))
self.style = int(line1[2])
self.closed = line1[3] == '#TRUE#'
self.filled = line1[4] == '#TRUE#'
self.legendvisible = line1[5] == '#TRUE#'
self.transparent = line1[6] == '#TRUE#'
self.alpha = int(line1[7])
self.flash = line1[8] == '#TRUE#'
self.legendtext = line2[0]
self.legendrelpos = int(line2[1])
self.legendx = float(line2[2])
self.legendy = float(line2[3])
self.legendbold = line2[4] == '#TRUE#'
self.legenditalic = line2[5] == '#TRUE#'
self.legendfontname = str(line2[6])
self.legendfontsize = int(line2[7])
self.legendcolor = int(line2[8])
self.legendunderlined = line2[9] == '#TRUE#'
self.used = lines[2] == '#TRUE#'
except Exception as e:
logging.warning('Problem with vector properties - line format not correct - {}'.format(e))
self.init_extra()
# ----------------------------------------------------------------
# Item access (key-value storage)
# ----------------------------------------------------------------
def __getitem__(self, key: str) -> Union[int, float, bool, str]:
"""Get a named value."""
if key in self._values:
return self._values[key]
return None
def __setitem__(self, key: str, value: Union[int, float, bool, str]):
"""Set a named value."""
self._values[key] = value
[docs]
def add_values(self, values: dict):
"""Merge a mapping of values into the internal storage."""
if values is None:
return
self._values.update(values)
[docs]
def set_shared_values(self, shared_values: dict):
"""Attach an existing values mapping by reference.
This is useful when multiple vectors must share identical attribute
values (e.g. split parts of a multipart geometry).
"""
if shared_values is None:
shared_values = {}
self._values = shared_values
[docs]
def value_keys(self) -> set:
"""Return the set of stored value keys."""
return set(self._values.keys())
[docs]
def values_dict(self) -> dict:
"""Return a shallow copy of stored values."""
return dict(self._values)
[docs]
def values_view(self) -> dict:
"""Return the internal values mapping for read-only access.
Callers must not mutate the returned dictionary directly.
"""
return self._values
# ----------------------------------------------------------------
# Width-profile helpers
# ----------------------------------------------------------------
[docs]
def _parse_width_profile(self, raw) -> list | None:
"""Convert a width profile string to a list or ``None``.
Two formats are accepted:
* **Dense CSV** – ``"1.0,1.5,2.0"`` — one multiplier per vertex,
stored as ``list[float]``.
* **Sparse** – ``"0:1.0,10:2.0,20:0.5"`` — ``index:value`` pairs,
stored as ``list[tuple[int, float]]``. Unspecified vertices are
filled by linear interpolation (constant beyond the endpoints).
Both formats are auto-detected from the presence of ``:`` tokens.
"""
if raw is None:
return None
if isinstance(raw, list):
if len(raw) == 0:
return None
if isinstance(raw[0], (tuple, list)):
return [(int(p[0]), float(p[1])) for p in raw]
return [float(v) for v in raw]
txt = str(raw).strip()
if txt == '' or txt.lower() == 'none':
return None
tokens = [t.strip() for t in txt.split(',') if t.strip()]
if not tokens:
return None
if any(':' in t for t in tokens):
result = []
for token in tokens:
if ':' in token:
left, right = token.split(':', 1)
result.append((int(left.strip()), float(right.strip())))
return result if result else None
else:
vals = [float(t) for t in tokens]
return vals if vals else None
# ----------------------------------------------------------------
# Color helpers
# ----------------------------------------------------------------
[docs]
def set_color_from_value(self, key: str, cmap: wolfpalette | str | Colormap | cm.ScalarMappable,
vmin: float = 0., vmax: float = 1.):
"""Derive the drawing color from a stored value and a colormap."""
if key in self._values:
val = self._values[key]
try:
val = float(val)
except (TypeError, ValueError):
logging.warning('Value {} is not a numeric value'.format(val))
self.color = 0
return
val_adim = (val - vmin) / (vmax - vmin)
if isinstance(cmap, wolfpalette):
self.color = getIfromRGB(cmap.get_rgba_oneval(val))
elif isinstance(cmap, str):
cmap = cm.get_cmap(cmap)
self.color = getIfromRGB(cmap(val_adim, bytes=True))
elif isinstance(cmap, Colormap):
self.color = getIfromRGB(cmap(val_adim, bytes=True))
elif isinstance(cmap, cm.ScalarMappable):
self.color = getIfromRGB(cmap.to_rgba(val, bytes=True)[:3])
else:
self.color = 0
[docs]
def set_color(self, color: int | RGB | str | list | tuple):
"""Set the drawing color from various input formats."""
if isinstance(color, int):
self.color = color
elif isinstance(color, (list, tuple)) and len(color) == 3:
self.color = getIfromRGB(color)
elif isinstance(color, RGB):
self.color = color.int_format
elif isinstance(color, str):
if color.startswith('#'):
self.color = Colors.hexa_to_int(color)
else:
self.color = Colors()[color].int_format
else:
logging.warning('Invalid color type: {}'.format(type(color)))
self.color = 0
# ----------------------------------------------------------------
# Extra properties (legend geometry, attached image path, …)
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# File I/O — VB6/Fortran format
# ----------------------------------------------------------------
[docs]
def save(self, f: io.TextIOWrapper):
"""Write properties in VB6/Fortran-compatible format."""
line1 = str(self.color) + ',' + str(self.width) + ',' + str(self.style)
line1 += ',#TRUE#' if self.closed else ',#FALSE#'
line1 += ',#TRUE#' if self.filled else ',#FALSE#'
line1 += ',#TRUE#' if self.legendvisible else ',#FALSE#'
line1 += ',#TRUE#' if self.transparent else ',#FALSE#'
line1 += ',' + str(self.alpha)
line1 += ',#TRUE#' if self.flash else ',#FALSE#'
f.write(line1 + '\n')
line1 = (self.legendtext + ',' + str(self.legendrelpos)
+ ',' + str(self.legendx) + ',' + str(self.legendy))
line1 += ',#TRUE#' if self.legendbold else ',#FALSE#'
line1 += ',#TRUE#' if self.legenditalic else ',#FALSE#'
line1 += ',' + self.legendfontname + ',' + str(self.legendfontsize) + ',' + str(self.legendcolor)
line1 += ',#TRUE#' if self.legendunderlined else ',#FALSE#'
f.write(line1 + '\n')
f.write(('#TRUE#' if self.used else '#FALSE#') + '\n')
# ----------------------------------------------------------------
# Image path property
# ----------------------------------------------------------------
@property
[docs]
def image_path(self) -> Path:
"""Path of the attached image."""
return Path(self.attachedimage) if self.attachedimage is not None else Path()
@image_path.setter
def image_path(self, value: Path):
self.attachedimage = Path(value)
# ----------------------------------------------------------------
# Font conversion helpers
# ----------------------------------------------------------------
[docs]
def _convert_fontname2int(self, fontname: str) -> int:
if fontname.lower() in 'arial.ttf':
return 1
elif fontname.lower() in 'helvetica.ttf':
return 2
elif fontname.lower() in 'sanserif.ttf':
return 3
return 1
[docs]
def _convert_int2fontname(self, id: int) -> str:
if id == 1:
return 'Arial'
elif id == 2:
return 'Helvetica'
elif id == 3:
return 'SanSerif'
return 'Arial'
# ----------------------------------------------------------------
# Copy properties
# ----------------------------------------------------------------
[docs]
def egalize(self, other):
"""Copy visual attributes from another properties object or vector."""
# Avoid circular import at module level
from ._vector import vector
if isinstance(other, vectorpropertiesModel):
src = other
elif isinstance(other, vector):
src = other.myprop
else:
return
self.color = src.color
self.width = src.width
self.style = src.style
self.closed = src.closed
self.filled = src.filled
self.legendvisible = src.legendvisible
self.transparent = src.transparent
self.alpha = src.alpha
self.flash = src.flash
# ----------------------------------------------------------------
# GUI hook (no-op in model)
# ----------------------------------------------------------------
[docs]
def update_myprops(self):
"""No-op in model. GUI subclass populates the Wolf_Param dialog."""
pass
[docs]
def update_image_texture(self):
"""No-op in model. GUI subclass updates texture coordinates."""
pass
# ====================================================================
# vectorModel
# ====================================================================
[docs]
class vectorModel:
"""Pure-data vector: vertices, properties, geometry, I/O.
Holds a list of :class:`wolfvertex` and a :class:`vectorpropertiesModel`,
plus geometry helpers (shapely, interpolation, etc.) without any
OpenGL, wx, or matplotlib dependency.
GUI subclass ``vector`` in :mod:`wolfhece.PyVertexvectors` adds
plotting, tree-widget integration and OpenGL display-list support.
"""
[docs]
parentzone: "zoneModel"
[docs]
myvertices: list[wolfvertex]
[docs]
myprop: "vectorpropertiesModel"
def __init__(self, lines: list = [], is2D=True, name='NoName',
parentzone: "zoneModel" = None, fromshapely=None,
fromnumpy: np.ndarray = None, fromlist: list = None) -> None:
self.myname = ''
self.parentzone = parentzone
self.xmin = -99999.
self.ymin = -99999.
self.xmax = -99999.
self.ymax = -99999.
[docs]
self._lengthparts2D = None
[docs]
self._lengthparts3D = None
[docs]
self._center_segments = None
[docs]
self._simplified_geometry = False
if type(lines) == list:
if len(lines) > 0:
self.myname = lines[0]
tmp_nbvertices = int(lines[1])
self.myvertices = []
if is2D:
for i in range(tmp_nbvertices):
try:
curx, cury = lines[2 + i].split()
except:
curx, cury = lines[2 + i].split(',')
curvert = wolfvertex(float(curx), float(cury))
self.myvertices.append(curvert)
else:
for i in range(tmp_nbvertices):
try:
curx, cury, curz = lines[2 + i].split()
except:
curx, cury, curz = lines[2 + i].split(',')
curvert = wolfvertex(float(curx), float(cury), float(curz))
self.myvertices.append(curvert)
self.myprop = self._make_properties(lines[tmp_nbvertices + 2:], parent=self)
if name != '' and self.myname == '':
self.myname = name
self.myvertices = []
self.myprop = self._make_properties(parent=self)
[docs]
self._cache_vertices = None
[docs]
self._move_start = None
[docs]
self._rotation_center = None
[docs]
self._rotation_step = None
[docs]
self._linestring = None
[docs]
self._vertex_kdtree = None
[docs]
self._vertex_kdtree_xy = None
[docs]
self.add_zdatum = False
[docs]
self.add_sdatum = False
if fromshapely is not None:
self.import_shapelyobj(fromshapely)
if fromnumpy is not None:
self.add_vertices_from_array(fromnumpy)
if fromlist is not None:
self.add_vertices_from_list(fromlist)
# ----------------------------------------------------------------
# Factory methods (overridden in GUI subclass)
# ----------------------------------------------------------------
[docs]
def _make_properties(self, *args, **kwargs) -> "vectorpropertiesModel":
"""Create a properties object. GUI subclass returns ``vectorproperties``."""
return vectorpropertiesModel(*args, **kwargs)
[docs]
def _make_vector(self, **kwargs) -> "vectorModel":
"""Create a sibling vector. GUI subclass returns ``vector``."""
return vectorModel(**kwargs)
[docs]
def _make_zone(self, **kwargs) -> "zoneModel":
"""Create a sibling zone. GUI subclass returns ``zone``."""
return zoneModel(**kwargs)
@classmethod
[docs]
def make_from_shapely(cls, shapelyobj, **kwargs) -> "vectorModel":
"""Factory method to create a vector from a Shapely geometry."""
newvec = cls(fromshapely=shapelyobj, **kwargs)
return newvec
@classmethod
[docs]
def make_from_numpy(cls, array: np.ndarray, **kwargs) -> "vectorModel":
"""Factory method to create a vector from a NumPy array of shape (N, 2) or (N, 3)."""
newvec = cls(fromnumpy=array, **kwargs)
return newvec
@classmethod
[docs]
def make_from_list(cls, lst: list, **kwargs) -> "vectorModel":
"""Factory method to create a vector from a list of coordinate tuples."""
newvec = cls(fromlist=lst, **kwargs)
return newvec
# ----------------------------------------------------------------
# GUI hooks (no-op in model)
# ----------------------------------------------------------------
[docs]
def _on_vertices_changed(self):
"""Hook called when vertices are modified.
Invalidates the cached Shapely geometries. The GUI subclass
overrides this (calling ``super()``) to also invalidate the
parent zone's OpenGL display list.
"""
self._invalidate_spatial_index()
self.reset_linestring()
# ----------------------------------------------------------------
# Properties
# ----------------------------------------------------------------
@property
[docs]
def closed(self) -> bool:
"""Whether the vector is closed (first vertex == last vertex)."""
return self.myprop.closed
@closed.setter
def closed(self, value: bool):
"""Set the closed state of the vector.
:param value: ``True`` to mark the vector as closed.
"""
self.myprop.closed = value
@property
[docs]
def nbvertices(self) -> int:
"""Number of vertices in the vector."""
return len(self.myvertices)
@property
[docs]
def used(self) -> bool:
"""Whether the vector is marked as used/active."""
return self.myprop.used
@used.setter
def used(self, value: bool):
"""Set the used/active state of the vector.
:param value: ``True`` to mark the vector as used.
"""
self.myprop.used = value
@property
[docs]
def polygon(self) -> Polygon:
"""Shapely ``Polygon`` built from the vertices (lazily created)."""
if self._polygon is None:
self.prepare_shapely(linestring=False)
return self._polygon
@property
[docs]
def linestring(self) -> LineString:
"""Shapely ``LineString`` built from the vertices (lazily created)."""
if self._linestring is None:
self.prepare_shapely(polygon=False)
return self._linestring
@property
[docs]
def has_interior(self) -> bool:
"""Whether any vertex is flagged as *not in use*, indicating interior boundaries."""
not_in_use = [curvert for curvert in self.myvertices if not curvert.in_use]
return len(not_in_use) > 0
@property
[docs]
def nb_interiors(self) -> int:
"""Number of interior boundaries (holes) derived from unused vertices."""
not_in_use = [curvert for curvert in self.myvertices if not curvert.in_use]
return len(not_in_use) // 2 if self.myprop.filled else len(not_in_use)
@property
[docs]
def centroid(self) -> Point:
"""Centroid of the polygon formed by the vertices (Shapely ``Point``)."""
return self.polygon.centroid
@property
[docs]
def surface(self) -> float:
"""Area of the polygon (0 if the vector is not closed)."""
if self.closed:
return self.polygon.area
else:
return 0.
@property
[docs]
def area(self) -> float:
"""Alias for :pyattr:`surface`."""
return self.surface
# ----------------------------------------------------------------
# Coordinate array properties
# ----------------------------------------------------------------
@property
[docs]
def z(self) -> np.ndarray:
"""Z coordinates of all vertices as a 1-D array.
If ``add_zdatum`` is ``True``, the ``zdatum`` offset is added.
"""
z = np.asarray([curvert.z for curvert in self.myvertices])
if self.add_zdatum:
z += self.zdatum
return z
@property
[docs]
def x(self) -> np.ndarray:
"""X coordinates of all vertices as a 1-D array."""
return np.asarray([curvert.x for curvert in self.myvertices])
@property
[docs]
def y(self) -> np.ndarray:
"""Y coordinates of all vertices as a 1-D array."""
return np.asarray([curvert.y for curvert in self.myvertices])
@property
[docs]
def xy(self) -> np.ndarray:
"""XY coordinates as a ``(N, 2)`` array."""
return np.asarray([[curvert.x, curvert.y] for curvert in self.myvertices])
@property
[docs]
def xz(self) -> np.ndarray:
"""XZ coordinates as a ``(N, 2)`` array."""
return np.asarray([[curvert.x, curvert.z] for curvert in self.myvertices])
@property
[docs]
def xyz(self) -> np.ndarray:
"""XYZ coordinates as a ``(N, 3)`` array (respects ``zdatum``)."""
return self.asnparray3d()
@property
[docs]
def i(self) -> np.ndarray:
"""In-use flags of all vertices as a 1-D boolean array."""
return np.asarray([curvert.in_use for curvert in self.myvertices])
@property
[docs]
def xyzi(self) -> np.ndarray:
"""XYZI coordinates as a ``(N, 4)`` array (x, y, z, in_use)."""
x = self.x; y = self.y; z = self.z; i = self.i
return np.column_stack((x, y, z, i))
@property
[docs]
def xyi(self) -> np.ndarray:
"""XYI coordinates as a ``(N, 3)`` array (x, y, in_use)."""
return np.asarray([[curvert.x, curvert.y, curvert.in_use] for curvert in self.myvertices])
@property
[docs]
def sz_curvi(self) -> tuple[np.ndarray, np.ndarray]:
"""Tuple ``(s, z)`` of curvilinear abscissa and altitude arrays."""
return self.get_sz()
@property
[docs]
def s_curvi(self) -> np.ndarray:
"""Curvilinear abscissa (cumulative 2-D distance from first vertex)."""
sz = self.get_sz()
return sz[0]
@x.setter
def x(self, new_x: np.ndarray | list):
"""Set X coordinates for all vertices.
:param new_x: Array or list of new X values (length must equal ``nbvertices``).
"""
if isinstance(new_x, list):
new_x = np.array(new_x)
if len(new_x) != self.nbvertices:
logging.warning(_('New x values have not the same length as the number of vertices'))
return
for curvert, newx in zip(self.myvertices, new_x):
curvert.x = newx
self._on_vertices_changed()
@y.setter
def y(self, new_y: np.ndarray | list):
"""Set Y coordinates for all vertices.
:param new_y: Array or list of new Y values (length must equal ``nbvertices``).
"""
if isinstance(new_y, list):
new_y = np.array(new_y)
if len(new_y) != self.nbvertices:
logging.warning(_('New y values have not the same length as the number of vertices'))
return
for curvert, newy in zip(self.myvertices, new_y):
curvert.y = newy
self._on_vertices_changed()
@z.setter
def z(self, new_z: np.ndarray | float | list):
"""Set Z coordinates for all vertices.
:param new_z: Array, list, scalar float, or ``WolfArray`` of new Z values.
If a ``WolfArray`` is given the Z values are sampled at each vertex position.
Respects ``zdatum`` when ``add_zdatum`` is ``True``.
"""
from ..wolf_array import WolfArray
if isinstance(new_z, (int, float)):
new_z = np.full(self.nbvertices, new_z, dtype=float)
elif isinstance(new_z, WolfArray):
wa = new_z
new_z = []
for curvert in self.myvertices:
i, j = wa.xy2ij(curvert.x, curvert.y)
if i > 0 and j > 0 and i < wa.nbx and j < wa.nby:
new_z.append(wa.array[i, j])
if isinstance(new_z, list):
new_z = np.array(new_z)
if len(new_z) != self.nbvertices:
logging.warning(_('New z values have not the same length as the number of vertices'))
return
if self.add_zdatum:
for curvert, newz in zip(self.myvertices, new_z):
curvert.z = newz - self.zdatum
else:
for curvert, newz in zip(self.myvertices, new_z):
curvert.z = newz
self._on_vertices_changed()
@xyz.setter
def xyz(self, new_xyz: np.ndarray | list):
"""Set XYZ coordinates for all vertices.
:param new_xyz: ``(N, 3)`` array or list of ``(x, y, z)`` tuples.
"""
if isinstance(new_xyz, list):
new_xyz = np.array(new_xyz)
if len(new_xyz) != self.nbvertices:
logging.warning(_('New xyz values have not the same length as the number of vertices'))
return
if self.add_zdatum:
for curvert, newxyz in zip(self.myvertices, new_xyz):
curvert.x = newxyz[0]
curvert.y = newxyz[1]
curvert.z = newxyz[2] - self.zdatum
else:
for curvert, newxyz in zip(self.myvertices, new_xyz):
curvert.x = newxyz[0]
curvert.y = newxyz[1]
curvert.z = newxyz[2]
self._on_vertices_changed()
@xy.setter
def xy(self, new_xy: np.ndarray | list):
"""Set XY coordinates for all vertices.
:param new_xy: ``(N, 2)`` array or list of ``(x, y)`` tuples.
"""
if isinstance(new_xy, list):
new_xy = np.array(new_xy)
if len(new_xy) != self.nbvertices:
logging.warning(_('New xy values have not the same length as the number of vertices'))
return
for curvert, newxy in zip(self.myvertices, new_xy):
curvert.x = newxy[0]
curvert.y = newxy[1]
self._on_vertices_changed()
@xz.setter
def xz(self, new_xz: np.ndarray | list):
"""Set XZ coordinates for all vertices.
:param new_xz: ``(N, 2)`` array or list of ``(x, z)`` tuples.
Respects ``zdatum`` when ``add_zdatum`` is ``True``.
"""
if isinstance(new_xz, list):
new_xz = np.array(new_xz)
if len(new_xz) != self.nbvertices:
logging.warning(_('New xz values have not the same length as the number of vertices'))
return
if self.add_zdatum:
for curvert, newxz in zip(self.myvertices, new_xz):
curvert.x = newxz[0]
curvert.z = newxz[1] - self.zdatum
else:
for curvert, newxz in zip(self.myvertices, new_xz):
curvert.x = newxz[0]
curvert.z = newxz[1]
self._on_vertices_changed()
@xyzi.setter
def xyzi(self, new_xyzi: np.ndarray | list):
"""Set XYZI coordinates for all vertices.
:param new_xyzi: ``(N, 4)`` array or list of ``(x, y, z, in_use)`` tuples.
"""
if isinstance(new_xyzi, list):
new_xyzi = np.array(new_xyzi)
if len(new_xyzi) != self.nbvertices:
logging.warning(_('New xyzi values have not the same length as the number of vertices'))
return
for curvert, newxyzi in zip(self.myvertices, new_xyzi):
curvert.x = newxyzi[0]
curvert.y = newxyzi[1]
curvert.z = newxyzi[2] - self.zdatum if self.add_zdatum else newxyzi[2]
curvert.in_use = newxyzi[3]
self._on_vertices_changed()
@xyi.setter
def xyi(self, new_xyi: np.ndarray | list):
"""Set XYI coordinates for all vertices.
:param new_xyi: ``(N, 3)`` array or list of ``(x, y, in_use)`` tuples.
"""
if isinstance(new_xyi, list):
new_xyi = np.array(new_xyi)
if len(new_xyi) != self.nbvertices:
logging.warning(_('New xyi values have not the same length as the number of vertices'))
return
for curvert, newxyi in zip(self.myvertices, new_xyi):
curvert.x = newxyi[0]
curvert.y = newxyi[1]
curvert.in_use = newxyi[2]
self._on_vertices_changed()
@i.setter
def i(self, new_i: np.ndarray | list):
"""Set in-use flags for all vertices.
:param new_i: Array or list of boolean/numeric in-use flags.
"""
if isinstance(new_i, list):
new_i = np.array(new_i)
if len(new_i) != self.nbvertices:
logging.warning(_('New i values have not the same length as the number of vertices'))
return
for curvert, newi in zip(self.myvertices, new_i):
curvert.in_use = newi
self._on_vertices_changed()
@sz_curvi.setter
def sz_curvi(self, sz_new: np.ndarray | list):
"""Set Z values by curvilinear interpolation.
:param sz_new: ``(M, 2)`` array of ``(s, z)`` pairs used as interpolation table.
"""
if isinstance(sz_new, list):
sz_new = np.array(sz_new)
f = interp1d(sz_new[:, 0], sz_new[:, 1], bounds_error=False, fill_value='extrapolate')
s = self.s_curvi
for idx, curvert in enumerate(self.myvertices):
curvert.z = f(s[idx])
self._on_vertices_changed()
@s_curvi.setter
def s_curvi(self, new_s: np.ndarray | list):
"""Redistribute vertices to new curvilinear abscissae along the linestring.
:param new_s: Array or list of new curvilinear abscissa values.
"""
if isinstance(new_s, list):
new_s = np.array(new_s)
if len(new_s) != self.nbvertices:
logging.warning(_('New s values have not the same length as the number of vertices'))
return
poly = self.linestring
for idx, curvert in enumerate(self.myvertices):
curvert.x, curvert.y = poly.interpolate(new_s[idx]).xy
self._on_vertices_changed()
# ----------------------------------------------------------------
# Dunder methods
# ----------------------------------------------------------------
def __add__(self, other: "vectorModel") -> "vectorModel":
"""Concatenate two vectors into a new one.
:param other: The other vector to append.
:returns: A new vector containing the vertices of both operands.
"""
if not isinstance(other, vectorModel):
raise TypeError("Can only add vector to vector")
new_vector = self._make_vector(name=self.myname + '_' + other.myname)
new_vector.myvertices = self.myvertices.copy() + other.myvertices.copy()
return new_vector
def __str__(self):
"""Return the vector name."""
return self.myname
def __len__(self):
"""Return the number of vertices."""
return self.nbvertices
def __iter__(self) -> wolfvertex:
"""Iterate over the vertices."""
return iter(self.myvertices)
def __getitem__(self, ndx: int) -> wolfvertex:
"""Return the vertex at index *ndx* (supports negative indexing).
:param ndx: Vertex index.
"""
if ndx >= 0 and ndx < self.nbvertices:
return self.myvertices[ndx]
elif ndx < 0 and abs(ndx) <= self.nbvertices:
return self.myvertices[ndx + self.nbvertices]
else:
logging.warning(_('Index out of range'))
def __setitem__(self, ndx: int, value: wolfvertex):
"""Replace the vertex at index *ndx*.
:param ndx: Vertex index (supports negative indexing).
:param value: New ``wolfvertex`` instance.
"""
if ndx >= 0 and ndx < self.nbvertices:
self.myvertices[ndx] = value
self._on_vertices_changed()
elif ndx < 0 and abs(ndx) <= self.nbvertices:
self.myvertices[ndx + self.nbvertices] = value
self._on_vertices_changed()
else:
logging.warning(_('Index out of range'))
def __delitem__(self, ndx: int):
"""Delete the vertex at index *ndx*.
:param ndx: Vertex index (supports negative indexing).
"""
if ndx >= 0 and ndx < self.nbvertices:
self.myvertices.pop(ndx)
self._on_vertices_changed()
elif ndx < 0 and abs(ndx) <= self.nbvertices:
self.myvertices.pop(ndx + self.nbvertices)
self._on_vertices_changed()
else:
logging.warning(_('Index out of range'))
def __del__(self):
"""Cleanup: release cached geometry."""
self.reset_linestring()
# ----------------------------------------------------------------
# Value / colour / legend helpers
# ----------------------------------------------------------------
[docs]
def add_value(self, key: str, value: Union[int, float, bool, str]):
"""Store an arbitrary key/value pair in the vector properties.
:param key: Property key.
:param value: Property value.
"""
self.myprop[key] = value
[docs]
def get_value(self, key: str) -> Union[int, float, bool, str]:
"""Retrieve a stored property value by key.
:param key: Property key.
:returns: The stored value.
"""
return self.myprop[key]
[docs]
def set_color_from_value(self, key: str, cmap: wolfpalette | Colormap | cm.ScalarMappable,
vmin: float = 0., vmax: float = 1.):
"""Set the vector colour from a stored property value using a colormap.
:param key: Property key whose numeric value drives the colour.
:param cmap: Colour map (``wolfpalette``, Matplotlib ``Colormap``, or ``ScalarMappable``).
:param vmin: Minimum value for colour mapping.
:param vmax: Maximum value for colour mapping.
"""
self.myprop.set_color_from_value(key, cmap, vmin, vmax)
[docs]
def set_color(self, color: int | RGB | str | list | tuple):
"""Set the vector rendering colour.
:param color: Colour as an integer, ``RGB`` instance, name string, or ``(r, g, b)`` tuple/list.
"""
self.myprop.set_color(color)
[docs]
def set_linewidth(self, linewidth: int | float):
"""Set the rendering line width.
:param linewidth: Line width in pixels.
"""
self.myprop.width = int(linewidth)
[docs]
def set_alpha(self, alpha: int):
"""Set the transparency level.
:param alpha: Alpha value (0 = fully transparent, 255 = opaque).
"""
self.myprop.alpha = alpha
self.myprop.transparent = alpha != 0
[docs]
def set_filled(self, filled: bool):
"""Enable or disable polygon filling.
:param filled: ``True`` to render as a filled polygon.
"""
self.myprop.filled = filled
[docs]
def set_props_from_other(self, other):
"""Copy visual properties from another vector's properties.
:param other: Source ``vectorpropertiesModel`` instance.
"""
self.myprop.egalize(other)
[docs]
def set_legend_text(self, text: str):
if text is _('Not used'):
pass
elif text == _('name'):
self.myprop.legendtext = self.myname
elif text == _('first z'):
if self.nbvertices > 0:
self.myprop.legendtext = str(self.myvertices[0].z)
else:
self.myprop.legendtext = ''
elif text == _('length2D'):
self.myprop.legendtext = str(self.length2D)
elif text == _('length3D'):
self.myprop.legendtext = str(self.length3D)
elif text == _('id'):
if self.parentzone is not None:
self.myprop.legendtext = str(self.parentzone.myvectors.index(self))
else:
self.myprop.legendtext = ''
else:
self.myprop.legendtext = str(text)
self.myprop.update_myprops()
[docs]
def set_legend_color(self, color: int | RGB | str | list | tuple):
if isinstance(color, int):
self.myprop.legendcolor = color
elif isinstance(color, RGB):
self.myprop.legendcolor = color.int_format
elif isinstance(color, str):
self.myprop.legendcolor = Colors()[color].int_format
elif isinstance(color, (list, tuple)):
self.myprop.legendcolor = getIfromRGB(color)
self.myprop.update_myprops()
[docs]
def set_legend_text_from_value(self, key: str, visible: bool = True):
"""Set the legend text from a stored property value.
:param key: Property key whose value becomes the label.
:param visible: Whether to make the legend visible.
"""
if key in self.myprop._values:
self.myprop.legendtext = str(self.myprop._values[key])
self.myprop.legendvisible = visible
self.myprop.update_myprops()
[docs]
def set_legend_position(self, x: str | float, y: str | float):
"""Set the legend display position.
:param x: X position as a float or keyword (*'median'*, *'mean'*, *'min'*,
*'max'*, *'first'*, *'last'*).
:param y: Y position as a float or keyword (same options as *x*).
"""
if isinstance(x, str):
if x in [_('Not used'), 'Not used']:
pass
elif x.lower() in [_('median'), 'median']:
xy = self.xy
self.myprop.legendx = np.median(xy[:, 0])
elif x.lower() in [_('mean'), 'mean']:
xy = self.xy
self.myprop.legendx = np.mean(xy[:, 0])
elif x.lower() in [_('min'), 'min']:
xy = self.xy
self.myprop.legendx = np.min(xy[:, 0])
elif x.lower() in [_('max'), 'max']:
xy = self.xy
self.myprop.legendx = np.max(xy[:, 0])
elif x.lower() in [_('first'), 'first']:
self.myprop.legendx = self.myvertices[0].x
elif x.lower() in [_('last'), 'last']:
self.myprop.legendx = self.myvertices[-1].x
elif x.lower() in [_('centroid'), 'centroid']:
centroid = self.centroid
self.myprop.legendx = centroid.x
else:
try:
self.myprop.legendx = float(x)
except ValueError:
logging.warning(_('Invalid x position for legend: {}').format(x))
elif isinstance(x, float):
self.myprop.legendx = x
if isinstance(y, str):
if y in [_('Not used'), 'Not used']:
pass
elif y.lower() in [_('median'), 'median']:
xy = self.xy
self.myprop.legendy = np.median(xy[:, 1])
elif y.lower() in [_('mean'), 'mean']:
xy = self.xy
self.myprop.legendy = np.mean(xy[:, 1])
elif y.lower() in [_('min'), 'min']:
xy = self.xy
self.myprop.legendy = np.min(xy[:, 1])
elif y.lower() in [_('max'), 'max']:
xy = self.xy
self.myprop.legendy = np.max(xy[:, 1])
elif y.lower() in [_('first'), 'first']:
self.myprop.legendy = self.myvertices[0].y
elif y.lower() in [_('last'), 'last']:
self.myprop.legendy = self.myvertices[-1].y
elif y.lower() in [_('centroid'), 'centroid']:
centroid = self.centroid
self.myprop.legendy = centroid.y
else:
try:
self.myprop.legendy = float(y)
except ValueError:
logging.warning(_('Invalid y position for legend: {}').format(y))
elif isinstance(y, float):
self.myprop.legendy = y
self.myprop.update_myprops()
[docs]
def set_legend_to_centroid(self, text: str = '', visible: bool = True):
"""Place the legend at the polygon centroid.
:param text: Legend text (defaults to the vector name).
:param visible: Whether to make the legend visible.
"""
self.myprop.legendvisible = visible
centroid = self.centroid
self.myprop.legendx = centroid.x
self.myprop.legendy = centroid.y
self.myprop.legendtext = text if text else self.myname
[docs]
def set_legend_visible(self, visible: bool = True):
"""Show or hide the legend.
:param visible: ``True`` to display the legend.
"""
self.myprop.legendvisible = visible
[docs]
def set_legend_position_to_centroid(self):
"""Move the legend position to the polygon centroid."""
centroid = self.centroid
self.myprop.legendx = centroid.x
self.myprop.legendy = centroid.y
# ----------------------------------------------------------------
# Shader-property helpers
# ----------------------------------------------------------------
[docs]
def set_glow(self, enabled: bool = True, width: float = 0.4,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.4):
"""Configure line glow effect.
:param enabled: Enable or disable glow.
:param width: Glow width as fraction of half-width (0–1).
:param color: RGB colour tuple with components in [0, 1].
:param alpha: Glow opacity in [0, 1].
"""
self.myprop.glow_enabled = enabled
self.myprop.glow_width = width
self.myprop.glow_color = (color[0], color[1], color[2], alpha)
self.myprop.update_myprops()
[docs]
def set_dash(self, enabled: bool = True, dash_length: float = 10.0,
gap_length: float = 5.0):
"""Configure dashed line rendering.
:param enabled: Enable or disable dashes.
:param dash_length: Dash length in world units.
:param gap_length: Gap length in world units.
"""
self.myprop.dash_enabled = enabled
self.myprop.dash_length = dash_length
self.myprop.gap_length = gap_length
self.myprop.update_myprops()
[docs]
def set_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure line animation.
:param mode: Animation mode (0=none, 1=pulse, 2=marching ants).
:param speed: Speed multiplier.
"""
self.myprop.anim_mode = mode
self.myprop.anim_speed = speed
self.myprop.update_myprops()
[docs]
def set_fill_animation(self, mode: int = 1, speed: float = 1.0,
center_index: int = 0,
start_angle: float = 90.0):
"""Configure filled-polygon animation.
:param mode: Animation mode (0=none, 1=pulse, 2=sweep, 3=ripple,
4=radial progressive, 5=clockwise clock, 6=counter-clockwise clock).
:param speed: Speed multiplier.
:param center_index: Vertex index used as the radial/clock centre.
:param start_angle: Clock fill start angle in degrees.
"""
self.myprop.fill_anim_mode = mode
self.myprop.fill_anim_speed = speed
self.myprop.fill_anim_center_index = int(center_index)
self.myprop.fill_anim_start_angle = float(start_angle)
self.myprop.update_myprops()
[docs]
def set_fill_color(self, color=None):
"""Set the independent fill colour for filled polygons.
:param color: RGB tuple ``(r, g, b)`` with values in ``[0, 255]``,
packed int, or ``None`` to inherit the vector's main colour.
"""
if color is None:
self.myprop.fill_color = None
elif isinstance(color, int):
self.myprop.fill_color = color
else:
self.myprop.fill_color = getIfromRGB(color)
self.myprop.update_myprops()
[docs]
def set_contour_color(self, color=None, width: float = 1.0,
enabled: bool = True):
"""Set the contour/stroke colour for filled polygons.
Activates a second polyline pass drawn on top of the filled area.
Set *enabled=False* or *color=None* to disable the contour.
:param color: RGB tuple ``(r, g, b)`` with values in ``[0, 255]``,
packed int, or ``None`` to inherit the vector's main colour.
:param width: Contour line width (same unit as the vector width).
:param enabled: Whether to draw the contour pass at all.
"""
if color is None:
self.myprop.contour_color = None
elif isinstance(color, int):
self.myprop.contour_color = color
else:
self.myprop.contour_color = getIfromRGB(color)
self.myprop.contour_width = float(width)
self.myprop.contour_enabled = enabled
self.myprop.update_myprops()
[docs]
def set_join_style(self, style: int = 0, size: float = 0.5,
size_mode: int = 0):
"""Configure join rendering at polyline corners.
:param style: Join style (0=miter, 1=bevel, 2=round, 3=fillet).
:param size: Join size (fraction or world distance depending on *size_mode*).
:param size_mode: 0=fraction (adimensional), 1=world distance (metres).
"""
self.myprop.join_style = style
self.myprop.join_size = size
self.myprop.join_size_mode = size_mode
self.myprop.update_myprops()
[docs]
def set_width_in_pixels(self, pixels: bool = True):
"""Toggle pixel-based vs world-unit line width.
:param pixels: ``True`` for pixel widths, ``False`` for world units.
"""
self.myprop.width_in_pixels = pixels
self.myprop.update_myprops()
[docs]
def set_text_along(self, text: str, enabled: bool = True,
offset: float = 0.0, perp: float = 0.0,
alignment: str = 'center'):
"""Configure text displayed along the polyline.
:param text: Text to render along the path.
:param enabled: Enable or disable text-along.
:param offset: Shift start along the polyline (world units).
:param perp: Perpendicular offset (world units, + left).
:param alignment: Text alignment ('left', 'center', 'right').
"""
self.myprop.text_along_enabled = enabled
self.myprop.text_along_text = text
self.myprop.text_along_offset = offset
self.myprop.text_along_perp = perp
self.myprop.text_along_alignment = alignment
self.myprop.update_myprops()
[docs]
def set_legend_glow(self, enabled: bool = True, width: float = 0.15,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.5):
"""Configure legend text glow effect.
:param enabled: Enable or disable legend glow.
:param width: SDF threshold offset (e.g. 0.15).
:param color: RGB colour tuple with components in [0, 1].
:param alpha: Glow opacity in [0, 1].
"""
self.myprop.legend_glow_enabled = enabled
self.myprop.legend_glow_width = width
self.myprop.legend_glow_color = (color[0], color[1], color[2], alpha)
self.myprop.update_myprops()
[docs]
def set_legend_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure legend text animation.
:param mode: Animation mode (0=none, 1=pulse, 2=wave, 3=typewriter).
:param speed: Speed multiplier.
"""
self.myprop.legend_anim_mode = mode
self.myprop.legend_anim_speed = speed
self.myprop.update_myprops()
[docs]
def set_legend_style(self, smoothing: float = 1.0,
alignment: str = 'left',
line_spacing: float = 1.2):
"""Fine-tune legend SDF text rendering.
:param smoothing: SDF anti-aliasing multiplier.
:param alignment: Text alignment ('left', 'center', 'right').
:param line_spacing: Multiline spacing multiplier.
"""
self.myprop.legend_smoothing = smoothing
self.myprop.legend_alignment = alignment
self.myprop.legend_line_spacing = line_spacing
self.myprop.update_myprops()
[docs]
def highlighting(self, rgb=(255, 0, 0), linewidth=3):
"""Highlight the vector with a temporary colour and line width.
:param rgb: Highlight colour as ``(r, g, b)`` tuple.
:param linewidth: Highlight line width.
"""
self._oldcolor = self.myprop.color
self._oldwidth = self.myprop.color
self.myprop.color = getIfromRGB(rgb)
self.myprop.width = linewidth
[docs]
def withdrawal(self):
"""Restore the colour and line width saved by :pymeth:`highlighting`."""
try:
self.myprop.color = self._oldcolor
self.myprop.width = self._oldwidth
except:
self.myprop.color = 0
self.myprop.width = 1
# ----------------------------------------------------------------
# Vertex management
# ----------------------------------------------------------------
[docs]
def add_vertex(self, addedvert: Union[list[wolfvertex], wolfvertex]):
"""Append one or more vertices.
:param addedvert: A single ``wolfvertex`` or a list of vertices to append.
"""
if type(addedvert) is list:
for curvert in addedvert:
assert curvert is not None
assert isinstance(curvert, wolfvertex)
self.myvertices.append(curvert)
else:
assert(addedvert is not None)
assert isinstance(addedvert, wolfvertex)
self.myvertices.append(addedvert)
self._on_vertices_changed()
[docs]
def add_vertices_from_array(self, xyz: np.ndarray):
"""Append vertices from a NumPy array.
:param xyz: ``(N, 2)`` or ``(N, 3)`` array of coordinates.
"""
assert isinstance(xyz, np.ndarray), "xyz must be a numpy array of shape (nb_vert, 2 or 3)"
if xyz.dtype == np.int32:
xyz = xyz.astype(np.float64)
if xyz.shape[1] == 3:
for cur in xyz:
self.myvertices.append(wolfvertex(cur[0], cur[1], cur[2]))
elif xyz.shape[1] == 2:
for cur in xyz:
self.myvertices.append(wolfvertex(cur[0], cur[1]))
self._on_vertices_changed()
[docs]
def add_vertices_from_list(self, xyz: list):
"""Append vertices from a list of coordinate tuples.
:param xyz: List of ``(x, y)`` or ``(x, y, z)`` tuples.
"""
for cur in xyz:
if len(cur) == 3:
self.myvertices.append(wolfvertex(cur[0], cur[1], cur[2]))
elif len(cur) == 2:
self.myvertices.append(wolfvertex(cur[0], cur[1]))
self._on_vertices_changed()
[docs]
def import_shapelyobj(self, obj):
"""Replace vertices from a Shapely ``LineString`` or ``Polygon``.
:param obj: Shapely geometry object.
"""
self.myvertices = []
if isinstance(obj, LineString):
self.is2D = not obj.has_z
xyz = np.array(obj.coords)
self.add_vertices_from_array(xyz)
elif isinstance(obj, Polygon):
self.is2D = not obj.has_z
xyz = np.array(obj.exterior.coords)
self.add_vertices_from_array(xyz)
self.close_force()
else:
logging.warning(_('Object type {} not supported -- Update "import_shapelyobj"').format(type(obj)))
[docs]
def count(self):
return
[docs]
def close_force(self):
"""Close the vector by appending the first vertex if needed."""
is_open = not ((self.myvertices[-1] is self.myvertices[0]) or
(self.myvertices[-1].x == self.myvertices[0].x and
self.myvertices[-1].y == self.myvertices[0].y))
if not self.is2D:
is_open = is_open or self.myvertices[-1].z != self.myvertices[0].z
if is_open:
self.add_vertex(self.myvertices[0])
self.closed = True
[docs]
def force_to_close(self):
"""Alias for :pymeth:`close_force`."""
self.close_force()
[docs]
def reverse(self):
"""Reverse the order of vertices in-place."""
self.myvertices.reverse()
self._on_vertices_changed()
[docs]
def reset(self):
"""Remove all vertices and reset caches."""
self.myvertices = []
self._on_vertices_changed()
[docs]
def append(self, other: "vectorModel", merge_type: Literal['link', 'copy'] = 'link'):
"""Append vertices from another vector.
:param other: Source vector.
:param merge_type: *'link'* to share vertex references, *'copy'* for deep copies.
"""
if merge_type == 'link':
self.myvertices.extend(other.myvertices)
elif merge_type == 'copy':
self.myvertices.extend(other.myvertices.copy())
else:
logging.warning(_('Merge type not supported'))
self.update_lengths()
self._on_vertices_changed()
# ----------------------------------------------------------------
# I/O
# ----------------------------------------------------------------
[docs]
def _set_limits(self):
"""Compute and store the bounding-box limits."""
self.find_minmax()
self._mylimits = ((self.xmin, self.xmax), (self.ymin, self.ymax))
[docs]
def _nblines(self):
return self.nbvertices + 5
[docs]
def save(self, f):
"""Write the vector definition (name, vertices, properties) to an open file.
:param f: Writable text file handle.
"""
f.write(self.myname + '\n')
f.write(str(self.nbvertices) + '\n')
force3D = False
if self.parentzone is not None:
force3D = self.parentzone.parent.force3D
if self.is2D and not force3D:
for curvert in self.myvertices:
f.write(f'{curvert.x},{curvert.y}' + '\n')
else:
for curvert in self.myvertices:
f.write(f'{curvert.x},{curvert.y},{curvert.z}' + '\n')
self.myprop.save(f)
[docs]
def set_z(self, new_z: np.ndarray):
"""Deprecated — use the ``z`` property setter instead."""
warnings.warn(_('This method is deprecated, use the z property instead.'), DeprecationWarning, stacklevel=2)
self.z = new_z
# ----------------------------------------------------------------
# Cache / move / rotate
# ----------------------------------------------------------------
[docs]
def set_cache(self):
"""Cache the current vertex coordinates and bounds for incremental transforms."""
self._cache_vertices = self.asnparray3d()
self._cache_bounds = np.array([[self.xmin, self.xmax], [self.ymin, self.ymax]])
[docs]
def clear_cache(self):
"""Clear the cached vertex coordinates and bounds."""
self._cache_vertices = None
self._cache_bounds = None
[docs]
def move(self, deltax: float, deltay: float, use_cache: bool = True, inplace: bool = True):
"""Translate the vector.
:param deltax: Displacement along X.
:param deltay: Displacement along Y.
:param use_cache: Use cached coordinates for incremental moves.
:param inplace: Modify in place; if ``False`` return a moved copy.
:returns: The moved vector (self or a new copy).
"""
if self._move_step is not None:
deltax = np.round(deltax / self._move_step) * self._move_step
deltay = np.round(deltay / self._move_step) * self._move_step
if inplace:
if use_cache and self._cache_vertices is None:
self.set_cache()
if use_cache:
self.xy = self._cache_vertices[:, :2] + np.array([deltax, deltay])
self.xmin = self._cache_bounds[0, 0] + deltax
self.xmax = self._cache_bounds[0, 1] + deltax
self.ymin = self._cache_bounds[1, 0] + deltay
self.ymax = self._cache_bounds[1, 1] + deltay
else:
for curvert in self.myvertices:
curvert.x += deltax
curvert.y += deltay
self.xmin += deltax
self.xmax += deltax
self.ymin += deltay
self.ymax += deltay
if self.myprop.textureimage is not None:
self.myprop.textureimage.xmin = self.xmin
self.myprop.textureimage.xmax = self.xmax
self.myprop.textureimage.ymin = self.ymin
self.myprop.textureimage.ymax = self.ymax
return self
else:
new_vector = self.deepcopy_vector(self.myname + '_moved')
new_vector.move(deltax, deltay, inplace=True)
return new_vector
[docs]
def rotate(self, angle: float, center: tuple[float, float] = (0., 0.),
use_cache: bool = True, inplace: bool = True):
"""Rotate the vector around a centre point.
:param angle: Rotation angle in degrees.
:param center: Centre of rotation ``(x, y)``.
:param use_cache: Use cached coordinates for incremental rotations.
:param inplace: Modify in place; if ``False`` return a rotated copy.
:returns: The rotated vector (self or a new copy).
"""
if inplace:
if use_cache and self._cache_vertices is None:
self.set_cache()
if use_cache:
locxy = self._cache_vertices[:, :2] - np.array(center)
rotation_matrix = np.array([
[np.cos(np.radians(angle)), -np.sin(np.radians(angle))],
[np.sin(np.radians(angle)), np.cos(np.radians(angle))]])
self.xy = np.dot(locxy, rotation_matrix) + np.array(center)
else:
for curvert in self.myvertices:
curvert.rotate(angle, center)
self.find_minmax()
if self.myprop.textureimage is not None:
self.myprop.textureimage.xmin = self.xmin
self.myprop.textureimage.xmax = self.xmax
self.myprop.textureimage.ymin = self.ymin
self.myprop.textureimage.ymax = self.ymax
return self
else:
new_vector = self.deepcopy_vector(self.myname + '_rotated')
new_vector.rotate(angle, center, inplace=True)
return new_vector
[docs]
def rotate_xy(self, x: float, y: float, use_cache: bool = True, inplace: bool = True):
if self._rotation_center is None:
logging.error('No rotation center defined -- set it before rotating by this routine')
return self
angle = np.degrees(np.arctan2(-(y - self._rotation_center[1]), x - self._rotation_center[0]))
if self._rotation_step is not None:
angle = np.round(angle / self._rotation_step) * self._rotation_step
return self.rotate(angle, center=self._rotation_center, use_cache=use_cache, inplace=inplace)
[docs]
def get_mapviewer(self):
"""Return the parent map-viewer widget, or ``None``."""
if self.parentzone is not None:
return self.parentzone.get_mapviewer()
else:
return None
# ----------------------------------------------------------------
# Spatial queries / bounds
# ----------------------------------------------------------------
[docs]
def find_minmax(self, only_firstlast: bool = False):
"""Compute the bounding box of the vertices.
:param only_firstlast: If ``True``, consider only the first and last vertex.
"""
if self.nbvertices > 0:
if only_firstlast:
self.xmin = min(self.myvertices[0].x, self.myvertices[-1].x)
self.ymin = min(self.myvertices[0].y, self.myvertices[-1].y)
self.xmax = max(self.myvertices[0].x, self.myvertices[-1].x)
self.ymax = max(self.myvertices[0].y, self.myvertices[-1].y)
else:
self.xmin = min(vert.x for vert in self.myvertices)
self.ymin = min(vert.y for vert in self.myvertices)
self.xmax = max(vert.x for vert in self.myvertices)
self.ymax = max(vert.y for vert in self.myvertices)
self.update_image_texture()
else:
self.xmin = -99999.
self.ymin = -99999.
self.xmax = -99999.
self.ymax = -99999.
[docs]
def update_image_texture(self):
"""Synchronise the image-texture bounds with the current bounding box."""
self.myprop.update_image_texture()
[docs]
def get_bounds(self, grid: float = None):
"""Return the bounding box as ``((xmin, ymin), (xmax, ymax))``.
:param grid: If given, snap bounds to the nearest grid step.
"""
if grid is None:
return ((self.xmin, self.ymin), (self.xmax, self.ymax))
else:
xmin = float(np.int32(self.xmin / grid) * grid)
ymin = float(np.int32(self.ymin / grid) * grid)
xmax = float(np.ceil(self.xmax / grid) * grid)
ymax = float(np.ceil(self.ymax / grid) * grid)
return ((xmin, ymin), (xmax, ymax))
[docs]
def get_bounds_xx_yy(self, grid: float = None):
"""Return the bounding box as ``((xmin, xmax), (ymin, ymax))``.
:param grid: If given, snap bounds to the nearest grid step.
"""
if grid is None:
return ((self.xmin, self.xmax), (self.ymin, self.ymax))
else:
xmin = float(np.int32(self.xmin / grid) * grid)
ymin = float(np.int32(self.ymin / grid) * grid)
xmax = float(np.ceil(self.xmax / grid) * grid)
ymax = float(np.ceil(self.ymax / grid) * grid)
return ((xmin, xmax), (ymin, ymax))
[docs]
def intersects_bounds(self, xmin: float = None, ymin: float = None,
xmax: float = None, ymax: float = None) -> bool:
"""Return whether the vector bounding box intersects the given bbox.
If bounds are not provided, the vector is considered intersecting.
Sentinel unset bounds are treated as intersecting to avoid false negatives.
"""
if any(val is None for val in (xmin, ymin, xmax, ymax)):
return True
if min(self.xmin, self.ymin, self.xmax, self.ymax) <= -99999.0:
return True
return not (
self.xmax < xmin or self.xmin > xmax or
self.ymax < ymin or self.ymin > ymax
)
[docs]
def get_mean_vertex(self, asshapelyPoint=False):
"""Return the mean position of all vertices.
:param asshapelyPoint: If ``True``, return a Shapely ``Point``;
otherwise a ``wolfvertex``.
"""
if self.closed or (self.myvertices[0].x == self.myvertices[-1].x and self.myvertices[0].y == self.myvertices[-1].y):
x_mean = np.mean([curvert.x for curvert in self.myvertices[:-1]])
y_mean = np.mean([curvert.y for curvert in self.myvertices[:-1]])
z_mean = np.mean([curvert.z for curvert in self.myvertices[:-1]])
else:
x_mean = np.mean([curvert.x for curvert in self.myvertices])
y_mean = np.mean([curvert.y for curvert in self.myvertices])
z_mean = np.mean([curvert.z for curvert in self.myvertices])
if asshapelyPoint:
return Point(x_mean, y_mean, z_mean)
else:
return wolfvertex(x_mean, y_mean, z_mean)
[docs]
def isinside(self, x:float, y:float) -> bool:
"""Test whether the point ``(x, y)`` is inside the polygon.
:param x: X coordinate.
:param y: Y coordinate.
"""
if self.nbvertices == 0:
return False
poly = self.polygon
inside2 = poly.contains(Point([x, y]))
return inside2
[docs]
def contains(self, x: float, y: float) -> bool:
"""Alias for :pymeth:`isinside`."""
return self.isinside(x, y)
[docs]
def _invalidate_spatial_index(self):
"""Drop cached spatial indices built from the current vertices."""
self._vertex_kdtree = None
self._vertex_kdtree_xy = None
[docs]
def _get_vertex_kdtree(self):
"""Return the lazy-built KDTree used for nearest-vertex queries."""
if self._vertex_kdtree is None:
coords = self.asnparray()
if len(coords) == 0:
return None, None
self._vertex_kdtree_xy = coords
self._vertex_kdtree = cKDTree(coords)
return self._vertex_kdtree, self._vertex_kdtree_xy
[docs]
def _get_nearest_search_geometry(self):
"""Return the geometry used for nearest-geometry queries."""
if self.nbvertices == 0:
return None
if self.nbvertices == 1:
return Point(self.myvertices[0].x, self.myvertices[0].y)
if self.closed and self.nbvertices >= 3:
return self.polygon
return self.linestring
[docs]
def distance_to_geometry(self, x: float, y: float) -> float:
"""Return the distance from ``(x, y)`` to the vector geometry."""
geometry = self._get_nearest_search_geometry()
if geometry is None:
return np.inf
return float(geometry.distance(Point(x, y)))
[docs]
def find_nearest_vertex(self, x, y) -> wolfvertex:
"""Return the vertex nearest to ``(x, y)``.
:param x: X coordinate.
:param y: Y coordinate.
"""
tree, coords = self._get_vertex_kdtree()
if tree is None or coords is None:
return None
distance, index = tree.query(np.asarray([x, y], dtype=np.float64), k=1)
index = int(index)
if np.isfinite(distance):
duplicate_indices = np.where(np.all(coords == coords[index], axis=1))[0]
if duplicate_indices.size > 0:
index = int(duplicate_indices[0])
return self.myvertices[index]
[docs]
def insert_nearest_vert(self, x, y) -> wolfvertex:
"""Insert a new vertex at ``(x, y)`` between the two nearest existing vertices.
:param x: X coordinate of the new vertex.
:param y: Y coordinate of the new vertex.
:returns: The inserted ``wolfvertex``.
"""
xy = Point(x, y)
mynp = self.asnparray()
mp = MultiPoint(mynp)
ls = LineString(mynp)
nearmp = nearest_points(mp, xy)
nearls = nearest_points(ls, xy)
smp = ls.project(nearmp[0])
sls = ls.project(nearls[0])
indexmp = np.argwhere(mynp == np.asarray([nearmp[0].x, nearmp[0].y]))[0, 0]
if indexmp == 0 and self.closed:
if sls <= ls.length and sls >= ls.project(Point([self.myvertices[-2].x, self.myvertices[-2].y])):
smp = ls.length
indexmp = self.nbvertices - 1
else:
indexmp = 1
else:
if sls >= smp:
indexmp += 1
myvert = wolfvertex(x, y)
self.myvertices.insert(indexmp, myvert)
self._on_vertices_changed()
return self.myvertices[indexmp]
[docs]
def insert_vertex_at_s(self, s: float, z: float = None, tolerance: float = 1e-3) -> wolfvertex:
"""Insert a vertex at curvilinear abscissa *s*.
:param s: Distance along the linestring.
:param z: Optional Z value for the new vertex.
:param tolerance: Distance threshold to reuse an existing vertex.
:returns: The inserted (or reused) vertex.
"""
ls = self.linestring
if s < 0 or s > ls.length:
logging.error(_('Distance s={} is out of bounds (0, {}) -- cannot insert vertex').format(s, ls.length))
return None
point = ls.interpolate(s)
if z is None:
z = point.z
for curvert in self.myvertices:
dist = np.sqrt((curvert.x - point.x) ** 2 + (curvert.y - point.y) ** 2)
if dist <= tolerance:
curvert.z = z
return curvert
newvert = wolfvertex(point.x, point.y, z)
self.myvertices.insert(np.searchsorted([v for v in self.s_curvi], self.linestring.project(point)), newvert)
self._on_vertices_changed()
return newvert
[docs]
def project_vertex_onto_vector_and_insert(self, xy: Point | wolfvertex, tolerance: float = 1e-3) -> wolfvertex:
"""Project a point onto the linestring and insert a new vertex there.
:param xy: Point to project (Shapely ``Point`` or ``wolfvertex``).
:param tolerance: Distance threshold to reuse an existing vertex.
:returns: The inserted (or reused) vertex.
"""
ls = self.linestring
if isinstance(xy, wolfvertex):
point = Point(xy.x, xy.y)
else:
point = xy
s = ls.project(point)
projected_point = ls.interpolate(s)
for curvert in self.myvertices:
dist = np.sqrt((curvert.x - projected_point.x) ** 2 + (curvert.y - projected_point.y) ** 2)
if dist <= tolerance:
return curvert
newvert = wolfvertex(projected_point.x, projected_point.y, projected_point.z)
self.myvertices.insert(np.searchsorted(self.linestring.project(point), [v for v in self.myvertices]), newvert)
self._on_vertices_changed()
return newvert
[docs]
def get_normal_segments(self, update: bool = False) -> np.ndarray:
"""Return outward unit normals for each segment as a ``(N-1, 2)`` array.
:param update: Force recomputation even if cached.
"""
if self._normals is not None and not update:
return self._normals
if self.nbvertices < 2:
return None
xy = self.asnparray()
delta = xy[1:] - xy[:-1]
norm = np.sqrt(delta[:, 0] ** 2 + delta[:, 1] ** 2)
normals = np.zeros_like(delta)
notnull = np.where(norm != 0)
normals[notnull, 0] = -delta[notnull, 1] / norm[notnull]
normals[notnull, 1] = delta[notnull, 0] / norm[notnull]
self._normals = normals
return normals
[docs]
def get_center_segments(self, update: bool = False) -> np.ndarray:
"""Return segment mid-points as a ``(N-1, 2)`` array.
:param update: Force recomputation even if cached.
"""
if self._center_segments is not None and not update:
return self._center_segments
if self.nbvertices < 2:
return None
xy = self.asnparray()
centers = (xy[1:] + xy[:-1]) / 2.
self._center_segments = centers
return centers
[docs]
def select_points_inside(self, xy: cloud_vertices | np.ndarray) -> list[bool]:
"""Return a list of booleans indicating which points are inside the polygon.
:param xy: ``cloud_vertices`` or ``(N, 2)`` array of XY coordinates.
"""
self.prepare_shapely(True)
if isinstance(xy, cloud_vertices):
xy = xy.get_xyz()[:, 0:2]
inside = [self.polygon.contains(Point(curxy)) for curxy in xy]
return inside
[docs]
def get_first_point_inside(self, xy: cloud_vertices | np.ndarray):
"""Return the first point that lies inside the polygon.
:param xy: ``cloud_vertices`` or ``(N, 2)`` array of XY coordinates.
:returns: ``(x, y)`` tuple or ``None``.
"""
self.prepare_shapely(True)
if isinstance(xy, cloud_vertices):
xy = xy.get_xyz()[:, 0:2]
for curxy in xy:
if self.polygon.contains(Point(curxy)):
return float(curxy[0]), float(curxy[1])
return None
[docs]
def split_cloud(self, cloud_to_split: cloud_vertices):
"""Split a point cloud into *inside* and *outside* subsets.
:param cloud_to_split: Point cloud to partition.
:returns: ``(cloud_inside, cloud_outside)`` tuple.
"""
inside = self.select_points_inside(cloud_to_split)
cloud_inside = cloud_vertices(idx='inside_' + cloud_to_split.idx)
cloud_outside = cloud_vertices(idx='outside_' + cloud_to_split.idx)
vertices = cloud_to_split.get_vertices()
for idx, (locinside, curvert) in enumerate(zip(inside, vertices)):
if locinside:
cloud_inside.add_vertex(curvert)
else:
cloud_outside.add_vertex(curvert)
return cloud_inside, cloud_outside
[docs]
def check_if_closed(self) -> bool:
"""Return ``True`` if the vector is geometrically closed."""
return not self.check_if_open()
[docs]
def check_if_open(self) -> bool:
"""Return ``True`` if the vector is open; also updates ``self.closed``."""
is_open = not ((self.myvertices[-1] is self.myvertices[0]) or
(self.myvertices[-1].x == self.myvertices[0].x and
self.myvertices[-1].y == self.myvertices[0].y))
if not self.is2D:
is_open = is_open or self.myvertices[-1].z != self.myvertices[0].z
self.closed = not is_open
return is_open
[docs]
def check_if_interior_exists(self):
"""Detect repeated XY coordinates and flag them as interior (``in_use = False``)."""
xy = self.xy
if self.closed and (xy[0, 0] == xy[-1, 0] and xy[0, 1] == xy[-1, 1]):
xy = xy[:-1]
xy_unique, inverse, count = np.unique(xy, return_inverse=True, return_counts=True, axis=0)
duplicate_found = False
if xy.shape[0] != xy_unique.shape[0]:
duplicate_indices = np.where(count > 1)[0]
duplicate_indices = np.where(np.isin(inverse, duplicate_indices))[0]
diff = np.diff(duplicate_indices)
for i in range(len(diff)):
if diff[i] == 1:
self.myvertices[duplicate_indices[i + 1]].in_use = False
duplicate_found = True
if duplicate_found:
self._on_vertices_changed()
# ----------------------------------------------------------------
# Shapely conversions
# ----------------------------------------------------------------
[docs]
def asshapely_pol(self) -> Polygon:
"""Convert vertex list to a Shapely 2-D ``Polygon``."""
if self.nbvertices < 3:
return Polygon()
coords = self.asnparray()
return Polygon(coords)
[docs]
def asshapely_pol3D(self) -> Polygon:
"""Convert vertex list to a Shapely 3-D ``Polygon``."""
if self.nbvertices < 3:
return Polygon()
coords = self.asnparray3d()
return Polygon(coords)
[docs]
def asshapely_ls3d(self) -> LineString:
"""Convert vertex list to a Shapely 3-D ``LineString``."""
coords = self.asnparray3d()
return LineString(coords)
[docs]
def asshapely_ls(self) -> LineString:
"""Convert vertex list to a Shapely 3-D ``LineString`` (returns empty if < 2 vertices)."""
coords = self.asnparray3d()
if len(coords) < 2:
logging.warning(_('Not enough coordinates to create a LineString -- return an empty LineString'))
return LineString()
return LineString(coords)
[docs]
def asshapely_mp(self) -> MultiPoint:
"""Convert vertex list to a Shapely ``MultiPoint``."""
multi = self.asnparray3d()
return MultiPoint(multi)
[docs]
def asnparray(self) -> np.ndarray:
"""Return vertices as a ``(N, 2)`` NumPy array of ``(x, y)``."""
return np.asarray(list([vert.x, vert.y] for vert in self.myvertices))
[docs]
def asnparray3d(self) -> np.ndarray:
"""Return vertices as a ``(N, 3)`` NumPy array of ``(x, y, z)`` (respects ``zdatum``)."""
xyz = np.asarray(list([vert.x, vert.y, vert.z] for vert in self.myvertices))
if self.add_zdatum:
xyz[:, 2] += self.zdatum
return xyz
[docs]
def prepare_shapely(self, prepare_shapely: bool = True, linestring: bool = True, polygon: bool = True):
"""Build and optionally prepare Shapely linestring and/or polygon geometries.
:param prepare_shapely: Whether to call Shapely ``prepare()`` for faster operations.
:param linestring: Build the linestring geometry.
:param polygon: Build the polygon geometry.
"""
self.reset_linestring()
if linestring:
self._linestring = self.asshapely_ls()
if polygon:
self._polygon = self.asshapely_pol()
if prepare_shapely:
if linestring:
if not self._linestring.is_empty:
prepare(self._linestring)
if polygon:
if not self._polygon.is_empty:
prepare(self._polygon)
[docs]
def reset_linestring(self):
"""Destroy cached Shapely linestring and polygon (with prepared versions)."""
if hasattr(self, '_linestring'):
if self._linestring is not None:
if is_prepared(self._linestring):
destroy_prepared(self._linestring)
self._linestring = None
if hasattr(self, '_polygon'):
if self._polygon is not None:
if is_prepared(self._polygon):
destroy_prepared(self._polygon)
self._polygon = None
[docs]
def intersection(self, vec2: "vectorModel" = None, eval_dist=False, norm=False, force_single=False):
"""Compute the intersection with another vector.
:param vec2: Other vector to intersect with.
:param eval_dist: Also return projected distance(s) on this linestring.
:param norm: Normalise distances to ``[0, 1]``.
:param force_single: If multi-result, keep only the first geometry.
:returns: Intersection geometry, or ``(geom, dist)`` when *eval_dist* is ``True``.
"""
ls1 = self.linestring
ls2 = vec2.linestring
myinter = ls1.intersection(ls2)
if isinstance(myinter, MultiPoint):
if force_single:
myinter = myinter.geoms[0]
elif isinstance(myinter, MultiLineString):
if force_single:
myinter = Point(myinter.geoms[0].xy[0][0], myinter.geoms[0].xy[1][0])
if myinter.is_empty:
if eval_dist:
return None, None
else:
return None
if eval_dist:
mydists = ls1.project(myinter, normalized=norm)
return myinter, mydists
else:
return myinter
[docs]
def intersects(self, x: float, y: float) -> bool:
"""Test whether the linestring intersects the point ``(x, y)``.
:param x: X coordinate.
:param y: Y coordinate.
"""
point = Point(x, y)
return self.linestring.intersects(point)
[docs]
def aligned_with(self, x: float, y: float, tolerance: float = 1e-6) -> bool:
"""Test whether the point ``(x, y)`` lies on the linestring within *tolerance*.
:param x: X coordinate.
:param y: Y coordinate.
:param tolerance: Distance tolerance.
"""
return self.intersects(x, y) and self.linestring.distance(Point(x, y)) < tolerance
[docs]
def buffer(self, distance: float, resolution: int = 16, inplace: bool = True) -> "vectorModel":
"""Return the buffered polygon.
:param distance: Buffer distance.
:param resolution: Number of segments to approximate a quarter-circle.
:param inplace: Modify this vector; if ``False`` return a new one.
"""
poly = self.polygon
buffered = poly.buffer(distance, resolution=resolution)
if inplace:
self.import_shapelyobj(buffered)
self.reset_linestring()
return self
else:
newvec = self._make_vector(fromshapely=buffered, name=self.myname)
return newvec
[docs]
def projectontrace(self, trace: "vectorModel") -> "vectorModel":
"""Project each vertex onto a trace and return a new 2-D vector ``(s, z)``.
:param trace: Reference trace vector.
:returns: New vector with ``(s, z)`` coordinates.
"""
tracels: LineString
tracels = trace.linestring
xyz = self.asnparray3d()
all_s = [tracels.project(Point(cur[0], cur[1])) for cur in xyz]
newvec = self._make_vector(name=_('Projection on ') + trace.myname)
for s, (x, y, z) in zip(all_s, xyz):
newvec.add_vertex(wolfvertex(s, z))
return newvec
[docs]
def parallel_offset(self, distance=5., side: Literal['left', 'right'] = 'left') -> "vectorModel | None":
"""Create a parallel offset of the linestring.
:param distance: Offset distance.
:param side: Side of the offset (*'left'* or *'right'*).
:returns: New vector or ``None`` on failure.
"""
if self.nbvertices < 2:
return None
myls = self.asshapely_ls()
from shapely.geometry import JOIN_STYLE
mypar = myls.parallel_offset(distance=abs(distance), side=side, join_style=JOIN_STYLE.round)
if mypar.geom_type == 'MultiLineString':
return None
if len(mypar.coords) > 0:
newvec = self._make_vector(name='par' + str(distance) + '_' + self.myname, parentzone=self.parentzone)
for x, y in mypar.coords:
xy = Point(x, y)
curs = mypar.project(xy, True)
curz = myls.interpolate(curs, True).z
newvert = wolfvertex(x, y, curz)
newvec.add_vertex(newvert)
return newvec
else:
return None
# ----------------------------------------------------------------
# Length / interpolation
# ----------------------------------------------------------------
[docs]
def verify_s_ascending(self) -> tuple[bool, list[int]]:
"""Check and fix non-ascending curvilinear abscissae by swapping vertices.
:returns: ``(correction_applied, list_of_corrected_indices)``.
"""
s, z = self.get_sz(cumul=False)
correction = False
where = []
for i in range(self.nbvertices - 1):
if s[i] > s[i + 1]:
correction = True
where.append(i + 1)
x = self.myvertices[i].x
y = self.myvertices[i].y
self.myvertices[i].x = self.myvertices[i + 1].x
self.myvertices[i].y = self.myvertices[i + 1].y
self.myvertices[i + 1].x = x
self.myvertices[i + 1].y = y
return correction, where
[docs]
def get_s2d(self) -> np.ndarray:
"""Return cumulative 2-D curvilinear distances from the first vertex."""
s2d = np.zeros(self.nbvertices)
for k in range(1, self.nbvertices):
s2d[k] = s2d[k - 1] + self.myvertices[k - 1].dist2D(self.myvertices[k])
return s2d
[docs]
def get_s3d(self) -> np.ndarray:
"""Return cumulative 3-D curvilinear distances from the first vertex."""
s3d = np.zeros(self.nbvertices)
for k in range(1, self.nbvertices):
s3d[k] = s3d[k - 1] + self.myvertices[k - 1].dist3D(self.myvertices[k])
return s3d
[docs]
def get_sz(self, cumul=True) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(s, z)`` arrays of curvilinear abscissa and altitude.
:param cumul: If ``True``, cumulative inter-vertex distance;
otherwise direct distance from the first vertex.
"""
z = np.asarray([vert.z for vert in self.myvertices])
nb = len(z)
s = np.zeros(nb)
if cumul:
x1 = self.myvertices[0].x
y1 = self.myvertices[0].y
for i in range(nb - 1):
x2 = self.myvertices[i + 1].x
y2 = self.myvertices[i + 1].y
length = np.sqrt((x2 - x1) ** 2. + (y2 - y1) ** 2.)
s[i + 1] = s[i] + length
x1 = x2
y1 = y2
else:
for i in range(nb):
s[i] = self.myvertices[0].dist2D(self.myvertices[i])
if self.add_sdatum:
s += self.sdatum
if self.add_zdatum:
z += self.zdatum
return s, z
[docs]
def update_lengths(self):
"""Recompute 2-D and 3-D segment lengths and total lengths."""
if self.nbvertices < 2:
logging.warning(_('No enough vertices in vector to compute lenghts'))
return
self._lengthparts2D = np.zeros(self.nbvertices - 1)
self._lengthparts3D = np.zeros(self.nbvertices - 1)
for k in range(self.nbvertices - 1):
self._lengthparts2D[k] = self.myvertices[k].dist2D(self.myvertices[k + 1])
self._lengthparts3D[k] = self.myvertices[k].dist3D(self.myvertices[k + 1])
if self.closed and self.myvertices[0] != self.myvertices[-1]:
self._lengthparts2D[-1] = self.myvertices[-2].dist2D(self.myvertices[-1])
self._lengthparts3D[-1] = self.myvertices[-2].dist3D(self.myvertices[-1])
self.length2D = np.sum(self._lengthparts2D)
self.length3D = np.sum(self._lengthparts3D)
[docs]
def get_segment(self, s, is3D, adim=True, frombegin=True):
"""Find the segment containing the curvilinear abscissa *s*.
:param s: Curvilinear abscissa.
:param is3D: Use 3-D lengths.
:param adim: If ``True``, *s* is normalised to ``[0, 1]``.
:param frombegin: Search from the beginning; if ``False``, from the end.
:returns: ``(segment_index, cumulative_s, lengthparts)``.
"""
if self.length2D is None or self.length3D is None:
self.update_lengths()
else:
if len(self._lengthparts2D) != self.nbvertices - 1 or len(self._lengthparts3D) != self.nbvertices - 1:
self.update_lengths()
if is3D:
length = self.length3D
lengthparts = self._lengthparts3D
else:
length = self.length2D
lengthparts = self._lengthparts2D
cums = np.cumsum(lengthparts)
if adim:
if length == 0.:
logging.warning(_('Length of vector {} is zero, cannot compute segments').format(self.myname))
cums = cums.copy() / length
cums[-1] = 1.
lengthparts = lengthparts.copy() / length
if s > 1.:
s = 1.
if s < 0.:
s = 0.
else:
if s > length:
s = length
if s < 0.:
s = 0.
if frombegin:
k = 0
while s > cums[k] and k < self.nbvertices - 2:
k += 1
else:
k = self.nbvertices - 2
if k > 0:
while s < cums[k] and k > 0:
k -= 1
if (k == 0 and s < cums[0]) or s == cums[self.nbvertices - 2]:
pass
else:
k += 1
if k == len(cums):
k -= 1
return k, cums[k], lengthparts
[docs]
def interpolate(self, s: float, is3D: bool = True, adim: bool = True, frombegin: bool = True) -> wolfvertex:
"""Interpolate a point along the polyline at curvilinear abscissa *s*.
:param s: Curvilinear abscissa.
:param is3D: Use 3-D lengths.
:param adim: If ``True``, *s* is normalised to ``[0, 1]``.
:param frombegin: Search from the beginning.
:returns: An interpolated ``wolfvertex``.
"""
if self.length2D is None or self.length3D is None:
self.update_lengths()
if self.length2D == 0.:
pond = 0.
k = 0
else:
k, cums, lengthparts = self.get_segment(s, is3D, adim, frombegin)
pond = (cums - s) / lengthparts[k]
return wolfvertex(self.myvertices[k].x * pond + self.myvertices[k + 1].x * (1. - pond),
self.myvertices[k].y * pond + self.myvertices[k + 1].y * (1. - pond),
self.myvertices[k].z * pond + self.myvertices[k + 1].z * (1. - pond))
[docs]
def tangent_at_s(self, s: float, is3D: bool = True, adim: bool = True, frombegin: bool = True) -> wolfvertex:
"""Return the unit tangent vector at curvilinear abscissa *s*.
:param s: Curvilinear abscissa.
:param is3D: Use 3-D lengths.
:param adim: Normalised abscissa.
:param frombegin: Search direction.
"""
if self.length2D is None or self.length3D is None:
self.update_lengths()
if self.length2D == 0.:
return wolfvertex(0., 0., 0.)
else:
k, cums, lengthparts = self.get_segment(s, is3D, adim, frombegin)
dx = self.myvertices[k + 1].x - self.myvertices[k].x
dy = self.myvertices[k + 1].y - self.myvertices[k].y
dz = self.myvertices[k + 1].z - self.myvertices[k].z
length = np.sqrt(dx * dx + dy * dy + dz * dz)
if length == 0.:
return wolfvertex(0., 0., 0.)
return wolfvertex(dx / length, dy / length, dz / length)
[docs]
def normal_at_s(self, s: float, is3D: bool = True, adim: bool = True, frombegin: bool = True,
counterclockwise=True) -> wolfvertex:
"""Return the unit normal vector at curvilinear abscissa *s*.
:param s: Curvilinear abscissa.
:param is3D: Use 3-D lengths.
:param adim: Normalised abscissa.
:param frombegin: Search direction.
:param counterclockwise: Normal orientation.
"""
tangent = self.tangent_at_s(s, is3D, adim, frombegin)
normal = wolfvertex(-tangent.y, tangent.x, 0.)
length = np.sqrt(normal.x * normal.x + normal.y * normal.y + normal.z * normal.z)
if length == 0.:
return wolfvertex(0., 0., 0.)
if counterclockwise:
return wolfvertex(normal.x / length, normal.y / length, normal.z / length)
else:
return wolfvertex(-normal.x / length, -normal.y / length, -normal.z / length)
[docs]
def _refine2D(self, ds):
"""Refine the polyline in 2-D with a target spacing *ds*.
:param ds: Maximum distance between interpolated points.
:returns: List of Shapely ``Point`` objects.
"""
myls = self.asshapely_ls()
length = myls.length
nb = int(np.ceil(length / ds)) + 1
if length < 1:
length *= 1000.
alls = np.linspace(0, length, nb)
alls /= 1000.
else:
alls = np.linspace(0, length, nb, endpoint=True)
pts = [myls.interpolate(curs) for curs in alls]
return pts
[docs]
def substring(self, s1: float, s2: float, is3D: bool = True, adim: bool = True, eps: float = 1.e-2) -> "vectorModel":
"""Extract a sub-vector between curvilinear abscissae *s1* and *s2*.
:param s1: Start curvilinear abscissa.
:param s2: End curvilinear abscissa.
:param is3D: Use 3-D lengths.
:param adim: Normalised abscissae.
:param eps: Minimum length offset when *s1* == *s2*.
:returns: New sub-vector.
"""
if s1 == s2:
logging.debug(_('Substring with same start and end abscissa: s1={} s2={}').format(s1, s2))
s2 += eps
k1, cums1, lengthparts1 = self.get_segment(s1, is3D, adim, True)
k2, cums2, lengthparts2 = self.get_segment(s2, is3D, adim, False)
pond1 = max((cums1 - s1) / lengthparts1[k1], 0.) if lengthparts1[k1] > 0. else 0.
pond2 = min((cums2 - s2) / lengthparts2[k2], 1.) if lengthparts2[k2] > 0. else 1.
v1 = wolfvertex(self.myvertices[k1].x * pond1 + self.myvertices[k1 + 1].x * (1. - pond1),
self.myvertices[k1].y * pond1 + self.myvertices[k1 + 1].y * (1. - pond1),
self.myvertices[k1].z * pond1 + self.myvertices[k1 + 1].z * (1. - pond1))
v2 = wolfvertex(self.myvertices[k2].x * pond2 + self.myvertices[k2 + 1].x * (1. - pond2),
self.myvertices[k2].y * pond2 + self.myvertices[k2 + 1].y * (1. - pond2),
self.myvertices[k2].z * pond2 + self.myvertices[k2 + 1].z * (1. - pond2))
newvec = self._make_vector(name='substr')
newvec.add_vertex(v1)
if s1 <= s2:
if is3D:
for k in range(k1 + 1, k2 + 1):
if self.myvertices[k].dist3D(newvec.myvertices[-1]) != 0.:
newvec.add_vertex(self.myvertices[k])
else:
for k in range(k1 + 1, k2 + 1):
if self.myvertices[k].dist2D(newvec.myvertices[-1]) != 0.:
newvec.add_vertex(self.myvertices[k])
else:
if is3D:
for k in range(k1 + 1, k2 + 1, -1):
if self.myvertices[k].dist3D(newvec.myvertices[-1]) != 0.:
newvec.add_vertex(self.myvertices[k])
else:
for k in range(k1 + 1, k2 + 1, -1):
if self.myvertices[k].dist2D(newvec.myvertices[-1]) != 0.:
newvec.add_vertex(self.myvertices[k])
if [v2.x, v2.y, v2.z] != [newvec.myvertices[-1].x, newvec.myvertices[-1].y, newvec.myvertices[-1].z]:
newvec.add_vertex(v2)
return newvec
[docs]
def split(self, ds, new=True):
"""Split the vector by inserting vertices at regular 3-D spacing *ds*.
:param ds: Target 3-D spacing.
:param new: If ``True``, add the new vector to the parent zone;
otherwise replace vertices in place.
"""
newvec = self._make_vector(name=self.myname + '_split', parentzone=self.parentzone)
self.update_lengths()
locds = ds / self.length3D
dist3d = np.concatenate([np.arange(0., 1., locds), np.cumsum(self._lengthparts3D) / self.length3D])
dist3d = np.unique(dist3d)
for curs in dist3d:
newvec.add_vertex(self.interpolate(curs, is3D=True, adim=True))
if new:
curzone = self.parentzone
if curzone is not None:
curzone.add_vector(newvec)
curzone._fill_structure()
else:
self.myvertices = newvec.myvertices
self.update_lengths()
[docs]
def cut(self, s: float, is3D: bool = True, adim: bool = True, frombegin: bool = True) -> "vectorModel":
"""Cut the vector at curvilinear abscissa *s* and return the detached part.
:param s: Curvilinear abscissa of the cut.
:param is3D: Use 3-D lengths.
:param adim: Normalised abscissa.
:param frombegin: Cut direction.
:returns: New vector containing the cut-off portion.
"""
newvec = self._make_vector(name=self.myname + '_cut', parentzone=self.parentzone)
self.parentzone.add_vector(newvec, update_struct=True)
k, cums, lengthparts = self.get_segment(s, is3D, adim, frombegin)
if frombegin:
newvec.myvertices = self.myvertices[:k + 1]
self.myvertices = self.myvertices[k:]
else:
newvec.myvertices = self.myvertices[k:]
self.myvertices = self.myvertices[:k + 1]
self.update_lengths()
newvec.update_lengths()
self._on_vertices_changed()
return newvec
[docs]
def simplify(self, tolerance: float = 1.0, preserve_topology: bool = True) -> "vectorModel":
"""Simplify the geometry using Shapely's Douglas-Peucker algorithm.
:param tolerance: Distance tolerance for simplification.
:param preserve_topology: Prevent self-intersections.
:returns: New simplified vector.
"""
ls = self.linestring
simplified = ls.simplify(tolerance, preserve_topology=preserve_topology)
newvec = self._make_vector(name=self.myname + '_simplified', fromshapely=simplified)
return newvec
# ----------------------------------------------------------------
# Deep copy
# ----------------------------------------------------------------
[docs]
def deepcopy_vector(self, name: str = None, parentzone=None) -> "vectorModel":
"""Return a deep copy of this vector.
:param name: Name for the copy (defaults to ``<name>_copy``).
:param parentzone: Optional parent zone for the copy.
"""
if name is None:
name = self.myname + "_copy"
if parentzone is not None:
copied_vector = self._make_vector(name=name, parentzone=parentzone)
else:
copied_vector = self._make_vector(name=name)
copied_vector.myvertices = copy.deepcopy(self.myvertices)
copied_vector.closed = self.closed
copied_vector.zdatum = self.zdatum
copied_vector.add_zdatum = self.add_zdatum
copied_vector.sdatum = self.sdatum
copied_vector.add_sdatum = self.add_sdatum
return copied_vector
[docs]
def deepcopy(self, name: str = None, parentzone=None) -> "vectorModel":
"""Alias for :pymeth:`deepcopy_vector`."""
return self.deepcopy_vector(name, parentzone)
# ----------------------------------------------------------------
# Interior / parts / subpolygons
# ----------------------------------------------------------------
@property
[docs]
def _parts(self) -> "zoneModel":
"""Split the vector into sub-parts based on in-use flags.
Vertices flagged as *not in use* act as separators between parts.
The first part is always closed.
"""
self.find_minmax()
parts = self._make_zone()
current_part = self._make_vector()
for curvert in self.myvertices:
if curvert.in_use:
current_part.add_vertex(curvert)
else:
if current_part.nbvertices > 0:
if parts.nbvectors > 0:
current_part.force_to_close()
parts.add_vector(current_part)
current_part = self._make_vector()
for curvert in current_part.myvertices:
parts.myvectors[0].add_vertex(curvert)
parts.myvectors[0].force_to_close()
return parts
[docs]
def get_subpolygons(self) -> list[list[wolfvertex]]:
"""Return triangulated sub-polygons for rendering, handling interiors.
:returns: List of vertex lists, each forming a triangle or polygon.
"""
if self.nbvertices == 0:
logging.debug(_('Vector {} has no vertices -- cannot create subpolygons').format(self.myname))
return []
if self.myprop.filled:
if self.has_interior:
parts = self._parts
tri = parts.create_constrainedDelaunay(nb=0)
centroid_interiors = [part.centroid for part in parts.myvectors[1:]]
tri.unuse_triangles_containing_points(centroid_interiors)
return tri.get_triangles_as_listwolfvertices()
else:
if self._simplified_geometry:
if self.myprop.closed and (self.myvertices[0].x != self.myvertices[-1].x or self.myvertices[0].y != self.myvertices[-1].y):
return [self.myvertices + [self.myvertices[0]]]
else:
return [self.myvertices]
else:
poly = self.polygon
if poly.is_empty or poly.area == 0:
logging.debug(_('Vector {} has no area -- cannot create triangulation').format(self.myname))
return []
# Validate polygon before calling the C-level triangle library.
# Invalid polygons (self-intersecting, bowtie, etc.) can cause
# the triangle C extension to abort()/exit() — crashing the
# entire process without a Python exception.
if not poly.is_valid or self.nbvertices < 4:
# Fallback: simple fan triangulation from vertex 0
verts = self.myvertices
n = len(verts)
fan_tris = []
for i in range(1, n - 1):
fan_tris.append([
wolfvertex(verts[0].x, verts[0].y, verts[0].z),
wolfvertex(verts[i].x, verts[i].y, verts[i].z),
wolfvertex(verts[i + 1].x, verts[i + 1].y, verts[i + 1].z),
])
return fan_tris
xx, yy = poly.exterior.xy
tr_x = np.array(xx).min()
tr_y = np.array(yy).min()
xx = np.array(xx) - tr_x
yy = np.array(yy) - tr_y
geom = {'vertices': [[x, y] for x, y in zip(xx[:-1], yy[:-1])],
'segments': [[i, i + 1] for i in range(len(xx) - 2)] + [[len(xx) - 2, 0]]}
try:
delaunay = tri_lib.triangulate(geom, 'p')
tri = []
for curtri in delaunay['triangles']:
tri.append([wolfvertex(delaunay['vertices'][curtri[i]][0] + tr_x,
delaunay['vertices'][curtri[i]][1] + tr_y) for i in range(3)])
return tri
except Exception:
# Fallback: fan triangulation
verts = self.myvertices
n = len(verts)
fan_tris = []
for i in range(1, n - 1):
fan_tris.append([
wolfvertex(verts[0].x, verts[0].y, verts[0].z),
wolfvertex(verts[i].x, verts[i].y, verts[i].z),
wolfvertex(verts[i + 1].x, verts[i + 1].y, verts[i + 1].z),
])
return fan_tris if fan_tris else []
else:
if self.has_interior:
alls = []
new_poly = []
alls.append(new_poly)
for curvert in self.myvertices:
if curvert.in_use:
new_poly.append(curvert)
else:
new_poly = []
alls.append(new_poly)
new_poly.append(curvert)
if self.myprop.closed and (self.myvertices[0].x != self.myvertices[-1].x or self.myvertices[0].y != self.myvertices[-1].y):
alls[0].append(self.myvertices[0])
return alls
else:
if self.myprop.closed and (self.myvertices[0].x != self.myvertices[-1].x or self.myvertices[0].y != self.myvertices[-1].y):
return [self.myvertices + [self.myvertices[0]]]
else:
return [self.myvertices]
@staticmethod
[docs]
def _fillet_polyline(xs: np.ndarray, ys: np.ndarray, ws: np.ndarray,
join_size: float, join_size_mode: int,
is_closed: bool, n_subdiv: int = 8
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Replace sharp corners with quadratic Bézier fillets.
For each interior corner, the original vertex is replaced by a
smooth curve starting/ending at a controlled distance along the
adjacent segments.
:param xs, ys: coordinate arrays (float64).
:param ws: width multiplier array (float32).
:param join_size: Fraction (mode 0) or world distance (mode 1).
In fraction mode, 0 = sharp corner, 1 = fillet eats half of
each adjacent segment (clamped so adjacent fillets don't
overlap).
:param join_size_mode: 0 = fraction (adimensional), 1 = world.
:param is_closed: whether the polyline is closed.
:param n_subdiv: number of Bézier subdivision points per corner.
:returns: (new_xs, new_ys, new_ws) with inserted fillet points.
"""
n = len(xs)
if n < 3 or join_size <= 0.0:
return xs, ys, ws
# Compute segment lengths
seg_lens = np.sqrt(np.diff(xs)**2 + np.diff(ys)**2) # (n-1,)
# For each interior vertex, compute cut distance along each
# adjacent segment. Interior vertices are indices 1..n-2
# (plus index 0 and n-1 if closed).
# cut_before[i] = distance from vertex i back toward vertex i-1
# cut_after[i] = distance from vertex i forward toward vertex i+1
cut_before = np.zeros(n)
cut_after = np.zeros(n)
# Determine which vertices get a fillet
if is_closed:
corners = list(range(n - 1)) # last == first, skip it
else:
corners = list(range(1, n - 1))
for i in corners:
seg_before = seg_lens[i - 1] if i > 0 else seg_lens[-1] if is_closed else 0.0
seg_after = seg_lens[i] if i < n - 1 else seg_lens[0] if is_closed else 0.0
if seg_before < 1e-12 or seg_after < 1e-12:
continue
if join_size_mode == 0:
# Fraction mode: join_size * half_segment_length
d_before = join_size * 0.5 * seg_before
d_after = join_size * 0.5 * seg_after
else:
# World distance mode
d_before = min(join_size, 0.5 * seg_before)
d_after = min(join_size, 0.5 * seg_after)
cut_before[i] = d_before
cut_after[i] = d_after
# Clamp so adjacent fillets don't overlap on the same segment.
# For segment i (from vertex i to i+1): cut_after[i] + cut_before[i+1] <= seg_lens[i]
for seg_i in range(len(seg_lens)):
i_start = seg_i
i_end = (seg_i + 1) % n if is_closed else seg_i + 1
total = cut_after[i_start] + cut_before[i_end]
if total > seg_lens[seg_i] and total > 1e-12:
scale = seg_lens[seg_i] / total
cut_after[i_start] *= scale
cut_before[i_end] *= scale
# Build output arrays
out_x, out_y, out_w = [], [], []
for i in range(n):
is_corner = i in corners if not is_closed else (i < n - 1)
if not is_corner or (cut_before[i] < 1e-12 and cut_after[i] < 1e-12):
out_x.append(xs[i])
out_y.append(ys[i])
out_w.append(ws[i])
continue
# Previous and next vertex indices
i_prev = (i - 1) % (n - 1) if is_closed else i - 1
i_next = (i + 1) % (n - 1) if is_closed else i + 1
# Bézier control points:
# P0 = point on segment [i_prev, i] at distance cut_before from i
# P1 = original corner vertex (control point)
# P2 = point on segment [i, i_next] at distance cut_after from i
d_b = cut_before[i]
d_a = cut_after[i]
seg_b = seg_lens[i - 1] if i > 0 else seg_lens[-1]
seg_a = seg_lens[i] if i < n - 1 else seg_lens[0]
t_b = d_b / seg_b if seg_b > 1e-12 else 0.0
t_a = d_a / seg_a if seg_a > 1e-12 else 0.0
# P0: lerp from corner toward previous vertex
p0x = xs[i] + t_b * (xs[i_prev] - xs[i])
p0y = ys[i] + t_b * (ys[i_prev] - ys[i])
w0 = ws[i] + t_b * (ws[i_prev] - ws[i])
# P1: corner
p1x, p1y = xs[i], ys[i]
w1 = ws[i]
# P2: lerp from corner toward next vertex
p2x = xs[i] + t_a * (xs[i_next] - xs[i])
p2y = ys[i] + t_a * (ys[i_next] - ys[i])
w2 = ws[i] + t_a * (ws[i_next] - ws[i])
# Emit n_subdiv+1 points along quadratic Bézier (P0, P1, P2)
for k in range(n_subdiv + 1):
t = k / n_subdiv
# De Casteljau / quadratic formula: B(t) = (1-t)²P0 + 2(1-t)t·P1 + t²P2
omt = 1.0 - t
bx = omt * omt * p0x + 2.0 * omt * t * p1x + t * t * p2x
by = omt * omt * p0y + 2.0 * omt * t * p1y + t * t * p2y
bw = omt * omt * w0 + 2.0 * omt * t * w1 + t * t * w2
out_x.append(bx)
out_y.append(by)
out_w.append(bw)
return (np.array(out_x, dtype=np.float64),
np.array(out_y, dtype=np.float64),
np.array(out_w, dtype=np.float32))
@staticmethod
[docs]
def _expand_width_profile(profile: list, n: int) -> np.ndarray:
"""Expand a raw width profile to exactly *n* per-vertex multipliers.
Supports two storage formats:
* **Dense** – a plain ``list[float]``, one value per vertex.
If the length differs from *n* (e.g. a phantom closing vertex was
appended by :meth:`get_subpolygons`), the values are resampled with
linear interpolation so the profile always spans the full polyline.
* **Sparse** – a ``list`` of ``(int, float)`` tuples
``[(index, multiplier), ...]``.
Unspecified vertices are filled by linear interpolation between
adjacent control points; extrapolation uses the nearest endpoint
value (constant fill). A single control point produces a uniform
profile at that multiplier value.
:param profile: Dense or sparse width profile (as returned by
:meth:`~vectorproperties._parse_width_profile`).
:param n: Number of vertices the result must cover.
:returns: float32 array of length *n*, all values ≥ 0.
"""
if not profile or n == 0:
return np.ones(n, dtype=np.float32)
# Detect sparse format: list contains tuples/lists
if isinstance(profile[0], (tuple, list)):
pairs = sorted(profile, key=lambda p: p[0])
src_idx = np.array([p[0] for p in pairs], dtype=np.float64)
src_val = np.array([p[1] for p in pairs], dtype=np.float64)
if len(pairs) == 1:
return np.full(n, src_val[0], dtype=np.float32)
dst_idx = np.arange(n, dtype=np.float64)
result = np.interp(dst_idx, src_idx, src_val,
left=src_val[0], right=src_val[-1])
return result.astype(np.float32)
else:
# Dense format
m = len(profile)
if m == n:
return np.asarray(profile, dtype=np.float32)
if m == 0:
return np.ones(n, dtype=np.float32)
# Resample: stretch/compress the `m`-value profile over `n` slots
src_idx = np.linspace(0, n - 1, m)
dst_idx = np.arange(n, dtype=np.float64)
result = np.interp(dst_idx, src_idx,
np.asarray(profile, dtype=np.float64))
return result.astype(np.float32)
[docs]
def build_shader_vbo_data(self) -> tuple[np.ndarray, int]:
"""Build VBO data for shader rendering of this vector.
Each sub-polyline (from :meth:`get_subpolygons` in non-filled mode)
is converted to ``GL_LINES_ADJACENCY`` segments: four vertices per
segment (prev, start, end, next), each carrying
``(x, y, curvilinear_distance, width_multiplier)``.
:returns: ``(vbo_data, vertex_count)`` — contiguous float32 array
and total number of vertices emitted.
"""
if self.nbvertices < 2:
return np.array([], dtype=np.float32), 0
all_polys = self.get_subpolygons()
parts = []
total = 0
widths_profile = getattr(self.myprop, 'width_profile', None)
for poly in all_polys:
n = len(poly)
if n < 2:
continue
xs = np.array([v.x for v in poly], dtype=np.float64)
ys = np.array([v.y for v in poly], dtype=np.float64)
diffs = np.sqrt(np.diff(xs)**2 + np.diff(ys)**2)
dists = np.zeros(n, dtype=np.float32)
dists[1:] = np.cumsum(diffs).astype(np.float32)
if widths_profile is not None and len(widths_profile) > 0:
ws = self._expand_width_profile(widths_profile, n)
else:
ws = np.ones(n, dtype=np.float32)
# Detect closed polyline (first vertex ≈ last vertex)
is_closed = (n >= 3 and
abs(xs[0] - xs[-1]) < 1e-10 and
abs(ys[0] - ys[-1]) < 1e-10)
# Apply CPU-side fillet subdivision when join_style == 3
join_style = getattr(self.myprop, 'join_style', 0)
if join_style == 3 and n >= 3:
join_size = getattr(self.myprop, 'join_size', 0.5)
join_size_mode = getattr(self.myprop, 'join_size_mode', 0)
xs, ys, ws = self._fillet_polyline(
xs, ys, ws, join_size, join_size_mode, is_closed)
n = len(xs)
if n < 2:
continue
# Recompute cumulative distances after subdivision
diffs = np.sqrt(np.diff(xs)**2 + np.diff(ys)**2)
dists = np.zeros(n, dtype=np.float32)
dists[1:] = np.cumsum(diffs).astype(np.float32)
# Emit 4 vertices per segment for GL_LINES_ADJACENCY:
# [prev, start, end, next]
num_seg = n - 1
seg_data = np.empty((num_seg * 4, 4), dtype=np.float32)
for i in range(num_seg):
# Previous vertex (adjacency)
if i > 0:
pi = i - 1
elif is_closed:
pi = n - 2 # wrap to last real vertex
else:
pi = i # degenerate: same as start
# Next vertex (adjacency)
if i < num_seg - 1:
ni = i + 2
elif is_closed:
ni = 1 # wrap to first real vertex after start
else:
ni = i + 1 # degenerate: same as end
seg_data[4*i] = [xs[pi], ys[pi], dists[pi], ws[pi]]
seg_data[4*i + 1] = [xs[i], ys[i], dists[i], ws[i]]
seg_data[4*i + 2] = [xs[i+1], ys[i+1], dists[i+1], ws[i+1]]
seg_data[4*i + 3] = [xs[ni], ys[ni], dists[ni], ws[ni]]
parts.append(seg_data)
total += num_seg * 4
if total == 0:
return np.array([], dtype=np.float32), 0
return np.ascontiguousarray(np.concatenate(parts)), total
# ----------------------------------------------------------------
# Linked-array extraction
# ----------------------------------------------------------------
[docs]
def get_values_linked_polygon(self, linked_arrays: list, getxy=False) -> dict:
"""Extract values from arrays inside this polygon.
:param linked_arrays: List of arrays to query.
:param getxy: Also return XY coordinates.
:returns: Dict mapping array index to extracted values.
"""
vals = {}
for curarray in linked_arrays:
if curarray.plotted:
vals[curarray.idx] = curarray.get_values_insidepoly(self, getxy=getxy)
else:
vals[curarray.idx] = None
return vals
[docs]
def get_all_values_linked_polygon(self, linked_arrays, getxy=False) -> dict:
"""Extract all values from arrays inside this polygon.
:param linked_arrays: Arrays to query.
:param getxy: Also return XY coordinates.
:returns: Dict mapping array index to all extracted values.
"""
vals = {}
for curarray in linked_arrays:
if curarray.plotted:
vals[curarray.idx] = curarray.get_all_values_insidepoly(self, getxy=getxy)
else:
vals[curarray.idx] = None
return vals
[docs]
def get_all_values_linked_polyline(self, linked_arrays, getxy=True) -> dict:
"""Extract values from arrays under this polyline.
:param linked_arrays: Arrays to query.
:param getxy: Also return XY coordinates.
:returns: Dict mapping array index to extracted values.
"""
vals = {}
for curarray in linked_arrays:
if curarray.plotted:
vals[curarray.idx], xy = curarray.get_all_values_underpoly(self, getxy=getxy)
else:
vals[curarray.idx] = None
return vals
[docs]
def get_values_on_vertices(self, curarray):
"""Sample a raster array at each vertex position and store as Z.
:param curarray: Array from which to sample values.
"""
if not curarray.plotted:
logging.warning(_('Array {} is not plotted, cannot get values on vertices').format(curarray.idx))
return
for curpt in self.myvertices:
curpt.z = curarray.get_value(curpt.x, curpt.y)
[docs]
def get_values_linked(self, linked_arrays: dict, refine=True, filter_null=False):
"""Extract values from linked arrays along the polyline.
:param linked_arrays: Dict of ``{label: array}`` to query.
:param refine: Refine the polyline to the array resolution before sampling.
:param filter_null: Exclude null values (``-99999``) from results.
:returns: A ``zoneModel`` containing one vector per array.
"""
exit_flag = True
for curlabel, curarray in linked_arrays.items():
if curarray.plotted:
exit_flag = False
if exit_flag:
return
if refine:
myzone = self._make_zone(name='linked_arrays - fine step')
for curlabel, curarray in linked_arrays.items():
if curarray.plotted:
myvec = self._make_vector(name=curlabel, parentzone=myzone)
myzone.add_vector(myvec)
ds = curarray.get_dxdy_min()
pts = self._refine2D(ds)
allz = [curarray.get_value(curpt.x, curpt.y, nullvalue=-99999) for curpt in pts]
if filter_null:
for curpt, curz in zip(pts, allz):
if curz != -99999:
myvec.add_vertex(wolfvertex(curpt.x, curpt.y, curz))
else:
for curpt, curz in zip(pts, allz):
myvec.add_vertex(wolfvertex(curpt.x, curpt.y, curz))
else:
myzone = self._make_zone(name='linked_arrays')
for curlabel, curarray in linked_arrays.items():
if curarray.plotted:
myvec = self._make_vector(name=curlabel, parentzone=myzone)
myzone.add_vector(myvec)
if filter_null:
for curpt in self.myvertices:
locval = curarray.get_value(curpt.x, curpt.y, nullvalue=-99999)
if locval != -99999:
myvec.add_vertex(wolfvertex(curpt.x, curpt.y, locval))
else:
for curpt in self.myvertices:
locval = curarray.get_value(curpt.x, curpt.y, nullvalue=-99999)
myvec.add_vertex(wolfvertex(curpt.x, curpt.y, locval))
return myzone
# ----------------------------------------------------------------
# Miscellaneous
# ----------------------------------------------------------------
[docs]
def interpolate_coordinates(self):
"""Linearly interpolate Z values at vertices where Z is invalid ``(-99999, 99999, NaN)``."""
sz = self.get_sz()
s = sz[0]
z = sz[1]
valid_indices = np.where((z != -99999.) & (z != 99999.) & (z != '') & (np.isfinite(z)))[0]
if len(valid_indices) == 0:
logging.warning(_('No valid z values to interpolate'))
return
f = interp1d(s[valid_indices], z[valid_indices])
for k in range(self.nbvertices):
if k not in valid_indices:
zval = f(s[k])
self.myvertices[k].z = zval
self.update_lengths()
self._on_vertices_changed()
[docs]
def create_cs(self, wa) -> "vectorModel":
"""Create a cross-section profile by sampling a raster array.
:param wa: ``WolfArray`` to sample.
:returns: A refined vector with Z values from the array, or ``None``.
"""
if not self.used:
return None
if wa is None:
logging.warning(_('No array - Cannot create cross-sections !'))
return
ds = wa.get_dxdy_min()
new_vec = self.deepcopy()
new_vec.split(ds=ds, new=False)
new_vec.get_values_on_vertices(wa)
return new_vec
[docs]
def keep_only_trace(self):
"""Reduce the vector to only two vertices: first and last, with Z = 0."""
self.myvertices = [self.myvertices[0].copy(), self.myvertices[-1].copy()]
self.myvertices[0].z = 0.0
self.myvertices[1].z = 0.0
self.update_lengths()
[docs]
def extend_along_trace(self, distance: float | str):
"""Extend the trace (first/last vertices) by *distance* in both directions.
:param distance: Extension distance in metres, or a percentage string (e.g. ``'10%'``).
"""
if self.nbvertices < 2:
logging.warning(_('Vector has less than 2 vertices, cannot extend along trace'))
return
if isinstance(distance, str):
if distance.endswith('%'):
try:
perc = float(distance[:-1]) / 100.0
total_length = self.length2D
distance = perc * total_length
except ValueError:
logging.warning(_('Invalid distance format'))
return
else:
logging.warning(_('Invalid distance format'))
return
start_vec = self.myvertices[-1] - self.myvertices[0]
start_vec.x /= self.length2D
start_vec.y /= self.length2D
self.myvertices[0].x -= start_vec.x * distance
self.myvertices[0].y -= start_vec.y * distance
self.myvertices[-1].x += start_vec.x * distance
self.myvertices[-1].y += start_vec.y * distance
# ====================================================================
# zoneModel
# ====================================================================
[docs]
class zoneModel:
"""Pure-data zone: a list of vectorModel instances plus geometry helpers.
GUI subclass ``zone`` in :mod:`wolfhece.PyVertexvectors` adds
plotting, OpenGL display-list support and wx integration.
"""
[docs]
myvectors: list[vectorModel]
def __init__(self,
lines: list[str] = [],
name: str = 'NoName',
parent=None,
is2D: bool = True,
fromshapely: Union[LineString, Polygon, MultiLineString, MultiPolygon] = None) -> None:
self.myname = ''
[docs]
self.active_vector = None
self.xmin = -99999.
self.ymin = -99999.
self.xmax = -99999.
self.ymax = -99999.
[docs]
self.has_legend = False
[docs]
self._move_start = None
[docs]
self._rotation_center = None
[docs]
self._rotation_step = None
if len(lines) > 0:
self.myname = lines[0]
tmp_nbvectors = int(lines[1])
self.myvectors = []
curstart = 2
if tmp_nbvectors > 1000:
logging.info(_('Many vectors in zone -- {} -- Be patient !').format(tmp_nbvectors))
for i in range(tmp_nbvectors):
curvec = self._make_vector(lines=lines[curstart:], parentzone=self, is2D=is2D)
curstart += curvec._nblines()
self.myvectors.append(curvec)
if tmp_nbvectors > 1000:
if i % 100 == 0:
logging.info(_('{} vectors read').format(i))
if name != '' and self.myname == '':
self.myname = name
self.myvectors = []
[docs]
self.selected_vectors = []
[docs]
self.multils: MultiLineString = None
if fromshapely is not None:
self.import_shapelyobj(fromshapely)
# ----------------------------------------------------------------
# Factory methods (overridden in GUI subclass)
# ----------------------------------------------------------------
[docs]
def _make_vector(self, **kwargs) -> vectorModel:
"""Create a sibling vector. GUI subclass returns ``vector``."""
return vectorModel(**kwargs)
[docs]
def _make_zone(self, **kwargs) -> "zoneModel":
"""Create a sibling zone. GUI subclass returns ``zone``."""
return zoneModel(**kwargs)
[docs]
def _make_triangulation(self, **kwargs) -> TriangulationModel:
"""Create a triangulation. GUI subclass returns ``Triangulation``."""
return TriangulationModel(**kwargs)
# ----------------------------------------------------------------
# GUI hooks (no-op in model)
# ----------------------------------------------------------------
[docs]
def reset_listogl(self):
"""No-op in model. GUI subclass deletes gl list and sets idgllist=-99999."""
pass
[docs]
def _fill_structure(self):
"""No-op in model. GUI subclass calls parent.fill_structure()."""
pass
# ----------------------------------------------------------------
# Properties
# ----------------------------------------------------------------
@property
[docs]
def nbvectors(self) -> int:
"""Number of vectors in the zone."""
return len(self.myvectors)
@property
[docs]
def area(self) -> float:
"""Total area of all vectors in the zone."""
return sum(self.areas)
@property
[docs]
def areas(self) -> list[float]:
"""List of individual vector areas."""
return [curvec.surface for curvec in self.myvectors]
@property
[docs]
def vector_names(self) -> list[str]:
"""List of vector names in the zone."""
return [cur.myname for cur in self.myvectors]
# ----------------------------------------------------------------
# Vector management
# ----------------------------------------------------------------
[docs]
def add_vector(self, addedvect: vectorModel, index=-99999, forceparent=False, update_struct=False):
"""Add a vector to the zone.
:param addedvect: Vector to add.
:param index: Insertion position (``-99999`` or > count appends).
:param forceparent: Set this zone as the vector's parent.
:param update_struct: Trigger GUI structure update.
"""
if index == -99999 or index > self.nbvectors:
self.myvectors.append(addedvect)
else:
self.myvectors.insert(index, addedvect)
if forceparent:
addedvect.parentzone = self
if self.nbvectors == 1:
self.active_vector = addedvect
if update_struct:
self._fill_structure()
[docs]
def get_vector(self, keyvector: Union[int, str]) -> vectorModel:
"""Return a vector by index or name.
:param keyvector: Integer index or string name.
:returns: The matching ``vectorModel``, or ``None``.
"""
if isinstance(keyvector, int):
if keyvector < self.nbvectors:
return self.myvectors[keyvector]
return None
if isinstance(keyvector, str):
zone_names = [cur.myname for cur in self.myvectors]
if keyvector in zone_names:
return self.myvectors[zone_names.index(keyvector)]
return None
def __getitem__(self, ndx: Union[int, str]) -> vectorModel:
"""Shortcut for :pymeth:`get_vector`."""
return self.get_vector(ndx)
[docs]
def count(self):
"""No-op placeholder (kept for backward compatibility)."""
return
[docs]
def import_shapelyobj(self, obj):
"""Import a Shapely geometry (LineString, Polygon, Multi…) as vector(s).
:param obj: Shapely geometry object.
"""
if isinstance(obj, LineString):
curvec = self._make_vector(fromshapely=obj, parentzone=self, name=self.myname)
self.add_vector(curvec)
elif isinstance(obj, Polygon):
curvec = self._make_vector(fromshapely=obj, parentzone=self, name=self.myname)
self.add_vector(curvec)
elif isinstance(obj, MultiLineString):
for curls in list(obj.geoms):
curvec = self._make_vector(fromshapely=curls, parentzone=self, name=self.myname)
self.add_vector(curvec)
elif isinstance(obj, MultiPolygon):
for curpoly in list(obj.geoms):
curvec = self._make_vector(fromshapely=curpoly, parentzone=self, name=self.myname)
self.add_vector(curvec)
else:
logging.warning(_('Object type {} not supported -- Update "import_shapelyobj"').format(type(obj)))
# ----------------------------------------------------------------
# Colour / legend / property helpers
# ----------------------------------------------------------------
[docs]
def check_if_interior_exists(self):
"""Check all vectors for interior boundaries."""
list(map(lambda curvec: curvec.check_if_interior_exists(), self.myvectors))
[docs]
def add_values(self, key: str, values: np.ndarray):
"""Assign a value to each vector.
:param key: Property key.
:param values: 1-D array with one value per vector.
"""
if self.nbvectors != values.shape[0]:
logging.warning(_('Number of vectors and values do not match'))
return
list(map(lambda cur: cur[0].add_value(key, cur[1]), zip(self.myvectors, values)))
[docs]
def get_values(self, key: str) -> np.ndarray:
"""Retrieve values for a given key from all vectors.
:param key: Property key.
"""
return np.array([curvec.get_value(key) for curvec in self.myvectors])
[docs]
def set_colors_from_value(self, key: str, cmap: wolfpalette | Colormap | cm.ScalarMappable,
vmin: float = 0., vmax: float = 1.):
"""Colour all vectors according to a stored value and a colormap.
:param key: Property key.
:param cmap: Colour map.
:param vmin: Minimum value for mapping.
:param vmax: Maximum value for mapping.
"""
list(map(lambda curvec: curvec.set_color_from_value(key, cmap, vmin, vmax), self.myvectors))
[docs]
def set_color(self, color: int | RGB | str | list | tuple):
"""Set the same colour on all vectors.
:param color: Colour specification.
"""
list(map(lambda curvec: curvec.set_color(color), self.myvectors))
[docs]
def set_linewidth(self, linewidth: int | float):
"""Set the line width for all vectors.
:param linewidth: Line width in pixels.
"""
list(map(lambda curvec: curvec.set_linewidth(linewidth), self.myvectors))
[docs]
def set_alpha(self, alpha: int):
"""Set the transparency for all vectors.
:param alpha: Alpha value (0–255).
"""
list(map(lambda curvec: curvec.set_alpha(alpha), self.myvectors))
[docs]
def set_filled(self, filled: bool):
"""Enable or disable polygon filling for all vectors.
:param filled: ``True`` to fill.
"""
list(map(lambda curvec: curvec.set_filled(filled), self.myvectors))
[docs]
def check_if_open(self):
"""Check all vectors for open/closed state."""
list(map(lambda curvect: curvect.check_if_open(), self.myvectors))
[docs]
def set_legend_text(self, text: str):
"""Set legend text for all vectors.
:param text: Legend text.
"""
list(map(lambda curvect: curvect.set_legend_text(text), self.myvectors))
[docs]
def set_legend_color(self, color: int | RGB | str | list | tuple):
"""Set legend colour for all vectors.
:param color: Colour specification.
"""
list(map(lambda curvect: curvect.set_legend_color(color), self.myvectors))
[docs]
def set_legend_text_from_values(self, key: str):
"""Set legend text from stored values for all vectors.
:param key: Property key.
"""
list(map(lambda curvect: curvect.set_legend_text_from_value(key), self.myvectors))
[docs]
def set_legend_position(self, x, y):
"""Set legend position for all vectors.
:param x: X position (float or keyword string).
:param y: Y position (float or keyword string).
"""
list(map(lambda curvect: curvect.set_legend_position(x, y), self.myvectors))
[docs]
def set_legend_to_centroid(self):
"""Place legends at each vector's centroid."""
list(map(lambda curvec: curvec.set_legend_to_centroid(), self.myvectors))
[docs]
def set_legend_visible(self, visible: bool = True):
"""Show or hide legends for all vectors.
:param visible: ``True`` to display.
"""
list(map(lambda curvec: curvec.set_legend_visible(visible), self.myvectors))
[docs]
def keep_only_trace(self):
"""Reduce all vectors to their trace (first and last vertices)."""
list(map(lambda curvec: curvec.keep_only_trace(), self.myvectors))
[docs]
def extend_along_trace(self, distance: float | str):
"""Extend all vectors along their trace.
:param distance: Extension distance or percentage string.
"""
list(map(lambda curvec: curvec.extend_along_trace(distance), self.myvectors))
# ----------------------------------------------------------------
# Shader-property helpers (propagate to all vectors)
# ----------------------------------------------------------------
[docs]
def set_glow(self, enabled: bool = True, width: float = 0.4,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.4):
"""Configure line glow for all vectors.
:param enabled: Enable or disable glow.
:param width: Glow width as fraction of half-width (0–1).
:param color: RGB colour tuple [0, 1].
:param alpha: Glow opacity [0, 1].
"""
list(map(lambda v: v.set_glow(enabled, width, color, alpha), self.myvectors))
[docs]
def set_dash(self, enabled: bool = True, dash_length: float = 10.0,
gap_length: float = 5.0):
"""Configure dashed lines for all vectors.
:param enabled: Enable or disable dashes.
:param dash_length: Dash length in world units.
:param gap_length: Gap length in world units.
"""
list(map(lambda v: v.set_dash(enabled, dash_length, gap_length), self.myvectors))
[docs]
def set_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure line animation for all vectors.
:param mode: 0=none, 1=pulse, 2=marching ants.
:param speed: Speed multiplier.
"""
list(map(lambda v: v.set_animation(mode, speed), self.myvectors))
[docs]
def set_fill_animation(self, mode: int = 1, speed: float = 1.0,
center_index: int = 0,
start_angle: float = 90.0):
"""Configure filled-polygon animation for all vectors.
:param mode: 0=none, 1=pulse, 2=sweep, 3=ripple,
4=radial progressive, 5=clockwise clock, 6=counter-clockwise clock.
:param speed: Speed multiplier.
:param center_index: Vertex index used as the radial/clock centre.
:param start_angle: Clock fill start angle in degrees.
"""
list(map(lambda v: v.set_fill_animation(mode, speed, center_index, start_angle), self.myvectors))
[docs]
def set_fill_color(self, color=None):
"""Set independent fill colour for all vectors.
:param color: RGB tuple, packed int, or None to inherit main colour.
"""
list(map(lambda v: v.set_fill_color(color), self.myvectors))
[docs]
def set_contour_color(self, color=None, width: float = 1.0,
enabled: bool = True):
"""Set contour colour/width for all vectors.
:param color: RGB tuple, packed int, or None to inherit main colour.
:param width: Contour line width.
:param enabled: Whether to draw the contour pass.
"""
list(map(lambda v: v.set_contour_color(color, width, enabled), self.myvectors))
[docs]
def set_join_style(self, style: int = 0, size: float = 0.5,
size_mode: int = 0):
"""Configure join rendering for all vectors.
:param style: 0=miter, 1=bevel, 2=round, 3=fillet.
:param size: Join size.
:param size_mode: 0=fraction, 1=world distance.
"""
list(map(lambda v: v.set_join_style(style, size, size_mode), self.myvectors))
[docs]
def set_width_in_pixels(self, pixels: bool = True):
"""Toggle pixel vs world-unit line width for all vectors.
:param pixels: ``True`` for pixel widths.
"""
list(map(lambda v: v.set_width_in_pixels(pixels), self.myvectors))
[docs]
def set_text_along(self, text: str, enabled: bool = True,
offset: float = 0.0, perp: float = 0.0,
alignment: str = 'center'):
"""Configure text along polyline for all vectors.
:param text: Text to display.
:param enabled: Enable or disable.
:param offset: Shift start along polyline (world units).
:param perp: Perpendicular offset (world units).
:param alignment: 'left', 'center', 'right'.
"""
list(map(lambda v: v.set_text_along(text, enabled, offset, perp, alignment), self.myvectors))
[docs]
def set_legend_glow(self, enabled: bool = True, width: float = 0.15,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.5):
"""Configure legend glow for all vectors.
:param enabled: Enable or disable legend glow.
:param width: SDF threshold offset.
:param color: RGB colour tuple [0, 1].
:param alpha: Glow opacity [0, 1].
"""
list(map(lambda v: v.set_legend_glow(enabled, width, color, alpha), self.myvectors))
[docs]
def set_legend_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure legend animation for all vectors.
:param mode: 0=none, 1=pulse, 2=wave, 3=typewriter.
:param speed: Speed multiplier.
"""
list(map(lambda v: v.set_legend_animation(mode, speed), self.myvectors))
[docs]
def set_legend_style(self, smoothing: float = 1.0,
alignment: str = 'left',
line_spacing: float = 1.2):
"""Configure legend SDF text style for all vectors.
:param smoothing: SDF AA multiplier.
:param alignment: 'left', 'center', 'right'.
:param line_spacing: Multiline spacing multiplier.
"""
list(map(lambda v: v.set_legend_style(smoothing, alignment, line_spacing), self.myvectors))
# ----------------------------------------------------------------
# Cache / move / rotate
# ----------------------------------------------------------------
[docs]
def set_cache(self):
"""Cache coordinates and bounds for all vectors."""
list(map(lambda curvect: curvect.set_cache(), self.myvectors))
[docs]
def clear_cache(self):
"""Clear cached data for all vectors and reset transform state."""
list(map(lambda curvect: curvect.clear_cache(), self.myvectors))
self._move_start = None
self._move_step = None
self._rotation_center = None
self._rotation_step = None
self.find_minmax(update=True)
[docs]
def move(self, dx: float, dy: float, use_cache: bool = True, inplace: bool = True):
"""Translate all vectors in the zone.
:param dx: Displacement along X.
:param dy: Displacement along Y.
:param use_cache: Use cached coordinates.
:param inplace: Modify in place; if ``False`` return a moved copy.
"""
if self._move_step is not None:
dx = np.round(dx / self._move_step) * self._move_step
dy = np.round(dy / self._move_step) * self._move_step
if inplace:
for curvect in self.myvectors:
curvect.move(dx, dy, use_cache)
return self
else:
newzone = self.deepcopy()
newzone.move(dx, dy, use_cache=False)
return newzone
[docs]
def rotate(self, angle: float, center: tuple[float, float], use_cache: bool = True, inplace: bool = True):
"""Rotate all vectors around a centre point.
:param angle: Rotation angle in degrees.
:param center: Centre of rotation ``(x, y)``.
:param use_cache: Use cached coordinates.
:param inplace: Modify in place or return a copy.
"""
if inplace:
for curvect in self.myvectors:
curvect.rotate(angle, center, use_cache)
return self
else:
newzone = self.deepcopy()
newzone.rotate(angle, center, use_cache=False)
return newzone
[docs]
def rotate_xy(self, x: float, y: float, use_cache: bool = True, inplace: bool = True):
if self._rotation_center is None:
logging.error(_('No rotation center defined - Set it before rotating by this routine'))
return self
angle = np.degrees(-np.arctan2(y - self._rotation_center[1], x - self._rotation_center[0]))
if self._rotation_step is not None:
angle = np.round(angle / self._rotation_step) * self._rotation_step
return self.rotate(angle, self._rotation_center, use_cache, inplace)
[docs]
def use(self):
"""Mark all vectors and the zone as used/active."""
for curvect in self.myvectors:
curvect.used = True
self.used = True
self.reset_listogl()
[docs]
def unuse(self):
"""Mark all vectors and the zone as unused/inactive."""
for curvect in self.myvectors:
curvect.used = False
self.used = False
self.reset_listogl()
# ----------------------------------------------------------------
# Spatial queries / bounds
# ----------------------------------------------------------------
[docs]
def get_mapviewer(self):
"""Return the parent map-viewer widget."""
return self.parent.get_mapviewer()
[docs]
def find_minmax(self, update=False, only_firstlast: bool = False):
"""Compute the bounding box across all vectors.
:param update: Recompute individual vector bounds first.
:param only_firstlast: Consider only first/last vertices.
"""
if update:
for vect in self.myvectors:
vect.find_minmax(only_firstlast=only_firstlast)
if self.nbvectors == 0:
self.xmin = -99999.
self.ymin = -99999.
self.xmax = -99999.
self.ymax = -99999.
else:
minsx = np.asarray([vect.xmin for vect in self.myvectors if vect.xmin != -99999.])
minsy = np.asarray([vect.ymin for vect in self.myvectors if vect.ymin != -99999.])
maxsx = np.asarray([vect.xmax for vect in self.myvectors if vect.xmax != -99999.])
maxsy = np.asarray([vect.ymax for vect in self.myvectors if vect.ymax != -99999.])
if minsx.size == 0:
self.xmin = -99999.
self.ymin = -99999.
self.xmax = -99999.
self.ymax = -99999.
else:
self.xmin = float(minsx.min())
self.xmax = float(maxsx.max())
self.ymin = float(minsy.min())
self.ymax = float(maxsy.max())
[docs]
def find_nearest_vertex(self, x: float, y: float) -> wolfvertex | None:
"""Return the nearest vertex across all vectors in the zone."""
distmin = np.inf
minvert = None
for curvect in self.myvectors:
curvert = curvect.find_nearest_vertex(x, y)
if curvert is None:
continue
dist = (curvert.x - x) ** 2. + (curvert.y - y) ** 2.
if dist < distmin:
distmin = dist
minvert = curvert
return minvert
[docs]
def find_nearest_vector(self, x: float, y: float) -> vectorModel | None:
"""Return the vector whose geometry is nearest to ``(x, y)``."""
distmin = np.inf
minvec = None
for curvect in self.myvectors:
dist = curvect.distance_to_geometry(x, y)
if dist < distmin:
distmin = dist
minvec = curvect
return minvec
[docs]
def select_vectors_from_point(self, x: float, y: float, inside=True):
"""Select vectors containing or nearest to the point ``(x, y)``.
:param x: X coordinate.
:param y: Y coordinate.
:param inside: If ``True``, select vectors whose polygon contains the point;
otherwise select the nearest vector.
"""
if self.nbvectors == 0:
logging.warning(_('No vector in zone -- {}').format(self.myname))
return
self.selected_vectors.clear()
if inside:
for curvect in self.myvectors:
if curvect.isinside(x, y):
self.selected_vectors.append((curvect, 99999.))
else:
vectmin = self.find_nearest_vector(x, y)
if vectmin is not None:
self.selected_vectors.append((vectmin, vectmin.distance_to_geometry(x, y)))
[docs]
def get_selected_vectors(self, all=False):
if all:
mylist = []
if len(self.selected_vectors) > 0:
mylist.append(self.selected_vectors)
return mylist
else:
if len(self.selected_vectors) > 0:
return self.selected_vectors[0]
return None
[docs]
def asshapely_ls(self) -> MultiLineString:
"""Return all vectors as a Shapely ``MultiLineString``."""
mylines = []
for curvect in self.myvectors:
mylines.append(curvect.asshapely_ls())
return MultiLineString(mylines)
[docs]
def prepare_shapely(self):
"""Build a Shapely ``MultiLineString`` for the zone."""
self.multils = self.asshapely_ls()
# ----------------------------------------------------------------
# I/O
# ----------------------------------------------------------------
[docs]
def _nblines(self) -> int:
"""Return the number of lines required when saving this zone."""
nb = 2
for curvec in self.myvectors:
nb += curvec._nblines()
return nb
[docs]
def save(self, f: io.TextIOWrapper):
"""Write the zone definition to an open file.
:param f: Writable text file handle.
"""
f.write(self.myname + '\n')
f.write(str(self.nbvectors) + '\n')
for curvect in self.myvectors:
curvect.save(f)
[docs]
def export_shape(self, fn: str = ''):
"""Export vectors as an ESRI Shapefile (Lambert 72 / EPSG:31370).
:param fn: Output file path.
"""
from osgeo import osr, ogr
fn = str(fn)
srs = osr.SpatialReference()
srs.ImportFromEPSG(31370)
driver = ogr.GetDriverByName("ESRI Shapefile")
if not fn.endswith('.shp'):
fn += '.shp'
ds = driver.CreateDataSource(fn)
layer = ds.CreateLayer("poly", srs, ogr.wkbPolygon)
idFields = []
idFields.append(ogr.FieldDefn('curvi', ogr.OFTReal))
layer.CreateField(idFields[-1])
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
for curvec in self.myvectors:
ring = ogr.Geometry(ogr.wkbLinearRing)
for curvert in curvec.myvertices:
ring.AddPoint(curvert.x, curvert.y)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
feature.SetGeometry(poly)
feature.SetField('curvi', float(curvec.myvertices[0].z))
layer.CreateFeature(feature)
feature = None
ds = None
# ----------------------------------------------------------------
# Buffer / parallel / reverse
# ----------------------------------------------------------------
[docs]
def buffer(self, distance: float, resolution: int = 16, inplace: bool = False) -> "zoneModel":
"""Buffer all vectors in the zone.
:param distance: Buffer distance.
:param resolution: Segments per quarter-circle.
:param inplace: Modify in place or return a new zone.
"""
if inplace:
newzone = self
else:
newzone = self._make_zone(name=self.myname)
retmap = list(map(lambda x: x.buffer(distance, resolution, inplace), self.myvectors))
if inplace:
return self
for curvec in retmap:
newzone.add_vector(curvec, forceparent=True)
return newzone
[docs]
def reverse(self, inplace: bool = True) -> "zoneModel":
"""Reverse vertex order in all vectors.
:param inplace: Modify in place or return a copy.
"""
if inplace:
list(map(lambda curvec: curvec.reverse(), self.myvectors))
self.reset_listogl()
return self
else:
newzone = self.deepcopy()
list(map(lambda curvec: curvec.reverse(), newzone.myvectors))
return newzone
[docs]
def add_parallel(self, distance):
"""Add a parallel offset of the active vector to the zone.
:param distance: Offset distance (positive = right, negative = left).
"""
if distance > 0.:
mypl = self.active_vector.parallel_offset(distance, 'right')
elif distance < 0.:
mypl = self.active_vector.parallel_offset(distance, 'left')
else:
mypl = self._make_vector(name=self.active_vector.myname + "_duplicate")
mypl.myvertices = [wolfvertex(cur.x, cur.y, cur.z) for cur in self.active_vector.myvertices]
if mypl is None:
return
mypl.parentzone = self
self.add_vector(mypl)
[docs]
def parallel_active(self, distance):
"""Replace zone contents with the active vector plus left and right offsets.
:param distance: Offset distance.
"""
if self.nbvectors > 1:
self.myvectors = [curv for curv in self.myvectors if curv == self.active_vector]
mypl = self.active_vector.parallel_offset(distance, 'left')
mypr = self.active_vector.parallel_offset(distance, 'right')
if mypl is None or mypr is None:
return
self.add_vector(mypl, 0)
self.add_vector(mypr, 2)
[docs]
def create_cs(self, wa) -> "zoneModel":
"""Create cross-section profiles for each vector by sampling a raster.
:param wa: ``WolfArray`` to sample.
:returns: New zone containing the cross-section vectors.
"""
if wa is None:
logging.warning(_('No array - Cannot create cross-sections !'))
return
new_zone = self._make_zone(name=self.myname + '_cs')
for curvec in self.myvectors:
new_vec = curvec.create_cs(wa)
if new_vec is not None:
new_zone.add_vector(new_vec, forceparent=True)
return new_zone
# ----------------------------------------------------------------
# Deep copy
# ----------------------------------------------------------------
[docs]
def deepcopy_zone(self, name: str = None, parent=None) -> "zoneModel":
"""Return a deep copy of this zone.
:param name: Name for the copy (defaults to ``<name>_copy``).
:param parent: Optional parent for the copy.
"""
if name is None:
name = self.myname + '_copy'
if parent is not None:
copied_zone = self._make_zone(name=name, parent=parent)
else:
copied_zone = self._make_zone(name=name)
copied_zone.myvectors = []
for vec in self.myvectors:
copied_vec = vec.deepcopy_vector(parentzone=copied_zone)
copied_zone.add_vector(copied_vec, forceparent=True)
return copied_zone
[docs]
def deepcopy(self, name: str = None, parent=None) -> "zoneModel":
"""Alias for :pymeth:`deepcopy_zone`."""
return self.deepcopy_zone(name, parent)
# ----------------------------------------------------------------
# Triangulation / multibin (model: requires explicit parameters)
# ----------------------------------------------------------------
[docs]
def create_multibin(self, nb: int = None, nb2: int = 0) -> TriangulationModel:
"""Create a multi-bin triangulation between the zone's vectors.
:param nb: Number of interpolation points per vector.
:param nb2: Number of intermediate vectors between each pair.
:returns: ``TriangulationModel`` or ``None`` on error.
"""
myls = []
for curv in self.myvectors:
myls.append(curv.asshapely_ls())
if nb is None:
logging.warning(_('Parameter nb is required'))
return None
try:
nb = int(nb)
except:
logging.warning(_('Bad parameter nb'))
return None
try:
nb2 = int(nb2)
except:
logging.warning(_('Bad parameter nb2'))
return None
nbvectors = self.nbvectors
s = np.linspace(0., 1., num=nb, endpoint=True)
newls = []
for curls in myls:
newls.append(LineString([curls.interpolate(curs, True) for curs in s]))
if nb2 > 0:
finalls = []
ds = 1. / float(nb2 + 1)
sperp = np.arange(ds, 1., ds)
for j in range(len(newls) - 1):
myls1 = newls[j]
myls2 = newls[j + 1]
xyz1 = np.asarray(myls1.coords[:])
xyz2 = np.asarray(myls2.coords[:])
finalls.append(myls1)
for curds in sperp:
finalls.append(LineString(xyz1 * (1. - curds) + xyz2 * curds))
finalls.append(myls2)
newls = finalls
nbvectors = len(newls)
points = np.zeros((nb * nbvectors, 3), dtype=np.float64)
xyz = []
for curls in newls:
xyz.append(np.asarray(curls.coords[:]))
decal = 0
for i in range(len(xyz[0])):
for k in range(nbvectors):
points[k + decal, :] = xyz[k][i]
decal += nbvectors
triangles = []
nbpts = nbvectors
triangles.append([[i + 0, i + 1, i + nbpts] for i in range(nbpts - 1)])
triangles.append([[i + nbpts, i + 1, i + nbpts + 1] for i in range(nbpts - 1)])
for k in range(1, nb - 1):
decal = k * nbpts
triangles.append([[i + decal, i + decal + 1, i + decal + nbpts] for i in range(nbpts - 1)])
triangles.append([[i + decal + nbpts, i + decal + 1, i + decal + nbpts + 1] for i in range(nbpts - 1)])
triangles = np.asarray(triangles, dtype=np.uint32).reshape([(2 * nbpts - 2) * (nb - 1), 3])
mytri = self._make_triangulation(pts=points, tri=triangles)
mytri.find_minmax(True)
return mytri
[docs]
def create_constrainedDelaunay(self, nb: int = None) -> TriangulationModel:
"""Create a constrained Delaunay triangulation from the zone's vectors.
:param nb: Number of interpolation points per vector (0 = no resampling).
:returns: ``TriangulationModel`` or ``None`` on error.
"""
myls = []
for curv in self.myvectors:
myls.append(curv.linestring)
if nb is None:
logging.warning(_('Parameter nb is required'))
return None
try:
nb = int(nb)
except:
logging.warning(_('Bad parameter nb'))
return None
if nb == 0:
newls = myls
else:
s = np.linspace(0., 1., num=nb, endpoint=True)
newls = [LineString([curls.interpolate(curs, True) for curs in s]) for curls in myls if curls.length > 0.]
xyz = [np.asarray(curls.coords[:]) for curls in newls]
xyz = np.concatenate(xyz)
# Normalize coordinates to improve triangulation robustness (we will restore the original coordinates later)
xmin = xyz[:, 0].min()
ymin = xyz[:, 1].min()
xyz[:, 0] -= xmin
xyz[:, 1] -= ymin
# Remove duplicate vertices and build an index mapping for segments
xyz, indices = np.unique(xyz, axis=0, return_inverse=True)
segments = []
k = 0
for cur in newls:
for i in range(len(cur.coords) - 1):
segments.append([indices[k], indices[k + 1]])
k += 1
# Build the geometry dictionary for triangulation
geom = {'vertices': [[x, y] for x, y in xyz[:, :2]],
'segments': segments}
try:
delaunay = tri_lib.triangulate(geom, 'p')
xyz_dict = {(curxyz[0], curxyz[1]): curxyz[2] for curxyz in xyz}
allvert = []
for curvert in delaunay['vertices']:
x = curvert[0]
y = curvert[1]
z = xyz_dict.get((x, y), 0.)
allvert.append([x + xmin, y + ymin, z])
mytri = self._make_triangulation(
pts=allvert,
tri=[curtri for curtri in delaunay['triangles']]
)
mytri.find_minmax(True)
except Exception as e:
logging.error(_('Error in constrained Delaunay triangulation - {e}'))
return None
return mytri
[docs]
def createmultibin_proj(self, nb: int = None, nb2: int = 0) -> TriangulationModel:
"""Create a projected multi-bin triangulation using a support line.
:param nb: Number of interpolation points per vector.
:param nb2: Number of intermediate vectors between each pair.
:returns: ``TriangulationModel`` or ``None`` on error.
"""
myls = []
nbvectors = self.nbvectors
for curv in self.myvectors:
myls.append(curv.asshapely_ls())
if nb is None:
logging.warning(_('Parameter nb is required'))
return None
try:
nb = int(nb)
except:
logging.warning(_('Bad parameter nb'))
return None
try:
nb2 = int(nb2)
except:
nb2 = 0
s = np.linspace(0., 1., num=nb, endpoint=True)
newls = []
supportls = myls[int(len(myls) / 2)]
supportls = LineString([supportls.interpolate(curs, True) for curs in s])
for curls in myls:
news = [curls.project(Point(curpt[0], curpt[1])) for curpt in supportls.coords]
news.sort()
newls.append(LineString([curls.interpolate(curs) for curs in news]))
if nb2 > 0:
finalls = []
ds = 1. / float(nb2 + 1)
sperp = np.arange(ds, 1., ds)
for j in range(len(newls) - 1):
myls1 = newls[j]
myls2 = newls[j + 1]
xyz1 = np.asarray(myls1.coords[:])
xyz2 = np.asarray(myls2.coords[:])
finalls.append(myls1)
for curds in sperp:
finalls.append(LineString(xyz1 * (1. - curds) + xyz2 * curds))
finalls.append(myls2)
newls = finalls
nbvectors = len(newls)
points = np.zeros((nb * nbvectors, 3), dtype=np.float64)
xyz = []
for curls in newls:
xyz.append(np.asarray(curls.coords[:]))
decal = 0
for i in range(len(xyz[0])):
for k in range(nbvectors):
points[k + decal, :] = xyz[k][i]
decal += nbvectors
triangles = []
nbpts = nbvectors
triangles.append([[i, i + 1, i + nbpts] for i in range(nbpts - 1)])
triangles.append([[i + nbpts, i + 1, i + nbpts + 1] for i in range(nbpts - 1)])
for k in range(1, nb - 1):
decal = k * nbpts
triangles.append([[i + decal, i + decal + 1, i + decal + nbpts] for i in range(nbpts - 1)])
triangles.append([[i + decal + nbpts, i + decal + 1, i + decal + nbpts + 1] for i in range(nbpts - 1)])
triangles = np.asarray(triangles, dtype=np.uint32).reshape([(2 * nbpts - 2) * (nb - 1), 3])
mytri = self._make_triangulation(pts=points, tri=triangles)
mytri.find_minmax(True)
return mytri
# ----------------------------------------------------------------
# Polygon creation from parallels
# ----------------------------------------------------------------
[docs]
def create_polygon_from_parallel(self, ds: float, howmanypoly=1) -> None:
"""Create polygons between left, centre, and right vectors.
The zone must contain exactly 3 vectors.
:param ds: Spacing along the trace.
:param howmanypoly: Number of polygons.
"""
assert self.nbvectors == 3, _('The zone must contain 3 and only 3 vectors')
vecleft = self.myvectors[0]
veccenter = self.myvectors[1]
vecright = self.myvectors[2]
lsl = vecleft.asshapely_ls()
lsr = vecright.asshapely_ls()
lsc = veccenter.asshapely_ls()
nb = int(np.ceil(lsc.length / ds))
sloc = np.linspace(0., 1., nb, endpoint=True)
ptsc = [lsc.interpolate(curs, True) for curs in sloc]
sl = [lsl.project(curs) for curs in ptsc]
sr = [lsr.project(curs) for curs in ptsc]
sc = [lsc.project(curs) for curs in ptsc]
if howmanypoly == 1:
zonepoly = self._make_zone(name='polygons_' + self.myname, parent=self.parent)
self.parent.add_zone(zonepoly)
for i in range(len(sl) - 1):
smean = (sc[i] + sc[i + 1]) / 2.
curvec = self._make_vector(name='poly' + str(i + 1), parentzone=zonepoly)
sublsl = vecleft.substring(sl[i], sl[i + 1], False, False)
sublsr = vecright.substring(sr[i], sr[i + 1], False, False)
if isinstance(sublsl, vectorModel):
vr = sublsr.myvertices.copy()
vr.reverse()
curvec.myvertices = sublsl.myvertices.copy() + vr
for curv in curvec.myvertices:
curv.z = smean
else:
if sublsl.geom_type == 'Point':
curvec.add_vertex(wolfvertex(sublsl.x, sublsl.y, smean))
elif sublsl.geom_type == 'LineString':
xy = np.asarray(sublsl.coords)
for (x, y) in xy:
curvec.add_vertex(wolfvertex(x, y, smean))
if sublsr.geom_type == 'Point':
curvec.add_vertex(wolfvertex(sublsr.x, sublsr.y, smean))
elif sublsr.geom_type == 'LineString':
xy = np.asarray(sublsr.coords)
xy = np.flipud(xy)
for (x, y) in xy:
curvec.add_vertex(wolfvertex(x, y, smean))
curvec.close_force()
zonepoly.add_vector(curvec)
curvec.myprop.legendtext = '{:.2f}'.format(smean)
xy = curvec.asnparray()
curvec.myprop.legendx = np.mean(xy[:, 0])
curvec.myprop.legendy = np.mean(xy[:, 1])
zonepoly.find_minmax(True)
else:
zonepolyleft = self._make_zone(name='polygons_left_' + self.myname, parent=self.parent)
zonepolyright = self._make_zone(name='polygons_right_' + self.myname, parent=self.parent)
self.parent.add_zone(zonepolyleft)
self.parent.add_zone(zonepolyright)
for i in range(len(sl) - 1):
smean = (sc[i] + sc[i + 1]) / 2.
curvecleft = self._make_vector(name='poly' + str(i + 1), parentzone=zonepolyleft)
curvecright = self._make_vector(name='poly' + str(i + 1), parentzone=zonepolyright)
sublsl = vecleft.substring(sl[i], sl[i + 1], False, False)
sublsr = vecright.substring(sr[i], sr[i + 1], False, False)
sublsc = veccenter.substring(sc[i], sc[i + 1], False, False)
if isinstance(sublsl, vectorModel):
vr = sublsr.myvertices.copy()
vr.reverse()
vcr = sublsc.myvertices.copy()
vcr.reverse()
curvecleft.myvertices = sublsl.myvertices.copy() + vcr
curvecright.myvertices = sublsc.myvertices.copy() + vr
for curv in curvecleft.myvertices:
curv.z = smean
for curv in curvecright.myvertices:
curv.z = smean
else:
if sublsl.geom_type == 'Point':
curvecleft.add_vertex(wolfvertex(sublsl.x, sublsl.y, smean))
elif sublsl.geom_type == 'LineString':
xy = np.asarray(sublsl.coords)
for (x, y) in xy:
curvecleft.add_vertex(wolfvertex(x, y, smean))
if sublsc.geom_type == 'Point':
curvecleft.add_vertex(wolfvertex(sublsc.x, sublsc.y, smean))
elif sublsc.geom_type == 'LineString':
xy = np.asarray(sublsc.coords)
xy = np.flipud(xy)
for (x, y) in xy:
curvecleft.add_vertex(wolfvertex(x, y, smean))
if sublsc.geom_type == 'Point':
curvecright.add_vertex(wolfvertex(sublsc.x, sublsc.y, smean))
elif sublsc.geom_type == 'LineString':
xy = np.asarray(sublsc.coords)
for (x, y) in xy:
curvecright.add_vertex(wolfvertex(x, y, smean))
if sublsr.geom_type == 'Point':
curvecright.add_vertex(wolfvertex(sublsr.x, sublsr.y, smean))
elif sublsr.geom_type == 'LineString':
xy = np.asarray(sublsr.coords)
xy = np.flipud(xy)
for (x, y) in xy:
curvecright.add_vertex(wolfvertex(x, y, smean))
curvecleft.close_force()
curvecright.close_force()
curvecleft.myprop.legendtext = '{:.2f}'.format(smean)
curvecright.myprop.legendtext = '{:.2f}'.format(smean)
xy = curvecleft.asnparray()
curvecleft.myprop.legendx = np.mean(xy[:, 0])
curvecleft.myprop.legendy = np.mean(xy[:, 1])
xy = curvecright.asnparray()
curvecright.myprop.legendx = np.mean(xy[:, 0])
curvecright.myprop.legendy = np.mean(xy[:, 1])
zonepolyleft.add_vector(curvecleft)
zonepolyright.add_vector(curvecright)
zonepolyleft.find_minmax(True)
zonepolyright.find_minmax(True)
self._fill_structure()
[docs]
def create_sliding_polygon_from_parallel(self,
poly_length: float,
ds_sliding: float,
farthest_parallel: float,
interval_parallel: float = None,
intersect=None,
howmanypoly=1,
eps_offset: float = 0.25):
"""Create sliding polygons along a trace with parallel offsets.
:param poly_length: Length of each polygon along the trace.
:param ds_sliding: Sliding step between successive polygons.
:param farthest_parallel: Maximum offset distance.
:param interval_parallel: Interval between parallel offsets.
:param intersect: Optional intersecting vector to limit the trace.
:param howmanypoly: Number of polygons per position.
:param eps_offset: Small offset applied at the starting position.
"""
assert self.nbvectors == 1, _('The zone must contain 1 and only 1 vector')
vecleft: dict[str, vectorModel] = {}
vecright: dict[str, vectorModel] = {}
veccenter = self.myvectors[0]
veccenter.update_lengths()
logging.info(_('Length of the center vector: {}').format(veccenter.length2D))
myparallels = self._make_zone()
if interval_parallel is None:
logging.warning(_('Interval between parallels is not defined --> set to farthest_parallel'))
interval_parallel = farthest_parallel
if interval_parallel > farthest_parallel:
logging.warning(_('Interval between parallels is greater than farthest_parallel --> set to farthest_parallel'))
interval_parallel = farthest_parallel
all_par = np.arange(0, farthest_parallel, interval_parallel)[1:]
all_par = np.concatenate((all_par, [farthest_parallel]))
logging.info(_('All parallel distances: {}').format(all_par))
for curpar in tqdm(all_par):
vecleft[curpar] = veccenter.parallel_offset(curpar, 'left')
vecright[curpar] = veccenter.parallel_offset(curpar, 'right')
myparallels.add_vector(vecleft[curpar], forceparent=True)
myparallels.add_vector(vecright[curpar], forceparent=True)
if isinstance(intersect, zoneModel):
for curint in intersect.myvectors:
if not curint.used:
continue
curint1 = curint.parallel_offset(-eps_offset / 2., side='left')
curint2 = curint.parallel_offset(eps_offset / 2., side='right')
pt, dist = vecleft[curpar].intersection(curint1, eval_dist=True, force_single=True)
if pt is not None:
logging.debug(_('Intersection found on left parallel at distance {}').format(dist))
dist2 = curint1.linestring.project(pt)
subs = curint1.substring(0., dist2, is3D=False, adim=False)
subs.reverse()
subs2 = curint2.substring(0., dist2, is3D=False, adim=False)
vec1 = vecleft[curpar].substring(0., dist, is3D=False, adim=False)
vec2 = vecleft[curpar].substring(dist, vecleft[curpar].length2D, is3D=False, adim=False)
vecleft[curpar].myvertices = vec1.myvertices.copy() + subs.myvertices.copy() + subs2.myvertices.copy() + vec2.myvertices.copy()
vecleft[curpar].find_minmax()
vecleft[curpar].update_lengths()
vecleft[curpar].reset_linestring()
curint1.reset_linestring()
curint2.reset_linestring()
pt, dist = vecright[curpar].intersection(curint1, eval_dist=True, force_single=True)
if pt is not None:
logging.debug(_('Intersection found on right parallel at distance {}').format(dist))
dist2 = curint1.linestring.project(pt)
subs = curint1.substring(0., dist2, is3D=False, adim=False)
subs2 = curint2.substring(0., dist2, is3D=False, adim=False)
subs2.reverse()
vec1 = vecright[curpar].substring(0., dist, is3D=False, adim=False)
vec2 = vecright[curpar].substring(dist, vecright[curpar].length2D, is3D=False, adim=False)
vecright[curpar].myvertices = vec1.myvertices.copy() + subs2.myvertices.copy() + subs.myvertices.copy() + vec2.myvertices.copy()
vecright[curpar].find_minmax()
vecright[curpar].update_lengths()
vecright[curpar].reset_linestring()
curint1.reset_linestring()
curint2.reset_linestring()
lsl = {key: vec.asshapely_ls() for key, vec in vecleft.items()}
lsr = {key: vec.asshapely_ls() for key, vec in vecright.items()}
lsc = veccenter.asshapely_ls()
nb = int(np.ceil(lsc.length / float(ds_sliding)))
sloc = np.asarray([float(ds_sliding) * cur for cur in range(nb)])
sloc2 = sloc + float(poly_length)
sloc2[sloc2 > veccenter.length2D] = veccenter.length2D
ptsc = [veccenter.interpolate(curs, is3D=False, adim=False) for curs in sloc]
ptsc2 = [veccenter.interpolate(curs, is3D=False, adim=False) for curs in sloc2]
sc = [lsc.project(Point(curs.x, curs.y)) for curs in ptsc]
sc2 = [lsc.project(Point(curs.x, curs.y)) for curs in ptsc2]
sl = {}; sr = {}; sl2 = {}; sr2 = {}
ptl = {}; ptl2 = {}; ptr = {}; ptr2 = {}
curpts = ptsc
for key, ls in lsl.items():
sl[key] = [ls.project(Point(curs.x, curs.y)) for curs in curpts]
ptl[key] = [ls.interpolate(curs) for curs in sl[key]]
curpts = ptl[key]
curpts = ptsc2
for key, ls in lsl.items():
sl2[key] = [ls.project(Point(curs.x, curs.y)) for curs in curpts]
ptl2[key] = [ls.interpolate(curs) for curs in sl2[key]]
curpts = ptl2[key]
curpts = ptsc
for key, ls in lsr.items():
sr[key] = [ls.project(Point(curs.x, curs.y)) for curs in curpts]
ptr[key] = [ls.interpolate(curs) for curs in sr[key]]
curpts = ptr[key]
curpts = ptsc2
for key, ls in lsr.items():
sr2[key] = [ls.project(Point(curs.x, curs.y)) for curs in curpts]
ptr2[key] = [ls.interpolate(curs) for curs in sr2[key]]
curpts = ptr2[key]
if howmanypoly == 1:
zonepoly = self._make_zone(name='polygons_' + self.myname, parent=self.parent)
self.parent.add_zone(zonepoly, forceparent=True)
for i in range(nb):
ptc1 = sc[i]; ptc2 = sc2[i]
pt1 = [cur[i] for cur in sl.values()]
pt2 = [cur[i] for cur in sl2.values()]
pt3 = [cur[i] for cur in sr.values()]
pt4 = [cur[i] for cur in sr2.values()]
smean = (ptc1 + ptc2) / 2.
curvec = self._make_vector(name='poly' + str(i), parentzone=zonepoly)
sublsl = vecleft[farthest_parallel].substring(pt1[-1], pt2[-1], is3D=False, adim=False)
sublsr = vecright[farthest_parallel].substring(pt3[-1], pt4[-1], is3D=False, adim=False)
sublsr.reverse()
sublsc = veccenter.substring(ptc1, ptc2, is3D=False, adim=False)
upl = [wolfvertex(pt[i].x, pt[i].y) for pt in ptl.values()]
upr = [wolfvertex(pt[i].x, pt[i].y) for pt in ptr.values()]
upr.reverse()
downl = [wolfvertex(pt[i].x, pt[i].y) for pt in ptl2.values()]
downl.reverse()
downr = [wolfvertex(pt[i].x, pt[i].y) for pt in ptr2.values()]
curvec.myvertices = sublsl.myvertices.copy() + downl[1:].copy() + [sublsc.myvertices[-1].copy()] + downr[:-1].copy() + sublsr.myvertices.copy() + upr[1:].copy() + [sublsc.myvertices[0].copy()] + upl[:-1].copy()
for curvert in curvec.myvertices:
curvert.z = smean
curvec.close_force()
zonepoly.add_vector(curvec, forceparent=True)
zonepoly.find_minmax(True)
else:
zonepolyleft = self._make_zone(name='polygons_left_' + self.myname, parent=self.parent)
zonepolyright = self._make_zone(name='polygons_right_' + self.myname, parent=self.parent)
self.parent.add_zone(zonepolyleft, forceparent=True)
self.parent.add_zone(zonepolyright, forceparent=True)
for i in range(nb):
ptc1 = sc[i]; ptc2 = sc2[i]
pt1 = [cur[i] for cur in sl.values()]
pt2 = [cur[i] for cur in sl2.values()]
pt3 = [cur[i] for cur in sr.values()]
pt4 = [cur[i] for cur in sr2.values()]
smean = (ptc1 + ptc2) / 2.
curvecleft = self._make_vector(name='poly' + str(i + 1), parentzone=zonepolyleft)
curvecright = self._make_vector(name='poly' + str(i + 1), parentzone=zonepolyright)
sublsl = vecleft[farthest_parallel].substring(pt1[-1], pt2[-1], is3D=False, adim=False)
sublsr = vecright[farthest_parallel].substring(pt3[-1], pt4[-1], is3D=False, adim=False)
sublsr.reverse()
sublsc = veccenter.substring(ptc1, ptc2, is3D=False, adim=False)
sublscr = sublsc.deepcopy()
sublscr.reverse()
upl = [wolfvertex(pt[i].x, pt[i].y) for pt in ptl.values()]
upr = [wolfvertex(pt[i].x, pt[i].y) for pt in ptr.values()]
upr.reverse()
downl = [wolfvertex(pt[i].x, pt[i].y) for pt in ptl2.values()]
downl.reverse()
downr = [wolfvertex(pt[i].x, pt[i].y) for pt in ptr2.values()]
curvecleft.myvertices = sublsl.myvertices.copy() + downl[1:-1].copy() + sublscr.myvertices.copy() + upl[1:-1].copy()
curvecright.myvertices = sublsc.myvertices.copy() + downr[1:-1].copy() + sublsr.myvertices.copy() + upr[1:-1].copy()
for curvert in curvecleft.myvertices:
curvert.z = smean
for curvert in curvecright.myvertices:
curvert.z = smean
curvecleft.close_force()
curvecright.close_force()
zonepolyleft.add_vector(curvecleft)
zonepolyright.add_vector(curvecright)
zonepolyleft.find_minmax(True)
zonepolyright.find_minmax(True)
self._fill_structure()
return myparallels
# ----------------------------------------------------------------
# Linked-array extraction / stats
# ----------------------------------------------------------------
[docs]
def get_values_linked_polygons(self, linked_arrays, stats=True) -> dict:
"""Extract raster values inside each vector's polygon.
:param linked_arrays: Arrays to query.
:param stats: Compute statistics on the results.
"""
exit_flag = True
for curarray in linked_arrays:
if curarray.plotted:
exit_flag = False
if exit_flag:
return None
vals = {idx: {'values': curpol.get_values_linked_polygon(linked_arrays)} for idx, curpol in enumerate(self.myvectors)}
if stats:
self._stats_values(vals)
return vals
[docs]
def get_all_values_linked_polygon(self, linked_arrays, stats=True,
key_idx_names: Literal['idx', 'name'] = 'idx', getxy=False) -> dict:
"""Extract all raster values inside each vector's polygon.
:param linked_arrays: Arrays to query.
:param stats: Compute statistics.
:param key_idx_names: Key type for result dict (*'idx'* or *'name'*).
:param getxy: Also return XY coordinates.
"""
exit_flag = True
for curarray in linked_arrays:
if curarray.plotted:
exit_flag = False
if exit_flag:
return None
if key_idx_names == 'idx':
vals = {idx: {'values': curpol.get_all_values_linked_polygon(linked_arrays, getxy=getxy)} for idx, curpol in enumerate(self.myvectors)}
else:
vals = {curpol.myname: {'values': curpol.get_all_values_linked_polygon(linked_arrays, getxy=getxy)} for idx, curpol in enumerate(self.myvectors)}
return vals
[docs]
def _stats_values(self, vals: dict):
"""Compute median, mean, min, max, p5, p95 statistics on extracted values.
:param vals: Dict of ``{polygon_id: {values: ...}}``.
"""
for curpol in vals.values():
medianlist = curpol['median'] = []
meanlist = curpol['mean'] = []
minlist = curpol['min'] = []
maxlist = curpol['max'] = []
p95 = curpol['p95'] = []
p5 = curpol['p5'] = []
for curval in curpol['values']:
if curval[1] is not None:
if curval[0] is not None and len(curval[0]) > 0:
medianlist.append((np.median(curval[0]), np.median(curval[1])))
meanlist.append((np.mean(curval[0]), np.mean(curval[1])))
minlist.append((np.min(curval[0]), np.min(curval[1])))
maxlist.append((np.max(curval[0]), np.max(curval[1])))
p95.append((np.percentile(curval[0], 95), np.percentile(curval[1], 95)))
p5.append((np.percentile(curval[0], 5), np.percentile(curval[1], 5)))
else:
medianlist.append((None, None))
meanlist.append((None, None))
minlist.append((None, None))
maxlist.append((None, None))
p95.append((None, None))
p5.append((None, None))
else:
if curval[0] is not None and len(curval[0]) > 0:
medianlist.append(np.median(curval[0]))
meanlist.append(np.mean(curval[0]))
minlist.append(np.min(curval[0]))
maxlist.append(np.max(curval[0]))
p95.append(np.percentile(curval[0], 95))
p5.append(np.percentile(curval[0], 5))
else:
medianlist.append(None)
meanlist.append(None)
minlist.append(None)
maxlist.append(None)
p95.append(None)
p5.append(None)
# ====================================================================
# ZonesModel
# ====================================================================
[docs]
class ZonesModel:
"""Pure-data Zones: a list of zoneModel instances plus I/O and geometry helpers.
GUI subclass ``Zones`` in :mod:`wolfhece.PyVertexvectors` adds
wx.Frame, OpenGL display-list support, tree-list UI and wx integration.
"""
[docs]
myzones: list[zoneModel]
def __init__(self,
filename: str | Path = '',
ox: float = 0.,
oy: float = 0.,
tx: float = 0.,
ty: float = 0.,
is2D: bool = True,
idx: str = '',
colname: str = None,
bbox: Polygon = None,
find_minmax_init: bool = True,
colors: dict = None) -> None:
"""
Data-only initialiser for vectorial information.
GUI subclass ``Zones`` in :mod:`wolfhece.PyVertexvectors` adds
wx.Frame, OpenGL display-list support and wx integration.
:param filename: file path to read (.vec, .vecz, .dxf, .shp, .gpkg, .gdb)
:param ox: X origin
:param oy: Y origin
:param tx: X translation
:param ty: Y translation
:param is2D: if True, vectors are 2D
:param idx: identifier
:param colname: column name for shapefile import
:param bbox: bounding box for filtering
:param find_minmax_init: if True, find min/max after loading
:param colors: dict of colors for zone colorization
"""
[docs]
self.active_vector = None
[docs]
self.active_zone = None
[docs]
self.last_active = None
[docs]
self.filename = str(filename)
[docs]
self.init_struct = True
[docs]
self._first_find_minmax: bool = True
self.tx = tx
self.ty = ty
self.myzones = []
[docs]
self._move_start = None
[docs]
self._rotation_center = None
[docs]
self._rotation_step = None
only_firstlast = False
if self.filename != '':
_filename = str(self.filename).strip()
if _filename.startswith('http:') or _filename.startswith('https:'):
try:
from ..pydownloader import download_file
self.filename = str(download_file(_filename))
except Exception as e:
logging.error(_('Error while downloading file: %s') % e)
return
if self.filename.endswith('.dxf'):
self.is2D = False
self.import_dxf(self.filename)
only_firstlast = True
elif self.filename.endswith('.shp'):
self.is2D = False
self.import_shapefile(self.filename, bbox=bbox, colname=colname)
only_firstlast = True
elif self.filename.endswith('.gpkg'):
self.is2D = False
self.import_gpkg(self.filename, bbox=bbox)
only_firstlast = True
elif Path(filename).is_dir() and self.filename.endswith('.gdb'):
self.is2D = False
self.import_gdb(self.filename, bbox=bbox)
only_firstlast = True
elif self.filename.endswith('.vec') or self.filename.endswith('.vecz'):
if self.filename.endswith('.vecz'):
self.is2D = False
f = open(self.filename, 'r')
lines = f.read().splitlines()
f.close()
try:
tx, ty = lines[0].split()
except:
tx, ty = lines[0].split(',')
self.tx = float(tx)
self.ty = float(ty)
tmp_nbzones = int(lines[1])
curstart = 2
for i in range(tmp_nbzones):
curzone = self._make_zone(lines=lines[curstart:], parent=self, is2D=self.is2D)
self.myzones.append(curzone)
curstart += curzone._nblines()
if Path(self.filename + '.extra').exists():
with open(self.filename + '.extra', 'r') as f:
lines = f.read().splitlines()
try:
nblines = len(lines)
i = 0
idx_zone = 0
while i < nblines:
curzone = self.myzones[idx_zone]
assert curzone.myname == lines[i], _('Error while reading extra properties of {}'.format(self.filename))
i += 1
ret = curzone.load_extra(lines[i:])
i += ret
idx_zone += 1
except:
logging.warning(_('Error while reading extra properties of {}'.format(self.filename)))
if find_minmax_init:
logging.info(_('Finding min and max values'))
self.find_minmax(True, only_firstlast)
if colors is not None:
self.colorize_data(colors, filled=True)
[docs]
def prep_listogl(self):
"""No-op in the model; GUI override prepares OpenGL display lists."""
pass
[docs]
def reset_listogl(self):
"""No-op in the model; GUI override resets OpenGL display lists."""
pass
[docs]
def _make_zone(self, **kwargs):
"""Factory: create a zoneModel. GUI override returns ``zone``."""
return zoneModel(**kwargs)
[docs]
def _make_zones(self, **kwargs):
"""Factory: create a ZonesModel. GUI override returns ``Zones``."""
return ZonesModel(**kwargs)
[docs]
def _make_vector(self, **kwargs):
"""Factory: create a vectorModel. GUI override returns ``vector``."""
return vectorModel(**kwargs)
def __getstate__(self):
""" Get the state of the object for pickling """
state = self.__dict__.copy()
# Remove unpicklable entries
if 'mapviewer' in state:
del state['mapviewer']
return state
def __setstate__(self, state):
""" Set the state of the object for unpickling """
self.__dict__.update(state)
@property
[docs]
def mynames(self) -> list[str]:
""" Return the names of all zones """
return [curzone.myname for curzone in self.myzones]
[docs]
def check_if_interior_exists(self):
""" Check if the zone has at least one vector with interior points """
for curzone in self.myzones:
curzone.check_if_interior_exists()
[docs]
def add_values(self, key:str, values:np.ndarray | dict):
"""Add values to the zones.
:param key: Value identifier.
:param values: Either a dict ``{zone_name: ndarray}`` or an ndarray
whose first axis length must match :pyattr:`nbzones`.
"""
if isinstance(values, dict):
for k, val in values.items():
if not isinstance(val, np.ndarray):
val = np.asarray(val)
if k in self.mynames:
self[k].add_values(key, val)
elif isinstance(values, np.ndarray):
if values.shape[0] != self.nbzones:
logging.warning(_('Number of values does not match the number of zones'))
return
for idx, curzone in enumerate(self.myzones):
curzone.add_values(key, values[idx])
[docs]
def get_values(self, key:str) -> np.ndarray:
"""Get values from the zones.
:param key: Value identifier.
:return: Array of values, one element per zone.
"""
return np.asarray([curzone.get_values(key) for curzone in self.myzones])
[docs]
def set_colors_from_value(self, key:str, cmap:wolfpalette | str | Colormap | cm.ScalarMappable, vmin:float = 0., vmax:float = 1.):
"""Set colours for all zones based on stored values.
:param key: Value identifier used to look up per-zone data.
:param cmap: Colour map (palette, name, Colormap or ScalarMappable).
:param vmin: Lower bound of the value range.
:param vmax: Upper bound of the value range.
"""
if isinstance(cmap, str):
cmap = cm.get_cmap(cmap)
for curzone in self.myzones:
curzone.set_colors_from_value(key, cmap, vmin, vmax)
self.prep_listogl()
[docs]
def set_color(self, color: int | RGB | str | list | tuple):
"""Set a uniform colour for all zones.
:param color: Colour specification (integer, RGB, name, or list/tuple).
"""
for curzone in self.myzones:
curzone.set_color(color)
self.prep_listogl()
[docs]
def set_linewidth(self, linewidth:int | float):
"""Set line width for all zones.
:param linewidth: Width in pixels.
"""
for curzone in self.myzones:
curzone.set_linewidth(linewidth)
self.prep_listogl()
[docs]
def set_alpha(self, alpha:int):
"""Set transparency for all zones.
:param alpha: Alpha value (0–255).
"""
for curzone in self.myzones:
curzone.set_alpha(alpha)
self.prep_listogl()
[docs]
def set_filled(self, filled:bool):
"""Set filled rendering for all zones.
:param filled: Whether polygons should be drawn filled.
"""
for curzone in self.myzones:
curzone.set_filled(filled)
self.prep_listogl()
[docs]
def check_if_open(self):
""" Check if the vectors in the zone are open """
for curzone in self.myzones:
curzone.check_if_open()
[docs]
def concatenate_all_vectors(self) -> list["vectorModel"]:
""" Concatenate all vectors in the zones """
ret = []
for curzone in self.myzones:
ret.extend(curzone.myvectors)
return ret
[docs]
def prepare_shapely(self):
""" Prepare shapely objects for all vectors in zones """
allvec = self.concatenate_all_vectors()
list(map(lambda x: x.prepare_shapely(True), allvec))
[docs]
def filter_contains(self, others:list["vectorModel"]) -> "ZonesModel":
"""Create a new ZonesModel with vectors from *others* contained in the zones.
:param others: List of vectors or a ZonesModel to filter.
:return: A new ``ZonesModel`` with the contained vectors.
"""
if isinstance(others, ZonesModel):
allvec = others.concatenate_all_vectors()
elif isinstance(others, list):
allvec = others
else:
logging.warning(_('Unknown type for others'))
return None
centroids = [curvec.centroid for curvec in allvec]
newzones = self._make_zones()
for curzone in self.myzones:
if curzone.nbvectors != 1:
logging.warning(_('Zone {} has more than one vector'.format(curzone.myname)))
continue
poly = curzone[0].polygon
# element-wise comparison
contains = list(map(lambda x: poly.contains(x), centroids))
newzone = self._make_zone(name=curzone.myname, parent=newzones)
newzone.myvectors = list(map(lambda x: allvec[x], np.where(contains)[0]))
newzones.add_zone(newzone, forceparent=True)
return newzones
@property
[docs]
def areas(self):
""" List of areas of all zones """
return [curzone.areas for curzone in self.myzones]
[docs]
def buffer(self, distance:float, resolution:int = 16, inplace:bool = True):
"""Buffer all zones.
:param distance: Buffer distance in world units.
:param resolution: Number of segments to approximate a quarter circle.
:param inplace: If True, modify in place; otherwise return a new object.
"""
if inplace:
newzones = self
else:
newzones = self._make_zones()
retmap = list(map(lambda x: x.buffer(distance, resolution, inplace=inplace), self.myzones))
for curzone in retmap:
newzones.add_zone(curzone, forceparent=True)
newzones.find_minmax(True)
if inplace:
return self
return newzones
[docs]
def set_cache(self):
""" Set cache for all zones """
for curzone in self.myzones:
curzone.set_cache()
[docs]
def clear_cache(self):
""" Clear cache for all zones """
for curzone in self.myzones:
curzone.clear_cache()
[docs]
def move(self, dx:float, dy:float, use_cache:bool = True, inplace:bool = True):
"""Move all zones.
:param dx: Translation along X.
:param dy: Translation along Y.
:param use_cache: If True, move relative to cached positions.
:param inplace: If True, modify in place; otherwise return a deep copy.
"""
if self._move_step is not None:
dx = np.round(dx/self._move_step)*self._move_step
dy = np.round(dy/self._move_step)*self._move_step
if inplace:
for curzone in self.myzones:
curzone.move(dx, dy, use_cache=use_cache, inplace=inplace)
return self
else:
newzones = self.deepcopy_zones()
newzones.move(dx, dy, use_cache=False, inplace=True)
newzones.find_minmax(True)
return newzones
[docs]
def rotate(self, angle:float, center:Point = None, use_cache:bool = True, inplace:bool = True):
"""Rotate all zones.
:param angle: Rotation angle in radians.
:param center: Centre of rotation (default: origin).
:param use_cache: If True, rotate from cached positions.
:param inplace: If True, modify in place; otherwise return a deep copy.
"""
if self._rotation_step is not None:
angle = np.round(angle/self._rotation_step)*self._rotation_step
if inplace:
for curzone in self.myzones:
curzone.rotate(angle, center, use_cache=use_cache, inplace=inplace)
return self
else:
newzones = self.deepcopy_zones()
newzones.rotate(angle, center, use_cache=False, inplace=True)
newzones.find_minmax(True)
return newzones
[docs]
def rotate_xy(self, angle:float, center:Point = None, use_cache:bool = True, inplace:bool = True):
"""Rotate all zones (XY plane only).
:param angle: Rotation angle in radians.
:param center: Centre of rotation (default: origin).
:param use_cache: If True, rotate from cached positions.
:param inplace: If True, modify in place; otherwise return a deep copy.
"""
if inplace:
for curzone in self.myzones:
curzone.rotate_xy(angle, center, use_cache=use_cache, inplace=inplace)
return self
else:
newzones = self.deepcopy_zones()
newzones.rotate_xy(angle, center, use_cache=False, inplace=True)
newzones.find_minmax(True)
return newzones
[docs]
def force_unique_zone_name(self):
"""
Check if all zones have a unique id
If not, the id will be set to the index of the zone in the list
"""
names = [curzone.myname for curzone in self.myzones]
unique_names = set(names)
if len(unique_names) != len(names):
for idx, curzone in enumerate(self.myzones):
if names.count(curzone.myname)>1:
curzone.myname += '_'+str(idx)
[docs]
def set_legend_visible(self, visible:bool=True):
"""Set the visibility of the legend for all zones.
:param visible: Whether the legend should be visible.
"""
for curzone in self.myzones:
curzone.set_legend_visible(visible)
[docs]
def set_legend_text(self, text:str):
"""Set the legend text for all zones.
:param text: Legend text string.
"""
for curzone in self.myzones:
curzone.set_legend_text(text)
[docs]
def set_legend_color(self, color:str | int | tuple | list | RGB):
"""Set the legend colour for all zones.
:param color: Colour specification.
"""
for curzone in self.myzones:
curzone.set_legend_color(color)
[docs]
def set_legend(self, x:float | str, y:float | str, text:str, visible:bool = True):
"""Set legend position, text and visibility for all zones.
:param x: X-coordinate of the legend.
:param y: Y-coordinate of the legend.
:param text: Legend text string.
:param visible: Whether the legend should be visible.
"""
self.set_legend_position(x, y)
self.set_legend_text(text)
self.set_legend_visible(visible)
[docs]
def set_legend_text_from_values(self, key:str):
"""Set legend text for all zones from stored values.
:param key: Value identifier.
"""
for curzone in self.myzones:
curzone.set_legend_text_from_values(key)
[docs]
def set_legend_to_centroid(self):
"""
Set the legend to the centroid of the zones
"""
for curzone in self.myzones:
curzone.set_legend_to_centroid()
[docs]
def set_legend_position(self, x, y):
"""Set legend position for all zones.
:param x: X-coordinate.
:param y: Y-coordinate.
"""
for curzone in self.myzones:
curzone.set_legend_position(x, y)
# ----------------------------------------------------------------
# Shader-property helpers (propagate to all zones)
# ----------------------------------------------------------------
[docs]
def set_glow(self, enabled: bool = True, width: float = 0.4,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.4):
"""Configure line glow for all zones.
:param enabled: Enable or disable glow.
:param width: Glow width as fraction of half-width (0–1).
:param color: RGB colour tuple [0, 1].
:param alpha: Glow opacity [0, 1].
"""
for curzone in self.myzones:
curzone.set_glow(enabled, width, color, alpha)
self.prep_listogl()
[docs]
def set_dash(self, enabled: bool = True, dash_length: float = 10.0,
gap_length: float = 5.0):
"""Configure dashed lines for all zones.
:param enabled: Enable or disable dashes.
:param dash_length: Dash length in world units.
:param gap_length: Gap length in world units.
"""
for curzone in self.myzones:
curzone.set_dash(enabled, dash_length, gap_length)
self.prep_listogl()
[docs]
def set_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure line animation for all zones.
:param mode: 0=none, 1=pulse, 2=marching ants.
:param speed: Speed multiplier.
"""
for curzone in self.myzones:
curzone.set_animation(mode, speed)
self.prep_listogl()
[docs]
def set_fill_animation(self, mode: int = 1, speed: float = 1.0,
center_index: int = 0,
start_angle: float = 90.0):
"""Configure filled-polygon animation for all zones.
:param mode: 0=none, 1=pulse, 2=sweep, 3=ripple,
4=radial progressive, 5=clockwise clock, 6=counter-clockwise clock.
:param speed: Speed multiplier.
:param center_index: Vertex index used as the radial/clock centre.
:param start_angle: Clock fill start angle in degrees.
"""
for curzone in self.myzones:
curzone.set_fill_animation(mode, speed, center_index, start_angle)
self.prep_listogl()
[docs]
def set_fill_color(self, color=None):
"""Set independent fill colour for all zones.
:param color: RGB tuple, packed int, or None to inherit main colour.
"""
for curzone in self.myzones:
curzone.set_fill_color(color)
self.prep_listogl()
[docs]
def set_contour_color(self, color=None, width: float = 1.0,
enabled: bool = True):
"""Set contour colour/width for all zones.
:param color: RGB tuple, packed int, or None to inherit main colour.
:param width: Contour line width.
:param enabled: Whether to draw the contour pass.
"""
for curzone in self.myzones:
curzone.set_contour_color(color, width, enabled)
self.prep_listogl()
[docs]
def set_join_style(self, style: int = 0, size: float = 0.5,
size_mode: int = 0):
"""Configure join rendering for all zones.
:param style: 0=miter, 1=bevel, 2=round, 3=fillet.
:param size: Join size.
:param size_mode: 0=fraction, 1=world distance.
"""
for curzone in self.myzones:
curzone.set_join_style(style, size, size_mode)
self.prep_listogl()
[docs]
def set_width_in_pixels(self, pixels: bool = True):
"""Toggle pixel vs world-unit line width for all zones.
:param pixels: ``True`` for pixel widths.
"""
for curzone in self.myzones:
curzone.set_width_in_pixels(pixels)
self.prep_listogl()
[docs]
def set_text_along(self, text: str, enabled: bool = True,
offset: float = 0.0, perp: float = 0.0,
alignment: str = 'center'):
"""Configure text along polyline for all zones.
:param text: Text to display.
:param enabled: Enable or disable.
:param offset: Shift start along polyline (world units).
:param perp: Perpendicular offset (world units).
:param alignment: 'left', 'center', 'right'.
"""
for curzone in self.myzones:
curzone.set_text_along(text, enabled, offset, perp, alignment)
[docs]
def set_legend_glow(self, enabled: bool = True, width: float = 0.15,
color: tuple = (1.0, 1.0, 1.0), alpha: float = 0.5):
"""Configure legend glow for all zones.
:param enabled: Enable or disable legend glow.
:param width: SDF threshold offset.
:param color: RGB colour tuple [0, 1].
:param alpha: Glow opacity [0, 1].
"""
for curzone in self.myzones:
curzone.set_legend_glow(enabled, width, color, alpha)
[docs]
def set_legend_animation(self, mode: int = 1, speed: float = 1.0):
"""Configure legend animation for all zones.
:param mode: 0=none, 1=pulse, 2=wave, 3=typewriter.
:param speed: Speed multiplier.
"""
for curzone in self.myzones:
curzone.set_legend_animation(mode, speed)
[docs]
def set_legend_style(self, smoothing: float = 1.0,
alignment: str = 'left',
line_spacing: float = 1.2):
"""Configure legend SDF text style for all zones.
:param smoothing: SDF AA multiplier.
:param alignment: 'left', 'center', 'right'.
:param line_spacing: Multiline spacing multiplier.
"""
for curzone in self.myzones:
curzone.set_legend_style(smoothing, alignment, line_spacing)
@property
[docs]
def nbzones(self):
"""Number of zones in the collection."""
return len(self.myzones)
[docs]
def import_shapefile(self, fn:str,
bbox:Polygon = None,
colname:str = None,
value_columns: str | list[str] | tuple[str, ...] | None = None,
share_multipart_values: bool = False):
"""Import a shapefile using geopandas.
The entire shapefile is loaded as a single zone.
:param fn: Path to the shapefile.
:param bbox: Optional bounding box to filter features.
:param colname: Column name used to build zone names.
:param value_columns: Columns to store into vector properties.
- ``None``: do not import attribute values (default).
- ``'all'``: import all non-geometry columns.
- list/tuple: import selected columns.
:param share_multipart_values: If True, multipart geometries share
the same underlying ``_values`` dictionary across split vectors.
"""
logging.info(_('Importing shapefile {}'.format(fn)))
content = gpd.read_file(fn, bbox=bbox)
self.import_GeoDataFrame(content=content,
colname=colname,
value_columns=value_columns,
share_multipart_values=share_multipart_values)
[docs]
def import_GeoDataFrame(self, content:gpd.GeoDataFrame,
bbox:Polygon = None,
colname:str = None,
value_columns: str | list[str] | tuple[str, ...] | None = None,
share_multipart_values: bool = False):
"""Import a GeoDataFrame.
:param content: GeoDataFrame to import.
:param bbox: Optional bounding box to filter features.
:param colname: Column name used to build zone names.
:param value_columns: Columns to store into vector properties.
- ``None``: do not import attribute values (default).
- ``'all'``: import all non-geometry columns.
- list/tuple: import selected columns.
:param share_multipart_values: If True, multipart geometries share
the same underlying ``_values`` dictionary across split vectors.
"""
logging.info(_('Importing GeoDataFrame'))
if bbox is not None:
# filter content
content = content.cx[bbox.bounds[0]:bbox.bounds[2], bbox.bounds[1]:bbox.bounds[3]]
logging.info(_('Converting DataFrame into zones'))
if colname is not None and colname not in content.columns:
logging.warning(_('Column {} not found in the DataFrame'.format(colname)))
logging.info(_('We are using the available known columns'))
colname = ''
columns = list(content.columns)
col_to_pos = {curcol: idx + 1 for idx, curcol in enumerate(columns)}
geometry_pos = col_to_pos.get('geometry')
if geometry_pos is None:
logging.warning(_('No geometry column found in the GeoDataFrame'))
self.myzones = []
return
# Resolve columns to import as vector-property values.
if value_columns is None:
selected_value_columns = []
elif isinstance(value_columns, str):
if value_columns.lower() == 'all':
selected_value_columns = [cur for cur in content.columns if cur != 'geometry']
elif value_columns.strip() == '':
selected_value_columns = []
else:
selected_value_columns = [value_columns]
else:
selected_value_columns = list(value_columns)
missing_columns = [cur for cur in selected_value_columns if cur not in content.columns]
if len(missing_columns) > 0:
logging.warning(_('Some requested columns are missing: {}').format(', '.join(map(str, missing_columns))))
selected_value_columns = [cur for cur in selected_value_columns
if cur in content.columns and cur != 'geometry']
selected_value_positions = [(curkey, col_to_pos[curkey]) for curkey in selected_value_columns]
def normalize_value(value):
if isinstance(value, np.generic):
return value.item()
return value
fallback_specs = [
(colname, lambda x: str(x)) if colname else None,
('NAME', lambda x: str(x)),
('location', lambda x: str(x)),
('Communes', lambda x: str(x)),
('name', lambda x: str(x)),
('MAJ_NIV3T', lambda x: str(x)),
('NATUR_DESC', lambda x: str(x)),
('mun_name_f', lambda x: str(x).replace('[', '').replace(']', '').replace("'", '')),
('mun_name_fr', lambda x: str(x)),
]
name_pos = None
name_formatter = None
for spec in fallback_specs:
if spec is None:
continue
curcol, formatter = spec
if curcol in col_to_pos:
name_pos = col_to_pos[curcol]
name_formatter = formatter
break
zones = []
for row in content.itertuples(index=True, name=None):
idx = row[0]
if name_pos is not None:
name = name_formatter(row[name_pos])
else:
name = str(idx)
poly = row[geometry_pos]
newzone = self._make_zone(name=name, parent = self, fromshapely = poly)
if len(selected_value_positions) > 0 and newzone.nbvectors > 0:
values_to_import = {curkey: normalize_value(row[pos]) for curkey, pos in selected_value_positions}
if share_multipart_values and newzone.nbvectors > 1:
shared_values = dict(values_to_import)
for curvec in newzone.myvectors:
curvec.myprop.set_shared_values(shared_values)
else:
for curvec in newzone.myvectors:
curvec.myprop.add_values(values_to_import)
zones.append(newzone)
self.myzones = zones
pass
[docs]
def export_GeoDataFrame(self, as_multipart: bool = False) -> gpd.GeoDataFrame:
"""
Export to a GeoDataFrame
:param as_multipart: If True, export one row per zone with Multi*
geometries when possible. Default False keeps legacy behavior.
"""
names=[]
geoms=[]
for curzone in self.myzones:
if curzone.nbvectors == 0:
logging.warning(_('Zone {} contains no vector'.format(curzone.myname)))
continue
if as_multipart:
curgeom = self._zone_multipart_geometry_for_export(curzone)
if curgeom is None:
continue
names.append(curzone.myname)
geoms.append(curgeom)
else:
# Legacy behavior: first vector only.
if curzone.nbvectors > 1:
logging.warning(_('Zone {} contains more than one vector -- only the first one will be exported'.format(curzone.myname)))
names.append(curzone.myname)
geoms.append(self._vector_geometry_for_export(curzone.myvectors[0]))
gdf = gpd.GeoDataFrame({'id':names,'geometry':geoms})
gdf.crs = 'EPSG:31370'
return gdf
[docs]
def _vector_geometry_for_export(self, curvect: vectorModel):
"""Return the most suitable shapely geometry for a vector export."""
if curvect.is2D:
if curvect.closed:
return curvect.polygon
return curvect.linestring
if curvect.closed:
return curvect.asshapely_pol3D()
return curvect.asshapely_ls3d()
[docs]
def _zone_multipart_geometry_for_export(self, curzone: zoneModel):
"""Return a Multi* geometry for a zone when possible.
If a zone mixes closed and open vectors, multipart export is ambiguous;
in that case ``None`` is returned and a warning is emitted.
"""
if curzone.nbvectors == 0:
return None
all_closed = all(curvec.closed for curvec in curzone.myvectors)
all_open = all(not curvec.closed for curvec in curzone.myvectors)
if not (all_closed or all_open):
logging.warning(_('Zone {} mixes open and closed vectors -- skipping multipart export for this zone').format(curzone.myname))
return None
geoms = [self._vector_geometry_for_export(curvec) for curvec in curzone.myvectors]
if all_closed:
return MultiPolygon(geoms)
return MultiLineString(geoms)
[docs]
def export_GeoDataFrame_with_values(self,
value_columns: str | list[str] | tuple[str, ...] = 'common',
missing_value=np.nan,
include_vector_meta: bool = True,
as_multipart: bool = False) -> gpd.GeoDataFrame:
"""Export vectors to a temporary GeoDataFrame including property values.
By default, one row is written per vector. With ``as_multipart=True``,
one row is written per zone with a Multi* geometry.
Property keys can be selected with:
- ``'common'``: keys present in every exported vector.
- ``'union'``: all keys found on at least one vector.
- list/tuple: explicit key selection.
:param value_columns: Column-selection mode.
:param missing_value: Fill value for missing keys.
:param include_vector_meta: Include helper metadata columns.
:param as_multipart: Export one Multi* geometry per zone.
:returns: A GeoDataFrame ready to be written with GeoPandas.
"""
records = []
for curzone in self.myzones:
if curzone.nbvectors == 0:
logging.warning(_('Zone {} contains no vector'.format(curzone.myname)))
continue
if as_multipart:
geom = self._zone_multipart_geometry_for_export(curzone)
if geom is None:
continue
values_list = [curvec.myprop.values_dict() for curvec in curzone.myvectors]
common_values = {}
union_keys = set()
for curvals in values_list:
union_keys.update(curvals.keys())
for curkey in union_keys:
present_values = [curvals[curkey] for curvals in values_list if curkey in curvals]
if len(present_values) == len(values_list):
first = present_values[0]
if all(cur == first for cur in present_values[1:]):
common_values[curkey] = first
record = {
'zone_name': curzone.myname,
'vector_name': None,
'vector_index': None,
'vector_count': curzone.nbvectors,
'all_keys': union_keys,
'values': common_values,
'geometry': geom,
}
records.append(record)
else:
for idx_vec, curvect in enumerate(curzone.myvectors):
geom = self._vector_geometry_for_export(curvect)
values = curvect.myprop.values_view()
record = {
'zone_name': curzone.myname,
'vector_name': curvect.myname,
'vector_index': idx_vec,
'vector_count': 1,
'values': values,
'geometry': geom,
}
records.append(record)
if len(records) == 0:
cols = ['geometry']
if include_vector_meta:
cols = ['zone_name', 'vector_name', 'vector_index', 'geometry']
gdf = gpd.GeoDataFrame(columns=cols, geometry='geometry')
gdf.crs = 'EPSG:31370'
return gdf
key_sets_common = [set(rec['values'].keys()) for rec in records]
key_sets_all = [set(rec.get('all_keys', rec['values'].keys())) for rec in records]
if isinstance(value_columns, str):
mode = value_columns.lower()
if mode == 'common':
selected_keys = set.intersection(*key_sets_common) if len(key_sets_common) > 0 else set()
elif mode == 'union':
selected_keys = set.union(*key_sets_all) if len(key_sets_all) > 0 else set()
else:
selected_keys = {value_columns}
else:
selected_keys = set(value_columns)
selected_keys = sorted(selected_keys)
rows = []
for rec in records:
row = {'geometry': rec['geometry']}
if include_vector_meta:
row['zone_name'] = rec['zone_name']
row['vector_name'] = rec['vector_name']
row['vector_index'] = rec['vector_index']
row['vector_count'] = rec['vector_count']
curvals = rec['values']
for curkey in selected_keys:
row[curkey] = curvals.get(curkey, missing_value)
rows.append(row)
gdf = gpd.GeoDataFrame(rows, geometry='geometry')
gdf.crs = 'EPSG:31370'
return gdf
[docs]
def export_to_shapefile(self,
filename:str,
include_values: bool = False,
value_columns: str | list[str] | tuple[str, ...] = 'common',
missing_value=np.nan,
include_vector_meta: bool = True,
as_multipart: bool = False):
"""Export to shapefile.
The first vector of each zone will be exported.
Use ``export_shape`` on the zone object to export all vectors.
:param filename: Output shapefile path.
:param include_values: If True, export per-vector property values.
:param value_columns: Property-key selection mode used when
``include_values=True``.
:param missing_value: Fill value for missing keys.
:param include_vector_meta: Include helper metadata columns.
:param as_multipart: If True, export one row per zone as Multi*.
"""
if include_values:
gdf = self.export_GeoDataFrame_with_values(
value_columns=value_columns,
missing_value=missing_value,
include_vector_meta=include_vector_meta,
as_multipart=as_multipart,
)
else:
gdf = self.export_GeoDataFrame(as_multipart=as_multipart)
gdf.to_file(filename)
[docs]
def export_to_geopackage(self,
filename: str,
layer: str = 'zones',
include_values: bool = True,
value_columns: str | list[str] | tuple[str, ...] = 'common',
missing_value=np.nan,
include_vector_meta: bool = True,
as_multipart: bool = False):
"""Export to GeoPackage using GeoPandas.
:param filename: Output GeoPackage path.
:param layer: Layer name inside the GeoPackage.
:param include_values: If True, export per-vector property values.
:param value_columns: Property-key selection mode used when
``include_values=True``.
:param missing_value: Fill value for missing keys.
:param include_vector_meta: Include helper metadata columns.
:param as_multipart: If True, export one row per zone as Multi*.
"""
if include_values:
gdf = self.export_GeoDataFrame_with_values(
value_columns=value_columns,
missing_value=missing_value,
include_vector_meta=include_vector_meta,
as_multipart=as_multipart,
)
else:
gdf = self.export_GeoDataFrame(as_multipart=as_multipart)
gdf.to_file(filename, layer=layer, driver='GPKG')
[docs]
def export_active_zone_to_shapefile(self, filename:str):
"""Export the active zone to a shapefile.
:param filename: Output shapefile path.
"""
if self.active_zone is None:
logging.warning(_('No active zone'))
return
self.active_zone.export_shape(filename)
[docs]
def import_gdb(self,
fn:str,
bbox:Polygon = None,
layers:list[str] = None,
colname:str = None,
value_columns: str | list[str] | tuple[str, ...] | None = None,
share_multipart_values: bool = False):
"""Import a GDB file using geopandas and Fiona.
:param fn: Path to the GDB file.
:param bbox: Optional bounding box to filter features.
:param layers: List of layer names to import (default: all).
:param colname: Column name used to build zone names.
:param value_columns: Columns to store into vector properties.
:param share_multipart_values: If True, multipart geometries share
the same underlying ``_values`` dictionary across split vectors.
"""
import fiona
if layers is None:
layers = fiona.listlayers(fn)
for curlayer in layers:
content = gpd.read_file(fn, bbox=bbox, layer=curlayer)
if len(content)>1000:
logging.warning(_('Layer {} contains more than 1000 elements -- it may take a while to import'.format(curlayer)))
temp_zones = self._make_zones()
temp_zones.import_GeoDataFrame(content=content,
colname=colname,
value_columns=value_columns,
share_multipart_values=share_multipart_values)
for newzone in temp_zones.myzones:
newzone.parent = self
self.myzones.extend(temp_zones.myzones)
if len(content) > 1000:
logging.info(_('Imported {} elements'.format(len(temp_zones.myzones))))
[docs]
def import_gpkg(self,
fn:str,
bbox:Polygon = None,
layers:list[str] = None,
colname:str = None,
value_columns: str | list[str] | tuple[str, ...] | None = None,
share_multipart_values: bool = False):
"""Import a GeoPackage file using geopandas and Fiona.
:param fn: Path to the GPKG file.
:param bbox: Optional bounding box to filter features.
:param layers: List of layer names to import (default: all).
:param colname: Column name used to build zone names.
:param value_columns: Columns to store into vector properties.
:param share_multipart_values: If True, multipart geometries share
the same underlying ``_values`` dictionary across split vectors.
"""
import fiona
if layers is None:
layers = fiona.listlayers(fn)
for curlayer in layers:
content = gpd.read_file(fn, bbox=bbox, layer=curlayer)
if len(content)>1000:
logging.warning(_('Number of elements in layer {} : {}'.format(curlayer, len(content))))
logging.warning(_('Layer {} contains more than 1000 elements -- it may take a while to import'.format(curlayer)))
temp_zones = self._make_zones()
temp_zones.import_GeoDataFrame(content=content,
colname=colname,
value_columns=value_columns,
share_multipart_values=share_multipart_values)
for newzone in temp_zones.myzones:
newzone.parent = self
self.myzones.extend(temp_zones.myzones)
if len(content) > 1000:
logging.info(_('Imported {} elements'.format(len(temp_zones.myzones))))
[docs]
def colorize_data(self, colors:dict[str:list[int]], filled:bool = False) -> None:
"""Colorize zones based on a dictionary of colours.
Zone names must match the dictionary keys.
:param colors: Mapping ``{zone_name: [R, G, B]}``.
:param filled: Whether polygons should be drawn filled.
"""
std_color = getIfromRGB([10, 10, 10])
for curzone in self.myzones:
if curzone.myname in colors:
curcolor = getIfromRGB(colors[curzone.myname])
else:
curcolor = std_color
for curvect in curzone.myvectors:
curvect.myprop.color = curcolor
curvect.myprop.alpha = 180
curvect.myprop.transparent = True
curvect.myprop.filled = filled and curvect.closed
[docs]
def set_width(self, width:int) -> None:
"""Set line width for all vectors in all zones.
:param width: Width in pixels.
"""
for curzone in self.myzones:
for curvect in curzone.myvectors:
curvect.myprop.width = width
self.prep_listogl()
[docs]
def get_zone(self,keyzone:Union[int, str])->"zoneModel":
"""Retrieve a zone by name or index.
If multiple zones share the same name, only the first is returned.
:param keyzone: Zone index (int) or zone name (str).
:return: The matching zone, or *None* if not found.
"""
if isinstance(keyzone,int):
if keyzone<self.nbzones:
return self.myzones[keyzone]
return None
if isinstance(keyzone,str):
zone_names = [cur.myname for cur in self.myzones]
if keyzone in zone_names:
return self.myzones[zone_names.index(keyzone)]
return None
def __getitem__(self, ndx:Union[int, str, tuple]) -> Union["zoneModel", "vectorModel"]:
"""
Retourne la zone sur base de son nom ou de sa position
:param ndx: Clé ou index de zone -- si tuple, alors (idx_zone, idx_vect) ou (keyzone, keyvect)
"""
if isinstance(ndx,tuple):
idx_zone = ndx[0]
idx_vect = ndx[1]
return self.get_zone(idx_zone)[idx_vect]
else:
return self.get_zone(ndx)
@property
[docs]
def zone_names(self) -> list[str]:
""" Return the list of zone names """
return [cur.myname for cur in self.myzones]
[docs]
def import_dxf(self, fn, imported_elts=['POLYLINE','LWPOLYLINE','LINE']):
"""
Import of a DXF file as a 'Zones'.
The DXF file is read and the elements are stored in zones based on their layers.
The supported elements are POLYLINE, LWPOLYLINE and LINE.
If you want to import other elements, you must upgrade this routine `import_dxf`.
:param fn: name of the DXF file to import
:param imported_elts: list of DXF elements to import. Default is ['POLYLINE','LWPOLYLINE','LINE'].
:return: None
"""
import ezdxf
if not path.exists(fn):
logging.warning(_('File not found !') + ' ' + fn)
return
for elt in imported_elts:
assert elt in ['POLYLINE', 'LWPOLYLINE', 'LINE'], _('Unsupported DXF element: {}').format(elt)
self.is2D = False # we assume it's a 3D DXF
# Lecture du fichier dxf et identification du modelspace
doc = ezdxf.readfile(fn)
msp = doc.modelspace()
# layers = doc.layers
used_layers = {}
notloaded = {}
# Bouclage sur les éléments du DXF pour identifier les layers utiles et ensuite créer les zones adhoc
for e in msp:
if doc.layers.get(e.dxf.layer).is_on():
if e.dxftype() in imported_elts:
if e.dxftype() == "POLYLINE":
if e.dxf.layer not in used_layers.keys():
curlayer = used_layers[e.dxf.layer]={}
else:
curlayer = used_layers[e.dxf.layer]
curlayer[e.dxftype().lower()]=0
elif e.dxftype() == "LWPOLYLINE":
if e.dxf.layer not in used_layers.keys():
curlayer = used_layers[e.dxf.layer]={}
else:
curlayer = used_layers[e.dxf.layer]
curlayer[e.dxftype().lower()]=0
elif e.dxftype() == "LINE": # dans ce cas spécifique, ce sont a priori les lignes composant les croix sur les points levés
if e.dxf.layer not in used_layers.keys():
curlayer = used_layers[e.dxf.layer]={}
else:
curlayer = used_layers[e.dxf.layer]
curlayer[e.dxftype().lower()]=0
else:
if not e.dxftype() in notloaded.keys():
notloaded[e.dxftype()] = 0
notloaded[e.dxftype()] += 1
logging.debug(_('DXF element not supported : ') + e.dxftype())
else:
logging.info(_('Layer {} is off'.format(e.dxf.layer)))
if len(notloaded)>0:
logging.warning(_('Not loaded DXF elements : '))
for curtype in notloaded.keys():
logging.warning(_(' {} : {}'.format(curtype, notloaded[curtype])))
# Création des zones
for curlayer in used_layers.keys():
for curtype in used_layers[curlayer].keys():
curzone = used_layers[curlayer][curtype] = self._make_zone(name = '{} - {}'.format(curlayer,curtype), is2D=self.is2D, parent=self)
self.add_zone(curzone)
# Nouveau bouclage sur les éléments du DXF pour remplissage
nbid=0
for e in msp:
if doc.layers.get(e.dxf.layer).is_on():
if e.dxftype() == "POLYLINE":
nbid+=1
# récupération des coordonnées
verts = [cur.dxf.location.xyz for cur in e.vertices]
curzone = used_layers[e.dxf.layer][e.dxftype().lower()]
curpoly = self._make_vector(is2D=False,name=e.dxf.handle,parentzone=curzone)
curzone.add_vector(curpoly)
for cur in verts:
myvert = wolfvertex(cur[0],cur[1],cur[2])
curpoly.add_vertex(myvert)
elif e.dxftype() == "LWPOLYLINE":
nbid+=1
# récupération des coordonnées
verts = np.array(e.lwpoints.values)
verts = verts.reshape([verts.size // 5,5])[:,:2] # in ezdxf 1.3.5, the lwpoints.values attribute is a np.ndarray [n,5]
# verts = verts.reshape([len(verts) // 5,5])[:,:2] # in ezdxf 1.2.0, the lwpoints.values attribute was flattened
verts = np.column_stack([verts,[e.dxf.elevation]*len(verts)])
curzone = used_layers[e.dxf.layer][e.dxftype().lower()]
curpoly = self._make_vector(is2D=False,name=e.dxf.handle,parentzone=curzone)
curzone.add_vector(curpoly)
for cur in verts:
myvert = wolfvertex(cur[0],cur[1],cur[2])
curpoly.add_vertex(myvert)
elif e.dxftype() == "LINE":
nbid+=1
curzone = used_layers[e.dxf.layer][e.dxftype().lower()]
curpoly = self._make_vector(is2D=False,name=e.dxf.handle,parentzone=curzone)
curzone.add_vector(curpoly)
# récupération des coordonnées
myvert = wolfvertex(e.dxf.start[0],e.dxf.start[1],e.dxf.start[2])
curpoly.add_vertex(myvert)
myvert = wolfvertex(e.dxf.end[0],e.dxf.end[1],e.dxf.end[2])
curpoly.add_vertex(myvert)
logging.info(_('Number of imported elements : ')+str(nbid))
[docs]
def find_nearest_vertex(self, x: float, y: float) -> wolfvertex | None:
"""Find the nearest vertex across all zones."""
distmin = np.inf
minvert = None
for curzone in self.myzones:
curvert = curzone.find_nearest_vertex(x, y)
if curvert is None:
continue
dist = (curvert.x - x) ** 2. + (curvert.y - y) ** 2.
if dist < distmin:
distmin = dist
minvert = curvert
return minvert
[docs]
def find_nearest_vector(self, x:float, y:float) -> "vectorModel":
"""Find the vector whose geometry is closest to a coordinate.
:param x: X-coordinate.
:param y: Y-coordinate.
:return: The nearest vector.
"""
distmin = np.inf
minvec = None
for curzone in self.myzones:
curvect = curzone.find_nearest_vector(x, y)
if curvect is None:
continue
dist = curvect.distance_to_geometry(x, y)
if dist < distmin:
minvec = curvect
distmin = dist
return minvec
[docs]
def find_vector_containing_point(self, x:float, y:float) -> "vectorModel":
"""Find the first vector whose polygon contains the point.
:param x: X-coordinate.
:param y: Y-coordinate.
:return: The containing vector, or *None*.
"""
xy = Point(x, y)
for curzone in self.myzones:
for curvect in curzone.myvectors:
if curvect.polygon.contains(xy):
return curvect
[docs]
def check_plot(self):
"""
L'objet doit être affiché
Fonction principalement utile pour l'objet WolfMapViewer et le GUI
"""
self.plotted = True
[docs]
def uncheck_plot(self, unload=True):
"""Mark the object as not plotted.
:param unload: Reserved for GUI subclasses.
"""
self.plotted = False
[docs]
def save(self):
"""
Sauvegarde sur disque, sans remise en cause du nom de fichier
"""
if self.filename =='':
logging.warning(_('No filename defined'))
return
self.saveas()
[docs]
def saveas(self, filename:str=''):
"""Save to disk, optionally under a new filename.
:param filename: New path; if empty, uses the current ``self.filename``.
"""
filename = str(filename)
if filename!='':
self.filename=filename
if self.filename.endswith('.shp'):
self.export_to_shapefile(self.filename)
else:
if self.filename.endswith('.vecz'):
self.force3D=True #on veut un fichier 3D --> forcage du paramètre
with open(self.filename, 'w') as f:
f.write(f'{self.tx} {self.ty}'+'\n')
f.write(str(self.nbzones)+'\n')
for curzone in self.myzones:
curzone.save(f)
with open(self.filename + '.extra', 'w') as f:
for curzone in self.myzones:
curzone.save_extra(f)
[docs]
def add_zone(self, addedzone:"zoneModel", forceparent=False):
"""Add a zone to the collection.
:param addedzone: The zone to add.
:param forceparent: If True, set this object as the zone's parent.
"""
self.myzones.append(addedzone)
if forceparent:
addedzone.parent = self
[docs]
def create_zone(self, name:str = '') -> "zoneModel":
"""Create a new zone and add it to the collection.
:param name: Name of the new zone.
:return: The newly created zone.
"""
newzone = self._make_zone(name=name, parent=self, is2D=self.is2D)
self.myzones.append(newzone)
return newzone
[docs]
def find_minmax(self, update=False, only_firstlast:bool=False):
"""
Trouve les bornes des vertices pour toutes les zones et tous les vecteurs
:param update : si True, force la MAJ des minmax dans chaque zone; si False, compile les minmax déjà présents
:param only_firstlast : si True, ne prend en compte que le premier et le dernier vertex de chaque vecteur
"""
if self.nbzones > 100:
with ThreadPoolExecutor() as executor:
# zone.find_minmax(update or self._first_find_minmax, only_firstlast)
futures = [executor.submit(zone.find_minmax, update or self._first_find_minmax, only_firstlast) for zone in self.myzones]
wait(futures)
else:
for curzone in self.myzones:
curzone.find_minmax(update or self._first_find_minmax, only_firstlast)
if self.nbzones > 0:
minsx = np.asarray([zone.xmin for zone in self.myzones if zone.xmin!=-99999.])
minsy = np.asarray([zone.ymin for zone in self.myzones if zone.ymin!=-99999.])
maxsx = np.asarray([zone.xmax for zone in self.myzones if zone.xmax!=-99999.])
maxsy = np.asarray([zone.ymax for zone in self.myzones if zone.ymax!=-99999.])
if minsx.size == 0:
self.xmin = 0.
self.ymin = 0.
self.xmax = 1.
self.ymax = 1.
else:
self.xmin = float(minsx.min())
self.xmax = float(maxsx.max())
self.ymin = float(minsy.min())
self.ymax = float(maxsy.max())
[docs]
def plot(self, sx=None, sy=None, xmin=None, ymin=None, xmax=None, ymax=None, size=None):
"""Draw all zones.
:param sx: Scale factor along X.
:param sy: Scale factor along Y.
:param xmin: Minimum X of the viewport.
:param ymin: Minimum Y of the viewport.
:param xmax: Maximum X of the viewport.
:param ymax: Maximum Y of the viewport.
:param size: Reference size for rendering.
"""
for curzone in self.myzones:
curzone.plot(sx=sx, sy=sy, xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax, size=size)
[docs]
def get_bounds(self, force_to_update:bool = True) -> tuple[tuple[float, float], tuple[float, float]]:
"""Return the bounds of all zones.
:param force_to_update: If True, recompute bounds before returning.
:return: ``((xmin, xmax), (ymin, ymax))``.
"""
self.find_minmax(update = force_to_update)
return (self.xmin, self.xmax), (self.ymin, self.ymax)
[docs]
def select_vectors_from_point(self, x:float, y:float, inside=True):
"""Select vectors in each zone from a coordinate.
Fills the ``selected_vectors`` list of each zone.
:param x: X-coordinate.
:param y: Y-coordinate.
:param inside: If True, select by containment; if False, select the
nearest vector.
"""
xmin=1e30
for curzone in self.myzones:
xmin = curzone.select_vectors_from_point(x,y,inside)
[docs]
def create_cs(self, wa) -> "ZonesModel":
"""Create cross-section zones from the active zone.
:param wa: Active WolfArray used to extract cross-section elevations.
:return: A new ZonesModel containing the cross-sections, or *None*.
"""
if wa is None:
logging.warning(_('No active array - Cannot create cross-sections !'))
return None
# The array is forced to be plotted to get the correct coordinates for the cross-section creation.
# It will be reset to its previous state at the end of the function.
old_plotted = wa.plotted
if not old_plotted:
logging.warning(_('The active array is not plotted. It will be temporarily plotted to create the cross-sections.'))
wa.plotted = True
new_zones = self._make_zones()
for curzone in self.myzones:
new_zone = curzone.create_cs(wa)
if new_zone is not None:
new_zones.add_zone(new_zone, forceparent=True)
wa.plotted = old_plotted
if not old_plotted:
logging.warning(_('The active array has been reset to not plotted after cross-section creation.'))
return new_zones
[docs]
def update_from_sz_direction(self, xy1:np.ndarray, xy2:np.ndarray, sz:np.ndarray):
"""Update the active vector from curvilinear s-z values along a direction.
:param xy1: Start point ``[x, y]`` of the direction.
:param xy2: End point ``[x, y]`` of the direction.
:param sz: Array of shape ``(n, 2)`` with columns ``[s, z]``.
"""
if self.active_vector is None:
logging.info(_('No active vector -- Nothing to do'))
return
curv = self.active_vector
# Creation of vertices
if sz.shape[1]==2 and xy1.shape==(2,) and xy2.shape==(2,):
if not np.array_equal(xy1,xy2):
curv.myvertices=[]
dx, dy = xy2[0]-xy1[0], xy2[1]-xy1[1]
norm = np.linalg.norm([dx,dy])
dx, dy = dx/norm, dy/norm
for cur in sz:
x, y = xy1[0] + dx*cur[0], xy1[1] + dy*cur[0]
curv.add_vertex(wolfvertex(x, y, float(cur[1])))
self.find_minmax(True)
[docs]
def get_selected_vectors(self, all=False):
"""Return selected vectors from all zones.
:param all: If True, return all selected vectors as a list; if False,
return only the nearest selected vector.
"""
if all:
mylist=[]
for curzone in self.myzones:
ret = curzone.get_selected_vectors(all)
if ret is not None:
mylist.append(ret)
return mylist
else:
distmin=99999.
vecmin:vectorModel = None
for curzone in self.myzones:
ret = curzone.get_selected_vectors()
if ret is not None:
if (ret[1]<distmin) or (vecmin is None):
distmin= ret[1]
vecmin = ret[0]
return vecmin
[docs]
def unuse(self):
"""
Rend inutilisé l'ensemble des zones
"""
for curzone in self.myzones:
curzone.unuse()
[docs]
def use(self):
"""
Rend utilisé l'ensemble des zones
"""
for curzone in self.myzones:
curzone.use()
[docs]
def _callback_destroy_props(self):
"""Callback invoked when the properties dialog is destroyed."""
self.myprops = None
[docs]
def deepcopy_zones(self, name:str = None) -> "ZonesModel":
"""Return a deep copy of this object.
:param name: Optional name for the copy.
:return: A new ``ZonesModel`` instance.
"""
copied_Zones = self._make_zones(idx=name)
copied_Zones.myzones = list(map(lambda zne: zne.deepcopy_zone(parent= copied_Zones), self.myzones))
# for zne in self.myzones:
# new_zne = zne.deepcopy_zone(parent= copied_Zones)
# copied_Zones.add_zone(new_zne,forceparent=True)
copied_Zones.find_minmax(True)
return copied_Zones
[docs]
def deepcopy(self, name:str = None) -> "ZonesModel":
"""Return a deep copy of this object (alias for :meth:`deepcopy_zones`).
:param name: Optional name for the copy.
:return: A new ``ZonesModel`` instance.
"""
return self.deepcopy_zones(name=name)
[docs]
def keep_only_trace(self):
"""
Keep only the trace (first and last vertices) of all vectors
for all zones in this Zones.
"""
for curzone in self.myzones:
curzone.keep_only_trace()
[docs]
def deepcopy_only_trace(self, name:str = None) -> "ZonesModel":
"""Return a deep copy keeping only trace vertices (first and last).
:param name: Optional name for the copy.
:return: A new ``ZonesModel`` instance with trace-only vectors.
"""
copied_Zones = self.deepcopy(name=name)
copied_Zones.keep_only_trace()
return copied_Zones
[docs]
def extend_along_trace(self, distance:float | str):
"""
Extend all vectors along their trace by the given distance.
:param distance: Distance to extend (float) or 'auto' to extend
by 10% of the vector length.
"""
for curzone in self.myzones:
curzone.extend_along_trace(distance)
# ====================================================================
# GridModel
# ====================================================================
[docs]
class GridModel(ZonesModel):
"""Pure-data Grid: one zone with gridx, gridy and contour vectors.
GUI subclass ``Grid`` in :mod:`wolfhece.PyVertexvectors` inherits
from this and adds wx/OpenGL support via ``Zones``.
"""
def __init__(self,
size: float = 1000.,
ox: float = 0.,
oy: float = 0.,
ex: float = 1000.,
ey: float = 1000.) -> None:
"""Initialise a grid with one zone containing gridx, gridy and contour vectors.
:param size: Grid cell size in world units.
:param ox: X-coordinate of the grid origin (lower-left).
:param oy: Y-coordinate of the grid origin (lower-left).
:param ex: X-coordinate of the grid extent (upper-right).
:param ey: Y-coordinate of the grid extent (upper-right).
"""
super().__init__(ox=ox, oy=oy)
mygrid = self._make_zone(name='grid', parent=self)
self.add_zone(mygrid)
gridx = self._make_vector(name='gridx')
gridy = self._make_vector(name='gridy')
contour = self._make_vector(name='contour')
mygrid.add_vector(gridx)
mygrid.add_vector(gridy)
mygrid.add_vector(contour)
self.creategrid(size, ox, oy, ex, ey)
[docs]
def creategrid(self,
size: float = 100.,
ox: float = 0.,
oy: float = 0.,
ex: float = 1000.,
ey: float = 1000.):
"""(Re)create the grid lines and contour.
Resets the gridx, gridy and contour vectors and rebuilds them
as a regular rectangular grid.
:param size: Grid cell size in world units.
:param ox: X-coordinate of the grid origin (lower-left).
:param oy: Y-coordinate of the grid origin (lower-left).
:param ex: X-coordinate of the grid extent (upper-right).
:param ey: Y-coordinate of the grid extent (upper-right).
"""
mygridx = self.myzones[0].myvectors[0]
mygridy = self.myzones[0].myvectors[1]
contour = self.myzones[0].myvectors[2]
mygridx.reset()
mygridy.reset()
contour.reset()
locox = int(ox / size) * size
locoy = int(oy / size) * size
locex = (np.ceil(ex / size)) * size
locey = (np.ceil(ey / size)) * size
nbx = int(np.ceil((locex - locox) / size))
nby = int(np.ceil((locey - locoy) / size))
dx = locex - locox
dy = locey - locoy
# grillage vertical
xloc = locox
yloc = locoy
for i in range(nbx):
newvert = wolfvertex(xloc, yloc)
mygridx.add_vertex(newvert)
yloc += dy
newvert = wolfvertex(xloc, yloc)
mygridx.add_vertex(newvert)
xloc += size
dy = -dy
# grillage horizontal
xloc = locox
yloc = locoy
for i in range(nby):
newvert = wolfvertex(xloc, yloc)
mygridy.add_vertex(newvert)
xloc += dx
newvert = wolfvertex(xloc, yloc)
mygridy.add_vertex(newvert)
yloc += size
dx = -dx
newvert = wolfvertex(locox, locoy)
contour.add_vertex(newvert)
newvert = wolfvertex(locex, locoy)
contour.add_vertex(newvert)
newvert = wolfvertex(locex, locey)
contour.add_vertex(newvert)
newvert = wolfvertex(locox, locey)
contour.add_vertex(newvert)
self.find_minmax(True)