# marche_a_lombre.marche_a_lombre_algorithm

View Source

      1# -*- coding: utf-8 -*-
      2
      3# /***************************************************************************
      4#  MarcheALOmbre
      5#                                  A QGIS plugin
      6#  This plugin calculates for a given hike the shady and sunny parts
      7#  Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
      8#                               -------------------
      9#         begin                : 2025-12-11
     10#         copyright            : (C) 2025 by Yolanda Seifert
     11#         email                : yolanda.seifert@gmx.de
     12#  ***************************************************************************/
     13
     14# /***************************************************************************
     15#  *                                                                         *
     16#  *   This program is free software; you can redistribute it and/or modify  *
     17#  *   it under the terms of the GNU General Public License as published by  *
     18#  *   the Free Software Foundation; either version 2 of the License, or     *
     19#  *   (at your option) any later version.                                   *
     20#  *                                                                         *
     21#  ***************************************************************************/
     22
     23__author__ = 'Yolanda Seifert'
     24__date__ = '2025-12-11'
     25__copyright__ = '(C) 2025 by Yolanda Seifert'
     26
     27# This will get replaced with a git SHA1 when you do a git archive
     28
     29__revision__ = '$Format:%H$'
     30
     31import os
     32import inspect
     33import math
     34import processing
     35from osgeo import gdal
     36from qgis.PyQt.QtGui import QIcon, QColor
     37from qgis.PyQt.QtCore import QCoreApplication, QVariant
     38from qgis.core import (QgsProcessing,
     39                        QgsFeatureSink,
     40                        QgsProcessingAlgorithm,
     41                        QgsProcessingParameterFeatureSource,
     42                        QgsProcessingParameterFeatureSink,
     43                        QgsProcessingParameterDateTime,
     44                        QgsProcessingParameterNumber,
     45                        QgsProcessingParameterPoint,
     46                        QgsProcessingParameterBoolean,
     47                        QgsPoint,
     48                        QgsProcessingException,
     49                        QgsProcessingParameterRasterDestination,
     50                        QgsProcessingParameterFileDestination,
     51                        QgsCoordinateReferenceSystem,
     52                        QgsFields,
     53                        QgsWkbTypes,
     54                        QgsFeature,
     55                        QgsGeometry,
     56                        QgsField,
     57                        QgsProcessingUtils,
     58                        QgsCategorizedSymbolRenderer,
     59                        QgsRendererCategory,
     60                        QgsSymbol,
     61                        QgsSimpleMarkerSymbolLayer,
     62                        QgsVectorLayer,
     63                        QgsProperty,
     64                        QgsSymbolLayer)
     65
     66
     67from .mns_downloader import MNSDownloader
     68from .trail import Trail
     69from .shadow_calculator import ShadowCalculator
     70
     71
     72class MarcheALOmbreAlgorithm(QgsProcessingAlgorithm):
     73    """
     74    This is an example algorithm that takes a vector layer and
     75    creates a new identical one.
     76
     77    It is meant to be used as an example of how to create your own
     78    algorithms and explain methods and variables used to do it. An
     79    algorithm like this will be available in all elements, and there
     80    is not need for additional work.
     81
     82    All Processing algorithms should extend the QgsProcessingAlgorithm
     83    class.
     84    """
     85
     86    # Constants used to refer to parameters and outputs. They will be
     87    # used when calling the algorithm from another algorithm, or when
     88    # calling from the QGIS console.
     89
     90    OUTPUT = 'OUTPUT'
     91    INPUT = 'INPUT'
     92    DEPARTURE_TIME = 'DEPARTURE_TIME'
     93    HIKING_SPEED = 'HIKING_SPEED'
     94    PICNIC_POINT = 'PICNIC_POINT'
     95    PICNIC_DURATION = 'PICNIC_DURATION'
     96    REVERSE_DIRECTION = 'REVERSE_DIRECTION'
     97    BUFFER_MODE = 'BUFFER_MODE'
     98    OUTPUT_POINTS = 'OUTPUT_POINTS'
     99    LOW_RES_MNS = 'LOW_RES_MNS'
    100    OUTPUT_CSV = 'OUTPUT_CSV'
    101
    102    def initAlgorithm(self, config):
    103        """
    104        Here we define the inputs and output of the algorithm, along
    105        with some other properties.
    106        """
    107
    108        # We add the input vector features source. It can have any kind of
    109        # geometry.
    110        self.addParameter(
    111            QgsProcessingParameterFeatureSource(
    112                self.INPUT,
    113                self.tr('Input layer - tracks'),
    114                [QgsProcessing.TypeVectorAnyGeometry]
    115            )
    116        )
    117
    118        self.addParameter(
    119            QgsProcessingParameterDateTime(
    120                self.DEPARTURE_TIME,
    121                self.tr('Departure Date and Time (Local Time)'),
    122                type=QgsProcessingParameterDateTime.DateTime  # Allows selecting both Date and Time
    123            )
    124        )
    125
    126        self.addParameter(
    127            QgsProcessingParameterNumber(
    128                self.HIKING_SPEED,
    129                self.tr('Average Hiking Speed (km/h)'),
    130                type=QgsProcessingParameterNumber.Double,
    131                defaultValue=5.0,  # A standard hiking speed
    132                minValue=0.0
    133            )
    134        )
    135
    136        self.addParameter(
    137            QgsProcessingParameterPoint(
    138                self.PICNIC_POINT,
    139                self.tr('Picnic Break Location (1h stop)'),
    140                optional=True  # Optional: user might not want a break
    141            )
    142        )
    143
    144        self.addParameter(
    145            QgsProcessingParameterNumber(
    146                self.PICNIC_DURATION,
    147                self.tr('Picnic Duration (minutes)'),
    148                type=QgsProcessingParameterNumber.Integer,
    149                defaultValue=60,
    150                minValue=0,
    151                optional=True
    152            )
    153        )
    154
    155        self.addParameter(
    156            QgsProcessingParameterBoolean(
    157                self.REVERSE_DIRECTION,
    158                self.tr('Reverse Trail Direction (Finish to Start)'),
    159                defaultValue=False
    160            )
    161        )
    162
    163        self.addParameter(
    164            QgsProcessingParameterBoolean(
    165                self.BUFFER_MODE,
    166                self.tr('Calculate with Buffer (Center, Left 5m, Right 5m)'),
    167                defaultValue=False
    168            )
    169        )
    170
    171
    172        # We add a feature sink in which to store our processed features (this
    173        # usually takes the form of a newly created vector layer when the
    174        # algorithm is run in QGIS).
    175        # self.addParameter(
    176        #     QgsProcessingParameterFeatureSink(
    177        #         self.OUTPUT,
    178        #         self.tr('Output layer')
    179        #     )
    180        # )
    181        self.addParameter(
    182            QgsProcessingParameterRasterDestination(
    183                self.OUTPUT,
    184                self.tr('MNS')
    185            )
    186        )
    187
    188        self.addParameter(
    189            QgsProcessingParameterFeatureSink(
    190                self.OUTPUT_POINTS,
    191                self.tr('Densified Trail Points'),
    192                type=QgsProcessing.TypeVectorPoint
    193            )
    194        )
    195
    196        self.addParameter(
    197            QgsProcessingParameterRasterDestination(
    198                self.LOW_RES_MNS,
    199                self.tr('Low Resolution MNS (Horizon)')
    200            )
    201        )
    202
    203        self.addParameter(
    204            QgsProcessingParameterFileDestination(
    205                self.OUTPUT_CSV,
    206                self.tr('Statistics Table'),
    207                fileFilter='CSV files (*.csv)'
    208            )
    209        )
    210
    211    def processAlgorithm(self, parameters, context, feedback):
    212        """
    213        Here is where the processing itself takes place.
    214        """
    215
    216        # Retrieve the feature source and sink. The 'dest_id' variable is used
    217        # to uniquely identify the feature sink, and must be included in the
    218        # dictionary returned by the processAlgorithm function.
    219        source = self.parameterAsSource(parameters, self.INPUT, context)
    220        departure_dt = self.parameterAsDateTime(parameters, self.DEPARTURE_TIME, context)
    221        speed = self.parameterAsDouble(parameters, self.HIKING_SPEED, context)
    222        picnic_point = self.parameterAsPoint(parameters, self.PICNIC_POINT, context)
    223        picnic_duration = self.parameterAsDouble(parameters, self.PICNIC_DURATION, context)
    224        picnic_point_crs = self.parameterAsPointCrs(parameters, self.PICNIC_POINT, context)
    225        reverse_direction = self.parameterAsBool(parameters, self.REVERSE_DIRECTION, context)
    226        buffer_mode = self.parameterAsBool(parameters, self.BUFFER_MODE, context)
    227        csv_path = self.parameterAsFileOutput(parameters, self.OUTPUT_CSV, context)
    228
    229        if not picnic_point or (picnic_point.x() == 0.0 and picnic_point.y() == 0.0):
    230            picnic_point_crs = None
    231
    232        # Compute the number of steps to display within the progress bar and
    233        # get features from source
    234        total = 100.0 / source.featureCount() if source.featureCount() else 0
    235        features = source.getFeatures()
    236
    237        ########################## TRAIL PROCESSING #########################
    238        feedback.pushInfo("Processing trail...")
    239        source_crs = source.sourceCrs()
    240        if not source_crs.isValid():
    241            raise QgsProcessingException(
    242                "Input layer has no valid CRS. Please define the layer CRS before running this algorithm."
    243            )
    244
    245        trail = Trail(
    246            max_sep=10,
    247            speed=speed,
    248            source_crs=source_crs,
    249            transform_context=context.transformContext(),
    250            feedback=feedback
    251        )
    252        trail.process_trail(source_tracks=source,
    253                            start_time=departure_dt,
    254                            break_point=picnic_point,
    255                            picnic_duration=picnic_duration,
    256                            reverse=reverse_direction,
    257                            buffer=buffer_mode,
    258                            project_crs=picnic_point_crs)
    259
    260        ########################## MNT DOWNLOAD (For Z Values) ##################
    261        # Generate a temporary file path for the MNT
    262        mnt_path = QgsProcessingUtils.generateTempFilename('mnt_elevation.tif')
    263        target_crs = trail.target_crs
    264        downloader = MNSDownloader(
    265            crs=target_crs,
    266            transform_context=context.transformContext(),
    267            feedback=feedback)
    268        target_resolution = 0.5
    269        # Download MNT (mns=False)
    270        success_mnt = downloader.read_tif(
    271            extent=trail.extent,
    272            resolution=target_resolution*2, # less variation in trail elevation so lower resolution necessary
    273            output_path=mnt_path,
    274            input_crs=target_crs,
    275            is_mns=False
    276        )
    277
    278        if not success_mnt:
    279            raise QgsProcessingException("Failed to download MNT data.")
    280        feedback.pushInfo(f"Downloading MNT for trail elevation data...")
    281
    282        # Integrate MNT into Trail (Sample Z values)
    283        feedback.pushInfo("Sampling elevation for trail points...")
    284        trail.sample_elevation(mnt_path)
    285
    286        ########################## MNS DOWNLOAD (for Shade) ############################
    287        output_path = self.parameterAsOutputLayer(parameters, self.OUTPUT, context)
    288        low_res_path = self.parameterAsOutputLayer(parameters, self.LOW_RES_MNS, context)
    289
    290        downloader = MNSDownloader(
    291            crs=target_crs,
    292            transform_context=context.transformContext(),
    293            feedback=feedback)
    294
    295        downloader.download_dual_quality_mns(
    296            trail_extent=trail.extent,
    297            high_res_path=output_path,
    298            low_res_path=low_res_path,
    299            trail_lat=trail.center_lat,
    300            input_crs=target_crs,
    301            high_res=target_resolution, # meters per pixel
    302        )
    303
    304        # Open MNS raster with GDAL
    305        ds = gdal.Open(output_path)
    306        if ds is None:
    307            raise QgsProcessingException("Could not open downloaded MNS.")
    308
    309        mns_band = ds.GetRasterBand(1)
    310        mns_array = mns_band.ReadAsArray() # Returns numpy array
    311        geo_transform = ds.GetGeoTransform()
    312
    313        ds = None # Close dataset
    314
    315        ########################## CALCULATE SHADOWS ############################
    316        feedback.pushInfo("Calculating shadows...")
    317
    318        calculator = ShadowCalculator(
    319            high_res_path=output_path,
    320            low_res_path=low_res_path
    321        )
    322
    323        shadow_results = calculator.calculate_shadows(
    324            trail_points=trail.trail_points,
    325            max_dist_m=20000
    326        )
    327
    328        ########################## CALCULATE STATISTICS ############################
    329        total_points = len(shadow_results)
    330        shady_points = sum(shadow_results)
    331        sunny_points = total_points - shady_points
    332
    333        percent_shady = (shady_points / total_points * 100) if total_points > 0 else 0.0
    334        percent_sunny = (sunny_points / total_points * 100) if total_points > 0 else 0.0
    335
    336        feedback.pushInfo(f"Statistics: {percent_shady:.1f}% Shady, {percent_sunny:.1f}% Sunny")
    337
    338        if csv_path:
    339            try:
    340                with open(csv_path, 'w') as f:
    341                    # Header
    342                    f.write("Duration (min), % Shady, % Sunny, Time spent in Shadow (min), Time spent in Sun (min), Shady Points, Sunny Points, Total Points\n")
    343
    344                    duration_min = 0
    345                    if trail.trail_points:
    346                        start = trail.trail_points[0].datetime
    347                        end = trail.trail_points[-1].datetime
    348                        duration_min = start.secsTo(end) / 60
    349                    time_shady = int(duration_min*percent_shady/100)
    350                    time_sunny = int(duration_min*percent_sunny/100)
    351
    352                    f.write(f"{duration_min:.1f},{percent_shady:.2f},{percent_sunny:.2f},{time_shady},{time_sunny},{shady_points},{sunny_points},{total_points}\n")
    353            except Exception as e:
    354                feedback.reportError(f"Failed to write CSV: {e}")
    355
    356        ########################## WRITE OUTPUT POINTS ##########################
    357        # Define attribute table columns
    358        fields = QgsFields()
    359        fields.append(QgsField("id", QVariant.Int))
    360        fields.append(QgsField("status", QVariant.String))     # Sunny / Shady
    361        fields.append(QgsField("is_shadow", QVariant.Int))     # 0 / 1
    362        fields.append(QgsField("x_proj", QVariant.Double))      # Lambert-93 X
    363        fields.append(QgsField("y_proj", QVariant.Double))      # Lambert-93 Y
    364        fields.append(QgsField("z_mnt", QVariant.Double))      # Altitude
    365        fields.append(QgsField("latitude", QVariant.Double))   # WGS84 Lat
    366        fields.append(QgsField("longitude", QVariant.Double))  # WGS84 Lon
    367        fields.append(QgsField("arrival_time", QVariant.DateTime)) # Time
    368        fields.append(QgsField("elevation_deg", QVariant.Double)) # Sun elevation angle
    369        fields.append(QgsField("azimuth_deg", QVariant.Double)) # Sun azimtuh angle
    370        fields.append(QgsField("course", QVariant.Double))
    371
    372        (point_sink, point_dest_id) = self.parameterAsSink(
    373            parameters,
    374            self.OUTPUT_POINTS,
    375            context,
    376            fields,
    377            QgsWkbTypes.PointZ,
    378            QgsCoordinateReferenceSystem(target_crs)
    379        )
    380        self.point_dest_id = point_dest_id
    381
    382        if point_sink is None:
    383            raise QgsProcessingException("Could not create point sink")
    384
    385        feedback.pushInfo("Writing results...")
    386        for i, tp in enumerate(trail.trail_points):
    387            if feedback.isCanceled(): break
    388
    389            # Calculate direction to next point for visualization
    390            course = 0.0
    391            if i < len(trail.trail_points) - 1:
    392                next_tp = trail.trail_points[i+1]
    393                course = QgsPoint(tp.x, tp.y).azimuth(QgsPoint(next_tp.x, next_tp.y))
    394            elif i > 0:
    395                course = prev_course
    396            prev_course = course
    397
    398            feat = QgsFeature(fields)
    399            geom = QgsGeometry.fromPoint(QgsPoint(tp.x, tp.y, tp.z))
    400            feat.setGeometry(geom)
    401
    402            is_shadow = shadow_results[i]
    403            status_str = "Shadow" if is_shadow == 1 else "Sun"
    404
    405            feat.setAttributes([
    406                i,
    407                status_str,
    408                int(is_shadow),
    409                tp.x,
    410                tp.y,
    411                tp.z,
    412                tp.lat,
    413                tp.lon,
    414                tp.datetime,
    415                math.degrees(tp.solar_pos[0]),
    416                math.degrees(tp.solar_pos[1]),
    417                course
    418            ])
    419            point_sink.addFeature(feat, QgsFeatureSink.FastInsert)
    420            feedback.setProgress(int(i * total))
    421
    422        # Return the results of the algorithm. In this case our only result is
    423        # the feature sink which contains the processed features, but some
    424        # algorithms may return multiple feature sinks, calculated numeric
    425        # statistics, etc. These should all be included in the returned
    426        # dictionary, with keys matching the feature corresponding parameter
    427        # or output names.
    428        self.results = {
    429            self.OUTPUT: output_path,
    430            self.OUTPUT_POINTS: point_dest_id,
    431            self.LOW_RES_MNS: low_res_path,
    432            self.OUTPUT_CSV: csv_path,
    433            'detected_crs': target_crs
    434        }
    435        return self.results
    436
    437    def postProcessAlgorithm(self, context, feedback):
    438        # Set project CRS to match the detected region
    439        detected_crs_str = self.results.get('detected_crs')
    440        if detected_crs_str:
    441            detected_crs = QgsCoordinateReferenceSystem(detected_crs_str)
    442            project = context.project()
    443
    444            if detected_crs.isValid() and project.crs().authid() != detected_crs.authid():
    445                feedback.pushInfo(f"Changing project CRS from {project.crs().authid()} to {detected_crs.authid()} to match detected region")
    446                project.setCrs(detected_crs)
    447        # style the output trail depending on Sun/Shadow
    448        layer = QgsProcessingUtils.mapLayerFromString(self.point_dest_id, context)
    449
    450        if layer:
    451            categories = []
    452            styles = [
    453                (0, "Sun", "gold"),
    454                (1, "Shadow", "#1f78b4")
    455            ]
    456
    457            for val, label, color in styles:
    458                sym_layer = QgsSimpleMarkerSymbolLayer.create({
    459                    'name': 'filled_arrowhead',
    460                    'color': color,
    461                    'size': '2.5',
    462                })
    463
    464                # Direction of arrow based on course field
    465                prop_angle = QgsProperty.fromExpression('"course" - 90')
    466                sym_layer.setDataDefinedProperty(QgsSymbolLayer.PropertyAngle, prop_angle)
    467
    468                sym = QgsSymbol.defaultSymbol(layer.geometryType())
    469                sym.changeSymbolLayer(0, sym_layer)
    470                categories.append(QgsRendererCategory(val, sym, label, True))
    471
    472            renderer = QgsCategorizedSymbolRenderer("is_shadow", categories)
    473            layer.setRenderer(renderer)
    474            layer.triggerRepaint()
    475
    476        csv_path = self.results.get(self.OUTPUT_CSV)
    477        if csv_path and os.path.exists(csv_path):
    478            # Create delimited text layer from csv
    479            path = csv_path.replace('\\', '/')
    480            uri = f"file:///{path}?type=csv&watchFile=no"
    481
    482            vlayer = QgsVectorLayer(uri, "Shadow Statistics", "delimitedtext")
    483            if vlayer.isValid():
    484                context.project().addMapLayer(vlayer)
    485        return self.results
    486
    487    def helpUrl(self):
    488        """
    489        Returns the URL to the help page
    490        """
    491        cmd_folder = os.path.split(inspect.getfile(inspect.currentframe()))[0]
    492        help_path = os.path.join(cmd_folder, 'help', 'build', 'html', 'index.html')
    493
    494        if os.path.exists(help_path):
    495            return "file://" + help_path
    496
    497        # Fallback if local help isn't built yet
    498        return "https://github.com/yolanda225/qgis-marche-a-lombre"
    499
    500    def shortHelpString(self):
    501        """
    502        Returns a brief description of the algorithm
    503        (appears in the right-hand panel of the dialog)
    504        """
    505        return self.tr("This algorithm calculates the shady and sunny portions of a trail "
    506                       "based on LiDAR HD data (MNS) and the sun's position at the time "
    507                       "you are physically at each point on the hike.")
    508
    509    def name(self):
    510        """
    511        Returns the algorithm name, used for identifying the algorithm. This
    512        string should be fixed for the algorithm, and must not be localised.
    513        The name should be unique within each provider. Names should contain
    514        lowercase alphanumeric characters only and no spaces or other
    515        formatting characters.
    516        """
    517        return "marche_a_lombre"
    518
    519    def displayName(self):
    520        """
    521        Returns the translated algorithm name, which should be used for any
    522        user-visible display of the algorithm name.
    523        """
    524        return self.tr(self.name())
    525
    526    def group(self):
    527        """
    528        Returns the name of the group this algorithm belongs to. This string
    529        should be localised.
    530        """
    531        return self.tr(self.groupId())
    532
    533    def groupId(self):
    534        """
    535        Returns the unique ID of the group this algorithm belongs to. This
    536        string should be fixed for the algorithm, and must not be localised.
    537        The group id should be unique within each provider. Group id should
    538        contain lowercase alphanumeric characters only and no spaces or other
    539        formatting characters.
    540        """
    541        return ''
    542
    543    def tr(self, string):
    544        return QCoreApplication.translate('Processing', string)
    545
    546    def icon(self):
    547        cmd_folder = os.path.split(inspect.getfile(inspect.currentframe()))[0]
    548        icon = QIcon(os.path.join(os.path.join(cmd_folder, 'logo.png')))
    549        return icon
    550
    551    def createInstance(self):
    552        return MarcheALOmbreAlgorithm()
    553
    554