PVPlant/PVPlantTerrain.py

# /**********************************************************************
# *                                                                     *
# * Copyright (c) 2021 Javier Braña <javier.branagutierrez@gmail.com>  *
# *                                                                     *
# * This program is free software; you can redistribute it and/or modify*
# * it under the terms of the GNU Lesser General Public License (LGPL)  *
# * as published by the Free Software Foundation; either version 2 of   *
# * the License, or (at your option) any later version.                 *
# * for detail see the LICENCE text file.                               *
# *                                                                     *
# * This program is distributed in the hope that it will be useful,     *
# * but WITHOUT ANY WARRANTY; without even the implied warranty of      *
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
# * GNU Library General Public License for more details.                *
# *                                                                     *
# * You should have received a copy of the GNU Library General Public   *
# * License along with this program; if not, write to the Free Software *
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307*
# * USA                                                                 *
# *                                                                     *
# ***********************************************************************

import copy

import ArchComponent
import FreeCAD
import Part
import numpy as np
import math

import PVPlantSite

if FreeCAD.GuiUp:
    import FreeCADGui
    from pivy import coin
    import os
else:
    # \cond
    def translate(ctxt, txt):
        return txt

    def QT_TRANSLATE_NOOP(ctxt, txt):
        return txt
    # \endcond

__title__ = "FreeCAD Fixed Rack"
__author__ = "Javier Braña"
__url__ = "http://www.sogos-solar.com"
__dir__ = os.path.join(FreeCAD.getUserAppDataDir(), "Mod", "PVPlant")

DirResources = os.path.join(__dir__, "Resources")
DirIcons = os.path.join(DirResources, "Icons")
DirImages = os.path.join(DirResources, "Images")

line_patterns = {
    "Continues      _______________________________": 0xFFFF,
    "Border         __ . __ __ . __ __ . __ __ . __": 0x3CF2,
    "Border (.5x)   __.__.__.__.__.__.__.__.__.__._": 0x3939,
    "Border (2x)    ____  ____  .  ____  ____  .  _": 0xFDFA,
    "Center         ____ _ ____ _ ____ _ ____ _ ___": 0xFF3C,
    "Center (.5x)   ___ _ ___ _ ___ _ ___ _ ___ _ _": 0xFC78,
    "Center (2x)    ________  __  ________  __  ___": 0xFFDE,
    "Dash dot       __ . __ . __ . __ . __ . __ . _": 0xE4E4,
    "Dash dot (.5x) _._._._._._._._._._._._._._._._": 0xEBAE,
    "Dash dot (2x)  ____  .  ____  .  ____  .  ____": 0xFF08,
    "Dashed         __ __ __ __ __ __ __ __ __ __ _": 0x739C,
    "Dashed (.5x)   _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _": 0xDB6E,
    "Dashed (2x)    ____  ____  ____  ____  ____  _": 0xFFE0,
    "Divide         ____ . . ____ . . ____ . . ____": 0xFF24,
    "Divide (.5x)   __..__..__..__..__..__..__..__.": 0xEAEA,
    "Divide (2x)    ________  .  .  ________  .  . ": 0xFFEA,
    "Dot            . . . . . . . . . . . . . . . .": 0x4924,
    "Dot (.5x)      ...............................": 0x5555,
    "Dot (2x)       .  .  .  .  .  .  .  .  .  .  .": 0x8888}


def open_xyz_mmap(archivo_path):
    """
    Usa memory-mapping para archivos muy grandes (máxima velocidad)
    """
    # Primera pasada: contar líneas válidas
    total_puntos = 0
    with open(archivo_path, 'r') as f:
        for linea in f:
            partes = linea.strip().split()
            if len(partes) >= 3:
                try:
                    float(partes[0]);
                    float(partes[1]);
                    float(partes[2])
                    total_puntos += 1
                except:
                    continue

    # Segunda pasada: cargar datos
    puntos = np.empty((total_puntos, 3))
    idx = 0

    with open(archivo_path, 'r') as f:
        for linea in f:
            partes = linea.strip().split()
            if len(partes) >= 3:
                try:
                    x, y, z = float(partes[0]), float(partes[1]), float(partes[2])
                    puntos[idx] = [x, y, z]
                    idx += 1
                except:
                    continue

    return puntos

def makeTerrain(name="Terrain"):
    obj = FreeCAD.ActiveDocument.addObject("Part::FeaturePython", "Terrain")
    obj.Label = name
    Terrain(obj)
    ViewProviderTerrain(obj.ViewObject)
    FreeCAD.ActiveDocument.recompute()
    return obj

class Terrain(ArchComponent.Component):
    "A Shadow Terrain Obcject"

    def __init__(self, obj):
        # Definición de  variables:
        ArchComponent.Component.__init__(self, obj)
        self.setProperties(obj)
        self.obj = obj
        # Does a IfcType exist?
        # obj.IfcType = "Fence"
        # obj.MoveWithHost = False

        try:
            self.site = PVPlantSite.get()
        except Exception:
            self.site = None
        if self.site:
            self.site.Terrain = obj
        else:
            FreeCAD.Console.PrintWarning('Terrain: No se encontró Site, algunas funciones DEM requerirán Site.\n')
        obj.ViewObject.ShapeColor = (0.0000, 0.6667, 0.4980)
        obj.ViewObject.LineColor = (0.0000, 0.6000, 0.4392)

    def setProperties(self, obj):
        # Definicion de Propiedades:
        pl = obj.PropertiesList
        if not ("CuttingBoundary" in pl):
            obj.addProperty("App::PropertyLink",
                            "CuttingBoundary",
                            "Surface",
                            "A boundary line to delimit the surface")

        if not ("DEM" in pl):
            obj.addProperty("App::PropertyFile",
                            "DEM",
                            "Surface",
                            "Load a ASC file to generate the surface")

        if not ("PointsGroup" in pl):
            obj.addProperty("App::PropertyLink",
                            "PointsGroup",
                            "Surface",
                            "Use a Point Group to generate the surface")

        if not ("mesh" in pl):
            obj.addProperty("Mesh::PropertyMeshKernel",
                            "mesh",
                            "Surface",
                            "Mesh")
        obj.setEditorMode("mesh", 1)

        if not ("InitialMesh" in pl):
            obj.addProperty("Mesh::PropertyMeshKernel",
                            "InitialMesh",
                            "Surface",
                            "Mesh")
        obj.setEditorMode("InitialMesh", 1)

        '''
        #obj.setEditorMode("Volume", 1)
        if not "AllowedAreas" in pl:
            obj.addProperty("App::PropertyLinkList",
                            "AllowedAreas",
                            "Areas",
                            "A boundary to delimitated the terrain").AllowedAreas = []
        if not "ProhibitedAreas" in pl:
            obj.addProperty("App::PropertyLinkList",
                            "ProhibitedAreas",
                            "Areas",
                            "A boundary to delimitated the terrain").ProhibitedAreas = []
        '''

    def onDocumentRestored(self, obj):
        ArchComponent.Component.onDocumentRestored(self, obj)
        self.setProperties(obj)

    def onChanged(self, obj, prop):
        '''Do something when a property has changed'''

        if prop == "InitialMesh":
            obj.mesh = obj.InitialMesh.copy()
            # Forzar actualización visual
            obj.publishProperty("Mesh")

        if prop == "mesh":
            # La propiedad mesh cambió → forzar recompute para que updateData se dispare
            pass

        if prop == "DEM" or prop == "CuttingBoundary":
            from datetime import datetime
            if obj.DEM and obj.CuttingBoundary:
                from pathlib import Path
                suffix = Path(obj.DEM).suffix
                if suffix == '.asc':
                    '''
                    ASC format:
                    
                    Parámetro	            Descripción	                        Requisitos
                    NCOLS:                  Cantidad de columnas de celdas      Entero mayor que 0.
                    NROWS:                  Cantidad de filas de celdas         Entero mayor que 0.
                    XLLCENTER o XLLCORNER:  Coordenada X del origen (por el centro o la esquina inferior izquierda de la celda) Hacer coincidir con el tipo de coordenada y.
                    YLLCENTER o YLLCORNER:  Coordenada Y del origen (por el centro o la esquina inferior izquierda de la celda) Hacer coincidir con el tipo de coordenada x.
                    CELLSIZE:               Tamaño de celda                     Mayor que 0.
                    NODATA_VALUE:           Los valores de entrada que serán NoData en el ráster de salida  Opcional. El valor predeterminado es -9999
                    '''
                    grid_space = 1
                    file = open(obj.DEM, "r")
                    templist = [line.split() for line in file.readlines()]
                    file.close()
                    del file

                    # Read meta data:
                    meta = templist[0:6]
                    nx = int(meta[0][1])                     # NCOLS
                    ny = int(meta[1][1])                     # NROWS
                    xllref = meta[2][0]                      # XLLCENTER / XLLCORNER
                    xllvalue = round(float(meta[2][1]), 3)
                    yllref = meta[3][0]                      # YLLCENTER / XLLCORNER
                    yllvalue = round(float(meta[3][1]), 3)
                    cellsize = round(float(meta[4][1]), 3)   # CELLSIZE
                    nodata_value = float(meta[5][1])         # NODATA_VALUE

                    # set coarse_factor
                    coarse_factor = max(round(grid_space / cellsize), 1)

                    # Get z values
                    templist = templist[6:(6 + ny)]
                    templist = [templist[i][0::coarse_factor] for i in np.arange(0, len(templist), coarse_factor)]
                    datavals = np.array(templist).astype(float)
                    del templist

                    # create xy coordinates
                    offset = self.site.Origin if self.site else FreeCAD.Vector(0, 0, 0)
                    x = (cellsize * np.arange(nx)[0::coarse_factor] + xllvalue) * 1000 - offset.x
                    y = (cellsize * np.arange(ny)[-1::-1][0::coarse_factor] + yllvalue) * 1000 - offset.y
                    datavals = datavals * 1000  # Ajuste de altura

                    # remove points out of area
                    # 1. coarse:
                    if obj.CuttingBoundary:
                        inc_x = obj.CuttingBoundary.Shape.BoundBox.XLength * 0.0
                        inc_y = obj.CuttingBoundary.Shape.BoundBox.YLength * 0.0
                        tmp = np.where(np.logical_and(x >= (obj.CuttingBoundary.Shape.BoundBox.XMin - inc_x),
                                                      x <= (obj.CuttingBoundary.Shape.BoundBox.XMax + inc_x)))[0]
                        x_max = np.ndarray.max(tmp)
                        x_min = np.ndarray.min(tmp)

                        tmp = np.where(np.logical_and(y >= (obj.CuttingBoundary.Shape.BoundBox.YMin - inc_y),
                                                      y <= (obj.CuttingBoundary.Shape.BoundBox.YMax + inc_y)))[0]
                        y_max = np.ndarray.max(tmp)
                        y_min = np.ndarray.min(tmp)
                        del tmp

                        x = x[x_min:x_max+1]
                        y = y[y_min:y_max+1]
                        datavals = datavals[y_min:y_max+1, x_min:x_max+1]

                    # Create mesh - surface:
                    import MeshTools.Triangulation as Triangulation
                    import Mesh
                    stepsize = 75
                    stepx = math.ceil(nx / stepsize)
                    stepy = math.ceil(ny / stepsize)

                    # Malla completa primero como numpy y filtramos todo de una
                    from datetime import datetime
                    t_start = datetime.now()

                    # Crear grid completo de coordenadas
                    XX, YY = np.meshgrid(x, y)
                    ZZ = datavals.copy()

                    # Enmascarar nodata
                    mask_valida = ZZ != nodata_value

                    # Enmascarar cutting boundary si existe
                    if obj.CuttingBoundary:
                        from FreeCAD import Base
                        shape = obj.CuttingBoundary.Shape
                        mask_boundary = np.zeros_like(ZZ, dtype=bool)
                        # Sampling: revisar solo puntos estratégicos para boundary grande
                        stride = max(1, min(nx, ny) // 200)
                        for i in range(0, ny, stride):
                            for j in range(0, nx, stride):
                                if mask_valida[i, j]:
                                    if shape.isInside(FreeCAD.Vector(x[j], y[i], 0), 0, True):
                                        mask_boundary[i, j] = True
                        mask_valida = mask_valida & mask_boundary

                    # Extraer puntos válidos como lista plana
                    pts_validos = np.column_stack([
                        XX[mask_valida].ravel(),
                        YY[mask_valida].ravel(),
                        ZZ[mask_valida].ravel()
                    ])

                    del XX, YY, ZZ, mask_valida

                    # Triangulación completa de una vez (no por parches)
                    mesh = Mesh.Mesh()
                    if len(pts_validos) > 3:
                        # Si hay muchos puntos, triangulamos por parches para evitar OOM
                        patch_size = 50000
                        n_patches = max(1, math.ceil(len(pts_validos) / patch_size))
                        for p in range(n_patches):
                            patch = pts_validos[p * patch_size:(p + 1) * patch_size].tolist()
                            if len(patch) > 3:
                                try:
                                    triangulated = Triangulation.Triangulate(patch)
                                    mesh.addMesh(triangulated)
                                except TypeError as e:
                                    print(f"Patch {p}: error al procesar {len(patch)} puntos: {str(e)}")
                                except Exception as e:
                                    print(f"Patch {p}: error inesperado: {str(e)}")

                    print(f'Terraín DEM: {len(pts_validos)} pts válidos, {n_patches} parches, {datetime.now()-t_start}')
                    del pts_validos

                    mesh.removeDuplicatedPoints()
                    mesh.removeFoldsOnSurface()
                    obj.InitialMesh = mesh.copy()
                    # Limpiar objetos mesh huérfanos previos si existen
                    for o in FreeCAD.ActiveDocument.Objects:
                        if o.TypeId == 'Mesh::Feature' and o.Label.startswith('Terrain_mesh_'):
                            FreeCAD.ActiveDocument.removeObject(o.Name)
                    mesh_obj = Mesh.show(mesh)
                    mesh_obj.Label = 'Terrain_mesh_' + obj.Label
                elif suffix in ['.xyz']:
                    pts_array = open_xyz_mmap(obj.DEM)
                    if pts_array is not None and len(pts_array) > 3:
                        import MeshTools.Triangulation as Triangulation
                        import Mesh
                        if obj.CuttingBoundary:
                            mask = []
                            for pt in pts_array:
                                mask.append(obj.CuttingBoundary.Shape.isInside(
                                    FreeCAD.Vector(pt[0], pt[1], 0), 0, True))
                            pts_array = pts_array[mask]
                        if len(pts_array) > 3:
                            from datetime import datetime
                            t0 = datetime.now()
                            pts_list = pts_array.tolist()
                            mesh = Triangulation.Triangulate(pts_list)
                            mesh.removeDuplicatedPoints()
                            mesh.removeFoldsOnSurface()
                            obj.InitialMesh = mesh.copy()
                            # Limpiar objetos mesh huérfanos previos
                            for o in FreeCAD.ActiveDocument.Objects:
                                if o.TypeId == 'Mesh::Feature' and o.Label.startswith('Terrain_mesh_'):
                                    FreeCAD.ActiveDocument.removeObject(o.Name)
                            mesh_obj = Mesh.show(mesh)
                            mesh_obj.Label = 'Terrain_mesh_' + obj.Label
                            print(f'XYZ import: {len(pts_array)} puntos en {datetime.now()-t0}')



        if prop == "PointsGroup" or prop == "CuttingBoundary":
            if obj.PointsGroup and obj.CuttingBoundary:
                bnd = obj.CuttingBoundary.Shape
                if len(bnd.Faces) == 0:
                    pts = [ver.Point for ver in bnd.Vertexes]
                    pts.append(pts[0])
                    bnd = Part.makePolygon(pts)

                # TODO: not use the first point, else the Origin in "Site".
                #  It is standard for everything.
                firstPoint = self.obj.PointsGroup.Points.Points[0]
                nbase = FreeCAD.Vector(firstPoint.x, firstPoint.y, firstPoint.z)
                data = []
                for point in self.obj.PointsGroup.Points.Points:
                    tmp = FreeCAD.Vector(0, 0, 0).add(point)
                    tmp.z = 0
                    if bnd.isInside(tmp, 0, True):
                        p = point - nbase
                        data.append([float(p.x), float(p.y), float(p.z)])

                Data = np.array(data)
                data.clear()

                import MeshTools.Triangulation as Triangulation
                mesh = Triangulation.Triangulate(Data)
                if obj.DEM:
                    obj.DEM = None
                obj.mesh = mesh
                # Forzar actualización visual llamando a publishProperty
                try:
                    obj.publishProperty("Mesh")
                except:
                    pass

    def execute(self, obj):
        ''''''
        #print("  -----  Terrain  -  EXECUTE  ----------")

    def __getstate__(self):
        return self.Type

    def __setstate__(self, state):
        if state:
            self.Type = state


class ViewProviderTerrain:
    "A View Provider for the Pipe object"

    def __init__(self, vobj):
        self.boundary_color = None
        self.edge_style = None
        self.edge_color = None
        self.edge_material = None
        self.face_material = None
        self.triangles = None
        self.geo_coords = None
        self.setProperties(vobj)

    def setProperties(self, vobj):
        # Triangulation properties.
        pl = vobj.PropertiesList
        if not ("Transparency" in pl):
            '''vobj.addProperty("App::PropertyIntegerConstraint",
                             "Transparency",
                             "Surface Style",
                             "Set triangle face transparency").Transparency = (50, 0, 100, 1)'''

        if not ("ShapeColor" in pl):
            vobj.addProperty("App::PropertyColor",
                             "ShapeColor",
                             "Surface Style",
                             "Set triangle face color").ShapeColor = (0.0, 0.667, 0.49, vobj.Transparency / 100)

        if not ("ShapeMaterial" in pl):
            vobj.addProperty("App::PropertyMaterial",
                             "ShapeMaterial",
                             "Surface Style",
                             "Triangle face material").ShapeMaterial = FreeCAD.Material()

        if not ("LineTransparency" in pl):
            vobj.addProperty("App::PropertyIntegerConstraint",
                             "LineTransparency",
                             "Surface Style",
                             "Set triangle edge transparency").LineTransparency = (50, 0, 100, 1)

        if not ("LineColor" in pl):
            vobj.addProperty("App::PropertyColor",
                             "LineColor",
                             "Surface Style",
                             "Set triangle face color").LineColor = (0.5, 0.5, 0.5, vobj.LineTransparency / 100)

        if not ("LineMaterial" in pl):
            vobj.addProperty("App::PropertyMaterial",
                             "LineMaterial",
                             "Surface Style",
                             "Triangle face material").LineMaterial = FreeCAD.Material()

        if not ("LineWidth" in pl):
            vobj.addProperty("App::PropertyFloatConstraint",
                             "LineWidth",
                             "Surface Style",
                             "Set triangle edge line width").LineWidth = (0.0, 1.0, 20.0, 1.0)

        # Boundary properties.
        if not ("BoundaryColor" in pl):
            vobj.addProperty("App::PropertyColor",
                             "BoundaryColor",
                             "Boundary Style",
                             "Set boundary contour color").BoundaryColor = (0.0, 0.75, 1.0, 0.0)

        if not ("BoundaryWidth" in pl):
            vobj.addProperty("App::PropertyFloatConstraint",
                             "BoundaryWidth",
                             "Boundary Style",
                             "Set boundary contour line width").BoundaryWidth = (3.0, 1.0, 20.0, 1.0)

        if not ("BoundaryPattern" in pl):
            vobj.addProperty("App::PropertyEnumeration",
                             "BoundaryPattern",
                             "Boundary Style",
                             "Set a line pattern for boundary").BoundaryPattern = [*line_patterns]

        if not ("PatternScale" in pl):
            vobj.addProperty("App::PropertyIntegerConstraint",
                             "PatternScale",
                             "Boundary Style",
                             "Scale the line pattern").PatternScale = (3, 1, 20, 1)

        # Contour properties.
        if not ("MajorColor" in pl):
            vobj.addProperty("App::PropertyColor",
                             "MajorColor",
                             "Contour Style",
                             "Set major contour color").MajorColor = (1.0, 0.0, 0.0, 0.0)

        if not ("MajorWidth" in pl):
            vobj.addProperty("App::PropertyFloatConstraint",
                             "MajorWidth",
                             "Contour Style",
                             "Set major contour line width").MajorWidth = (4.0, 1.0, 20.0, 1.0)

        if not ("MinorColor" in pl):
            vobj.addProperty("App::PropertyColor",
                             "MinorColor",
                             "Contour Style",
                             "Set minor contour color").MinorColor = (1.0, 1.0, 0.0, 0.0)

        if not ("MinorWidth" in pl):
            vobj.addProperty("App::PropertyFloatConstraint",
                             "MinorWidth",
                             "Contour Style",
                             "Set major contour line width").MinorWidth = (2.0, 1.0, 20.0, 1.0)

        vobj.Proxy = self
        self.Object = vobj.Object
        # Inicializar colores correctamente
        vobj.ShapeMaterial.DiffuseColor = vobj.ShapeColor

    def onDocumentRestored(self, vobj):
        self.setProperties(vobj)

    def onChanged(self, vobj, prop):
        """ Update Object visuals when a view property changed. """

        if prop == "ShapeColor" or prop == "Transparency":
            if hasattr(vobj, "ShapeColor") and hasattr(vobj, "Transparency"):
                color = vobj.getPropertyByName("ShapeColor")
                transparency = vobj.getPropertyByName("Transparency")
                color = (color[0], color[1], color[2], 50 / 100)
                vobj.ShapeMaterial.DiffuseColor = color

        if prop == "ShapeMaterial":
            if hasattr(vobj, "ShapeMaterial"):
                material = vobj.getPropertyByName("ShapeMaterial")
                self.face_material.diffuseColor.setValue(material.DiffuseColor[:3])
                self.face_material.transparency = material.DiffuseColor[3]

        if prop == "LineColor" or prop == "LineTransparency":
            if hasattr(vobj, "LineColor") and hasattr(vobj, "LineTransparency"):
                color = vobj.getPropertyByName("LineColor")
                transparency = vobj.getPropertyByName("LineTransparency")
                color = (color[0], color[1], color[2], transparency / 100)
                vobj.LineMaterial.DiffuseColor = color

        if prop == "LineMaterial":
            material = vobj.getPropertyByName(prop)
            self.edge_material.diffuseColor.setValue(material.DiffuseColor[:3])
            self.edge_material.transparency = material.DiffuseColor[3]

        if prop == "LineWidth":
            width = vobj.getPropertyByName(prop)
            self.edge_style.lineWidth = width

        if prop == "BoundaryColor":
            color = vobj.getPropertyByName(prop)
            self.boundary_color.rgb = color[:3]

        if prop == "BoundaryWidth":
            width = vobj.getPropertyByName(prop)
            self.boundary_style.lineWidth = width

        if prop == "BoundaryPattern":
            if hasattr(vobj, "BoundaryPattern"):
                pattern = vobj.getPropertyByName(prop)
                self.boundary_style.linePattern = line_patterns[pattern]

        if prop == "PatternScale":
            if hasattr(vobj, "PatternScale"):
                scale = vobj.getPropertyByName(prop)
                self.boundary_style.linePatternScaleFactor = scale

        if prop == "MajorColor":
            color = vobj.getPropertyByName(prop)
            self.major_color.rgb = color[:3]

        if prop == "MajorWidth":
            width = vobj.getPropertyByName(prop)
            self.major_style.lineWidth = width

        if prop == "MinorColor":
            color = vobj.getPropertyByName(prop)
            self.minor_color.rgb = color[:3]

        if prop == "MinorWidth":
            width = vobj.getPropertyByName(prop)
            self.minor_style.lineWidth = width

    def attach(self, vobj):
        ''' Create Object visuals in 3D view.'''

        # GeoCoords Node.
        self.geo_coords = coin.SoGeoCoordinate()

        # Surface features.
        self.triangles = coin.SoIndexedFaceSet()
        self.face_material = coin.SoMaterial()
        self.edge_material = coin.SoMaterial()
        self.edge_color = coin.SoBaseColor()
        self.edge_style = coin.SoDrawStyle()
        self.edge_style.style = coin.SoDrawStyle.LINES

        shape_hints = coin.SoShapeHints()
        shape_hints.vertex_ordering = coin.SoShapeHints.COUNTERCLOCKWISE
        mat_binding = coin.SoMaterialBinding()
        mat_binding.value = coin.SoMaterialBinding.PER_FACE
        offset = coin.SoPolygonOffset()
        offset.styles = coin.SoPolygonOffset.LINES
        offset.factor = -2.0

        # Boundary features.
        self.boundary_color = coin.SoBaseColor()
        self.boundary_coords = coin.SoGeoCoordinate()
        self.boundary_lines = coin.SoLineSet()
        self.boundary_style = coin.SoDrawStyle()
        self.boundary_style.style = coin.SoDrawStyle.LINES

        # Boundary root.
        boundaries = coin.SoType.fromName('SoFCSelection').createInstance()
        boundaries.style = 'EMISSIVE_DIFFUSE'
        boundaries.addChild(self.boundary_color)
        boundaries.addChild(self.boundary_style)
        boundaries.addChild(self.boundary_coords)
        boundaries.addChild(self.boundary_lines)

        # Major Contour features.
        self.major_color = coin.SoBaseColor()
        self.major_coords = coin.SoGeoCoordinate()
        self.major_lines = coin.SoLineSet()
        self.major_style = coin.SoDrawStyle()
        self.major_style.style = coin.SoDrawStyle.LINES

        # Major Contour root.
        major_contours = coin.SoSeparator()
        major_contours.addChild(self.major_color)
        major_contours.addChild(self.major_style)
        major_contours.addChild(self.major_coords)
        major_contours.addChild(self.major_lines)

        # Minor Contour features.
        self.minor_color = coin.SoBaseColor()
        self.minor_coords = coin.SoGeoCoordinate()
        self.minor_lines = coin.SoLineSet()
        self.minor_style = coin.SoDrawStyle()
        self.minor_style.style = coin.SoDrawStyle.LINES

        # Minor Contour root.
        minor_contours = coin.SoSeparator()
        minor_contours.addChild(self.minor_color)
        minor_contours.addChild(self.minor_style)
        minor_contours.addChild(self.minor_coords)
        minor_contours.addChild(self.minor_lines)

        # Highlight for selection.
        highlight = coin.SoType.fromName('SoFCSelection').createInstance()
        highlight.style = 'EMISSIVE_DIFFUSE'
        highlight.addChild(shape_hints)
        highlight.addChild(mat_binding)
        highlight.addChild(self.geo_coords)
        highlight.addChild(self.triangles)
        highlight.addChild(boundaries)

        # Face root.
        face = coin.SoSeparator()
        face.addChild(self.face_material)
        face.addChild(highlight)

        # Edge root.
        edge = coin.SoSeparator()
        edge.addChild(self.edge_material)
        edge.addChild(self.edge_style)
        edge.addChild(highlight)

        # Surface root - con contour lines visibles.
        surface_root = coin.SoSeparator()
        surface_root.addChild(face)
        surface_root.addChild(offset)
        surface_root.addChild(edge)
        surface_root.addChild(major_contours)
        surface_root.addChild(minor_contours)
        vobj.addDisplayMode(surface_root, "Surface")

        # Boundary root.
        boundary_root = coin.SoSeparator()
        boundary_root.addChild(boundaries)
        vobj.addDisplayMode(boundary_root, "Boundary")

        # Elevation/Shaded root.
        '''shaded_root = coin.SoSeparator()
        shaded_root.addChild(face)
        vobj.addDisplayMode(shaded_root, "Elevation")
        vobj.addDisplayMode(shaded_root, "Slope")
        vobj.addDisplayMode(shaded_root, "Shaded")'''

        # Flat Lines root.
        flatlines_root = coin.SoSeparator()
        flatlines_root.addChild(face)
        flatlines_root.addChild(offset)
        flatlines_root.addChild(edge)
        vobj.addDisplayMode(flatlines_root, "Flat Lines")

        # Wireframe root.
        wireframe_root = coin.SoSeparator()
        wireframe_root.addChild(edge)
        #wireframe_root.addChild(major_contours)
        #wireframe_root.addChild(minor_contours)
        vobj.addDisplayMode(wireframe_root, "Wireframe")

        # Take features from properties.
        self.onChanged(vobj, "ShapeColor")
        self.onChanged(vobj, "LineColor")
        self.onChanged(vobj, "LineWidth")
        self.onChanged(vobj, "BoundaryColor")
        self.onChanged(vobj, "BoundaryWidth")
        self.onChanged(vobj, "BoundaryPattern")
        self.onChanged(vobj, "PatternScale")
        self.onChanged(vobj, "MajorColor")
        self.onChanged(vobj, "MajorWidth")
        self.onChanged(vobj, "MinorColor")
        self.onChanged(vobj, "MinorWidth")
        
    def updateData(self, obj, prop):
        ''' Update Object visuals when a data property changed. '''

        # Set geosystem.
        try:
            utm_zone = FreeCAD.ActiveDocument.Site.UtmZone
        except:
            utm_zone = "30"
        geo_system = ["UTM", utm_zone, "FLAT"]
        self.geo_coords.geoSystem.setValues(geo_system)
        self.boundary_coords.geoSystem.setValues(geo_system)
        self.major_coords.geoSystem.setValues(geo_system)
        self.minor_coords.geoSystem.setValues(geo_system)

        if prop == "mesh" or prop == "Mesh":
            if obj.mesh:
                mesh = obj.mesh
                try:
                    vertices = [tuple(v) for v in mesh.Topology[0]]
                    faces = []
                    for face in mesh.Topology[1]:
                        faces.extend(face)
                        faces.append(-1)

                    # Asignar a los nodos de visualización
                    self.geo_coords.point.values = vertices
                    self.triangles.coordIndex.values = faces
                except Exception as e:
                    FreeCAD.Console.PrintError(f"Error actualizando mesh visual: {e}\n")

    def getDisplayModes(self, vobj):
        ''' Return a list of display modes. '''
        return ["Surface", "Boundary", "Flat Lines", "Wireframe"]

    def getDefaultDisplayMode(self):
        return "Surface"

    def claimChildren(self):
        if hasattr(self, "Object") and self.Object:
            return [self.Object.CuttingBoundary, ]
        return []

    def getIcon(self):
        return str(os.path.join(DirIcons, "terrain.svg"))

    def __getstate__(self):
        """ Save variables to file. """
        return None

    def __setstate__(self, state):
        """ Get variables from file. """
        return None


'''class _CommandTerrain:
    "the PVPlant Terrain command definition"

    def GetResources(self):
        return {'Pixmap': str(os.path.join(DirIcons, "terrain.svg")),
                'MenuText': "Terrain",
                'Accel': "S, T",
                'ToolTip': "Creates a Terrain object from setup dialog."}

    def IsActive(self):
        return (not (FreeCAD.ActiveDocument is None) and
                not (FreeCAD.ActiveDocument.getObject("Site") is None) and
                (FreeCAD.ActiveDocument.getObject("Terrain") is None))

    def Activated(self):
        makeTerrain()
        # task = _TerrainTaskPanel()
        # FreeCADGui.Control.showDialog(task)
        return


if FreeCAD.GuiUp:
    FreeCADGui.addCommand('Terrain', _CommandTerrain())'''