From b12c31c8c4937cceb15da0a5339c76e7a5e9ffe9 Mon Sep 17 00:00:00 2001
From: ptresson <paul.tresson@ird.fr>
Date: Fri, 30 Aug 2024 15:49:20 +0200
Subject: [PATCH] add main code (untested yet)
---
src/__init__.py | 0
src/iamap/__init__.py | 10 +
src/iamap/clustering.py | 494 +++++++++++
.../packages_installer_dialog.py | 382 ++++++++
.../packages_installer_dialog.ui | 65 ++
src/iamap/dialogs/resizable_message_box.py | 20 +
src/iamap/encoder.py | 817 ++++++++++++++++++
src/iamap/iamap.py | 234 +++++
src/iamap/icons/__init__.py | 16 +
.../icons/chart-scatter-3d-svgrepo-com.svg | 4 +
.../chart-scatterplot-solid-svgrepo-com.svg | 19 +
src/iamap/icons/encoder_tool.svg | 92 ++
src/iamap/icons/forest-svgrepo-com.svg | 26 +
.../icons/network-base-solid-svgrepo-com.svg | 22 +
.../icons/network-platform-svgrepo-com.svg | 15 +
...atter-plot-outline-alerted-svgrepo-com.svg | 6 +
src/iamap/icons/scatter-plot-svgrepo-com.svg | 6 +
src/iamap/provider.py | 46 +
src/iamap/random_forest.py | 599 +++++++++++++
src/iamap/reduction.py | 568 ++++++++++++
src/iamap/requirements.txt | 4 +
src/iamap/similarity.py | 473 ++++++++++
src/iamap/utils/geo.py | 118 +++
src/iamap/utils/torchgeo.py | 44 +
24 files changed, 4080 insertions(+)
create mode 100644 src/__init__.py
create mode 100644 src/iamap/__init__.py
create mode 100644 src/iamap/clustering.py
create mode 100644 src/iamap/dialogs/packages_installer/packages_installer_dialog.py
create mode 100644 src/iamap/dialogs/packages_installer/packages_installer_dialog.ui
create mode 100644 src/iamap/dialogs/resizable_message_box.py
create mode 100644 src/iamap/encoder.py
create mode 100644 src/iamap/iamap.py
create mode 100644 src/iamap/icons/__init__.py
create mode 100644 src/iamap/icons/chart-scatter-3d-svgrepo-com.svg
create mode 100644 src/iamap/icons/chart-scatterplot-solid-svgrepo-com.svg
create mode 100644 src/iamap/icons/encoder_tool.svg
create mode 100644 src/iamap/icons/forest-svgrepo-com.svg
create mode 100644 src/iamap/icons/network-base-solid-svgrepo-com.svg
create mode 100644 src/iamap/icons/network-platform-svgrepo-com.svg
create mode 100644 src/iamap/icons/scatter-plot-outline-alerted-svgrepo-com.svg
create mode 100644 src/iamap/icons/scatter-plot-svgrepo-com.svg
create mode 100644 src/iamap/provider.py
create mode 100644 src/iamap/random_forest.py
create mode 100644 src/iamap/reduction.py
create mode 100644 src/iamap/requirements.txt
create mode 100644 src/iamap/similarity.py
create mode 100644 src/iamap/utils/geo.py
create mode 100644 src/iamap/utils/torchgeo.py
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/iamap/__init__.py b/src/iamap/__init__.py
new file mode 100644
index 0000000..e8e5304
--- /dev/null
+++ b/src/iamap/__init__.py
@@ -0,0 +1,10 @@
+import os
+import inspect
+
+cmd_folder = os.path.split(inspect.getfile(inspect.currentframe()))[0]
+
+def classFactory(iface):
+ from .dialogs.packages_installer import packages_installer_dialog
+ packages_installer_dialog.check_required_packages_and_install_if_necessary(iface=iface)
+ from .iamap import IAMap
+ return IAMap(iface, cmd_folder)
diff --git a/src/iamap/clustering.py b/src/iamap/clustering.py
new file mode 100644
index 0000000..f2b19b0
--- /dev/null
+++ b/src/iamap/clustering.py
@@ -0,0 +1,494 @@
+import os
+import numpy as np
+from pathlib import Path
+from typing import Dict, Any
+import joblib
+
+import rasterio
+from rasterio import windows
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (Qgis,
+ QgsGeometry,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterEnum,
+ QgsCoordinateTransform,
+ QgsProcessingException,
+ QgsProcessingAlgorithm,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterBand,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterCrs,
+ QgsProcessingParameterDefinition,
+ )
+
+
+from sklearn.cluster import KMeans
+
+import json
+
+
+
+class ClusterAlgorithm(QgsProcessingAlgorithm):
+ """
+ """
+
+ INPUT = 'INPUT'
+ BANDS = 'BANDS'
+ EXTENT = 'EXTENT'
+ LOAD = 'LOAD'
+ OUTPUT = 'OUTPUT'
+ RESOLUTION = 'RESOLUTION'
+ CRS = 'CRS'
+ CLUSTERS = 'CLUSTERS'
+ SUBSET = 'SUBSET'
+ METHOD = 'METHOD'
+ SAVE_MODEL = 'SAVE_MODEL'
+
+
+ def initAlgorithm(self, config=None):
+ """
+ Here we define the inputs and output of the algorithm, along
+ with some other properties.
+ """
+ cwd = Path(__file__).parent.absolute()
+
+ self.addParameter(
+ QgsProcessingParameterRasterLayer(
+ name=self.INPUT,
+ description=self.tr(
+ 'Input raster layer or image file path'),
+ defaultValue=os.path.join(cwd,'rasters','test.tif'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBand(
+ name=self.BANDS,
+ description=self.tr(
+ 'Selected Bands (defaults to all bands selected)'),
+ defaultValue=None,
+ parentLayerParameterName=self.INPUT,
+ optional=True,
+ allowMultiple=True,
+ )
+ )
+
+ crs_param = QgsProcessingParameterCrs(
+ name=self.CRS,
+ description=self.tr('Target CRS (default to original CRS)'),
+ optional=True,
+ )
+
+ res_param = QgsProcessingParameterNumber(
+ name=self.RESOLUTION,
+ description=self.tr(
+ 'Target resolution in meters (default to native resolution)'),
+ type=QgsProcessingParameterNumber.Double,
+ optional=True,
+ minValue=0,
+ maxValue=100000
+ )
+
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ name=self.EXTENT,
+ description=self.tr(
+ 'Processing extent (default to the entire image)'),
+ optional=True
+ )
+ )
+ self.method_opt = ['K-means', '--Empty--']
+ self.addParameter (
+ QgsProcessingParameterEnum(
+ name = self.METHOD,
+ description = self.tr(
+ 'Method for the dimension reduction'),
+ defaultValue = 0,
+ options = self.method_opt,
+
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name=self.CLUSTERS,
+ description=self.tr(
+ 'Number of target clusters'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue = 5,
+ minValue=1,
+ maxValue=1024
+ )
+ )
+ subset_param = QgsProcessingParameterNumber(
+ name=self.SUBSET,
+ description=self.tr(
+ 'Select a subset of random pixels of the image to fit transform'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=None,
+ minValue=1,
+ maxValue=10_000,
+ optional=True,
+ )
+
+ save_param = QgsProcessingParameterBoolean(
+ self.SAVE_MODEL,
+ self.tr("Save projection model after fit."),
+ defaultValue=True
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT,
+ self.tr(
+ "Output directory (choose the location that the image features will be saved)"),
+ defaultValue=os.path.join(cwd,'features'),
+ )
+ )
+
+
+ for param in (crs_param, res_param, subset_param, save_param):
+ param.setFlags(
+ param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(param)
+
+
+
+ def processAlgorithm(self, parameters, context, feedback):
+ """
+ Here is where the processing itself takes place.
+ """
+ self.process_options(parameters, context, feedback)
+
+ input_bands = [i_band -1 for i_band in self.selected_bands]
+
+ if self.method == 'K-means':
+ proj = KMeans(int(self.nclusters))
+ save_file = 'kmeans_cluster.pkl'
+ params = proj.get_params()
+
+ out_path = os.path.join(self.output_dir, save_file)
+
+ with rasterio.open(self.rlayer_path) as ds:
+ transform = ds.transform
+ crs = ds.crs
+ win = windows.from_bounds(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ transform=transform
+ )
+ raster = ds.read(window=win)
+ transform = ds.window_transform(win)
+ raster = np.transpose(raster, (1,2,0))
+ raster = raster[:,:,input_bands]
+
+ feedback.pushInfo(f'{raster.shape}')
+ feedback.pushInfo(f'{raster.reshape(-1, raster.shape[0]).shape}')
+
+ if self.subset:
+ feedback.pushInfo(f'Using a random subset of {self.subset} pixels')
+ fit_raster = raster.reshape(-1, raster.shape[-1])
+ nsamples = fit_raster.shape[0]
+
+ # Generate random indices to select subset_size number of samples
+ np.random.seed(42)
+ random_indices = np.random.choice(nsamples, size=self.subset, replace=False)
+ fit_raster = fit_raster[random_indices,:]
+ feedback.pushInfo(f'Starting fit\n')
+ proj.fit(fit_raster)
+ if self.save_model:
+ joblib.dump(proj, out_path)
+
+ feedback.pushInfo(f'starting inference\n')
+ proj_img = proj.predict(raster.reshape(-1, raster.shape[-1]))
+
+
+ else:
+ proj_img = proj.fit_predict(raster.reshape(-1, raster.shape[-1]))
+ if self.save_model:
+ joblib.dump(proj, out_path)
+
+ proj_img = proj_img.reshape((raster.shape[0], raster.shape[1],-1))
+ height, width, channels = proj_img.shape
+
+
+ dst_path = os.path.join(self.output_dir,'cluster.tif')
+ params_file = os.path.join(self.output_dir, 'cluster_parameters.json')
+
+
+ if os.path.exists(dst_path):
+ i = 1
+ while True:
+ modified_output_file = os.path.join(self.output_dir, f"cluster_{i}.tif")
+ if not os.path.exists(modified_output_file):
+ dst_path = modified_output_file
+ break
+ i += 1
+
+ if os.path.exists(params_file):
+ i = 1
+ while True:
+ modified_output_file_params = os.path.join(self.output_dir, f"cluster_parameters_{i}.json")
+ if not os.path.exists(modified_output_file_params):
+ params_file = modified_output_file_params
+ break
+ i += 1
+
+
+ with rasterio.open(dst_path, 'w', driver='GTiff',
+ height=height, width=width, count=channels, dtype='int8',
+ crs=crs, transform=transform) as dst_ds:
+ dst_ds.write(np.transpose(proj_img, (2, 0, 1)))
+
+ with open(params_file, 'w') as f:
+ json.dump(params, f, indent=4)
+ feedback.pushInfo(f"Parameters saved to {params_file}")
+
+ parameters['OUTPUT_RASTER']=dst_path
+
+ return {'OUTPUT_RASTER':dst_path}
+
+ def process_options(self,parameters, context, feedback):
+ self.iPatch = 0
+
+ self.feature_dir = ""
+
+ feedback.pushInfo(
+ f'PARAMETERS :\n{parameters}')
+
+ feedback.pushInfo(
+ f'CONTEXT :\n{context}')
+
+ feedback.pushInfo(
+ f'FEEDBACK :\n{feedback}')
+
+ rlayer = self.parameterAsRasterLayer(
+ parameters, self.INPUT, context)
+
+ if rlayer is None:
+ raise QgsProcessingException(
+ self.invalidRasterError(parameters, self.INPUT))
+
+ self.selected_bands = self.parameterAsInts(
+ parameters, self.BANDS, context)
+
+ if len(self.selected_bands) == 0:
+ self.selected_bands = list(range(1, rlayer.bandCount()+1))
+
+ if max(self.selected_bands) > rlayer.bandCount():
+ raise QgsProcessingException(
+ self.tr("The chosen bands exceed the largest band number!")
+ )
+
+ self.nclusters = self.parameterAsInt(
+ parameters, self.CLUSTERS, context)
+ self.subset = self.parameterAsInt(
+ parameters, self.SUBSET, context)
+ method_idx = self.parameterAsEnum(
+ parameters, self.METHOD, context)
+ self.method = self.method_opt[method_idx]
+
+ res = self.parameterAsDouble(
+ parameters, self.RESOLUTION, context)
+ crs = self.parameterAsCrs(
+ parameters, self.CRS, context)
+ extent = self.parameterAsExtent(
+ parameters, self.EXTENT, context)
+ self.output_dir = self.parameterAsString(
+ parameters, self.OUTPUT, context)
+ self.save_model = self.parameterAsBoolean(
+ parameters, self.SAVE_MODEL, context)
+
+ rlayer_data_provider = rlayer.dataProvider()
+
+ # handle crs
+ if crs is None or not crs.isValid():
+ crs = rlayer.crs()
+ """
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}') # 0 for meters, 6 for degrees, 9 for unknown
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ if crs.mapUnits() == Qgis.DistanceUnit.Degrees:
+ crs = self.estimate_utm_crs(rlayer.extent())
+
+ # target crs should use meters as units
+ if crs.mapUnits() != Qgis.DistanceUnit.Meters:
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}')
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ raise QgsProcessingException(
+ self.tr("Only support CRS with the units as meters")
+ )
+ """
+
+ # 0 for meters, 6 for degrees, 9 for unknown
+ UNIT_METERS = 0
+ UNIT_DEGREES = 6
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ layer_units = 'degrees'
+ else:
+ layer_units = 'meters'
+ # if res is not provided, get res info from rlayer
+ if np.isnan(res) or res == 0:
+ res = rlayer.rasterUnitsPerPixelX() # rasterUnitsPerPixelY() is negative
+ target_units = layer_units
+ else:
+ # when given res in meters by users, convert crs to utm if the original crs unit is degree
+ if crs.mapUnits() != UNIT_METERS: # Qgis.DistanceUnit.Meters:
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ # estimate utm crs based on layer extent
+ crs = self.estimate_utm_crs(rlayer.extent())
+ else:
+ raise QgsProcessingException(
+ f"Resampling of image with the CRS of {crs.authid()} in meters is not supported.")
+ target_units = 'meters'
+ # else:
+ # res = (rlayer_extent.xMaximum() -
+ # rlayer_extent.xMinimum()) / rlayer.width()
+ self.res = res
+
+ # handle extent
+ if extent.isNull():
+ extent = rlayer.extent() # QgsProcessingUtils.combineLayerExtents(layers, crs, context)
+ extent_crs = rlayer.crs()
+ else:
+ if extent.isEmpty():
+ raise QgsProcessingException(
+ self.tr("The extent for processing can not be empty!"))
+ extent_crs = self.parameterAsExtentCrs(
+ parameters, self.EXTENT, context)
+ # if extent crs != target crs, convert it to target crs
+ if extent_crs != crs:
+ transform = QgsCoordinateTransform(
+ extent_crs, crs, context.transformContext())
+ # extent = transform.transformBoundingBox(extent)
+ # to ensure coverage of the transformed extent
+ # convert extent to polygon, transform polygon, then get boundingBox of the new polygon
+ extent_polygon = QgsGeometry.fromRect(extent)
+ extent_polygon.transform(transform)
+ extent = extent_polygon.boundingBox()
+ extent_crs = crs
+
+ # check intersects between extent and rlayer_extent
+ if rlayer.crs() != crs:
+ transform = QgsCoordinateTransform(
+ rlayer.crs(), crs, context.transformContext())
+ rlayer_extent = transform.transformBoundingBox(
+ rlayer.extent())
+ else:
+ rlayer_extent = rlayer.extent()
+ if not rlayer_extent.intersects(extent):
+ raise QgsProcessingException(
+ self.tr("The extent for processing is not intersected with the input image!"))
+
+ img_width_in_extent = round(
+ (extent.xMaximum() - extent.xMinimum())/self.res)
+ img_height_in_extent = round(
+ (extent.yMaximum() - extent.yMinimum())/self.res)
+
+ # Send some information to the user
+ feedback.pushInfo(
+ f'Layer path: {rlayer_data_provider.dataSourceUri()}')
+ # feedback.pushInfo(
+ # f'Layer band scale: {rlayer_data_provider.bandScale(self.selected_bands[0])}')
+ feedback.pushInfo(f'Layer name: {rlayer.name()}')
+ if rlayer.crs().authid():
+ feedback.pushInfo(f'Layer CRS: {rlayer.crs().authid()}')
+ else:
+ feedback.pushInfo(
+ f'Layer CRS in WKT format: {rlayer.crs().toWkt()}')
+ feedback.pushInfo(
+ f'Layer pixel size: {rlayer.rasterUnitsPerPixelX()}, {rlayer.rasterUnitsPerPixelY()} {layer_units}')
+
+ feedback.pushInfo(f'Bands selected: {self.selected_bands}')
+
+ if crs.authid():
+ feedback.pushInfo(f'Target CRS: {crs.authid()}')
+ else:
+ feedback.pushInfo(f'Target CRS in WKT format: {crs.toWkt()}')
+ feedback.pushInfo(f'Target resolution: {self.res} {target_units}')
+ feedback.pushInfo(
+ (f'Processing extent: minx:{extent.xMinimum():.6f}, maxx:{extent.xMaximum():.6f},'
+ f'miny:{extent.yMinimum():.6f}, maxy:{extent.yMaximum():.6f}'))
+ feedback.pushInfo(
+ (f'Processing image size: (width {img_width_in_extent}, '
+ f'height {img_height_in_extent})'))
+
+ self.rlayer_path = rlayer.dataProvider().dataSourceUri()
+
+ feedback.pushInfo(f'Selected bands: {self.selected_bands}')
+
+ ## passing parameters to self once everything has been processed
+ self.extent = extent
+ self.rlayer = rlayer
+ self.crs = crs
+
+
+ # used to handle any thread-sensitive cleanup which is required by the algorithm.
+ def postProcessAlgorithm(self, context, feedback) -> Dict[str, Any]:
+ return {}
+
+
+ def tr(self, string):
+ """
+ Returns a translatable string with the self.tr() function.
+ """
+ return QCoreApplication.translate('Processing', string)
+
+ def createInstance(self):
+ return ClusterAlgorithm()
+
+ def name(self):
+ """
+ Returns the algorithm name, used for identifying the algorithm. This
+ string should be fixed for the algorithm, and must not be localised.
+ The name should be unique within each provider. Names should contain
+ lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return 'cluster'
+
+ def displayName(self):
+ """
+ Returns the translated algorithm name, which should be used for any
+ user-visible display of the algorithm name.
+ """
+ return self.tr('Clustering')
+
+ def group(self):
+ """
+ Returns the name of the group this algorithm belongs to. This string
+ should be localised.
+ """
+ return self.tr('')
+
+ def groupId(self):
+ """
+ Returns the unique ID of the group this algorithm belongs to. This
+ string should be fixed for the algorithm, and must not be localised.
+ The group id should be unique within each provider. Group id should
+ contain lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return ''
+
+ def shortHelpString(self):
+ """
+ Returns a localised short helper string for the algorithm. This string
+ should provide a basic description about what the algorithm does and the
+ parameters and outputs associated with it..
+ """
+ return self.tr("Cluster a raster.")
+
+ def icon(self):
+ return 'E'
+
+
+
diff --git a/src/iamap/dialogs/packages_installer/packages_installer_dialog.py b/src/iamap/dialogs/packages_installer/packages_installer_dialog.py
new file mode 100644
index 0000000..fd128d4
--- /dev/null
+++ b/src/iamap/dialogs/packages_installer/packages_installer_dialog.py
@@ -0,0 +1,382 @@
+"""
+This QGIS plugin requires some Python packages to be installed and available.
+This tool allows to install them in a local directory, if they are not installed yet.
+"""
+
+import importlib
+import logging
+import os
+import subprocess
+import sys
+import traceback
+import urllib
+from dataclasses import dataclass
+from pathlib import Path
+from threading import Thread
+from typing import List
+
+from qgis.PyQt import QtCore, uic
+from qgis.PyQt.QtCore import pyqtSignal
+from qgis.PyQt.QtGui import QCloseEvent
+from qgis.PyQt.QtWidgets import QDialog, QMessageBox, QTextBrowser
+
+# from deepness.common.defines import PLUGIN_NAME
+
+PLUGIN_NAME = "iamap"
+
+PYTHON_VERSION = sys.version_info
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+PLUGIN_ROOT_DIR = os.path.realpath(os.path.abspath(os.path.join(SCRIPT_DIR, '..', '..')))
+PACKAGES_INSTALL_DIR = os.path.join(PLUGIN_ROOT_DIR, f'python{PYTHON_VERSION.major}.{PYTHON_VERSION.minor}')
+
+
+FORM_CLASS, _ = uic.loadUiType(os.path.join(
+ os.path.dirname(__file__), 'packages_installer_dialog.ui'))
+
+_ERROR_COLOR = '#ff0000'
+
+
+@dataclass
+class PackageToInstall:
+ name: str
+ version: str
+ import_name: str # name while importing package
+
+ def __str__(self):
+ return f'{self.name}{self.version}'
+
+
+# REQUIREMENTS_PATH = os.path.join(PLUGIN_ROOT_DIR, 'python_requirements/requirements.txt')
+REQUIREMENTS_PATH = os.path.join(PLUGIN_ROOT_DIR, 'requirements.txt')
+
+with open(REQUIREMENTS_PATH, 'r') as f:
+ raw_txt = f.read()
+
+libraries_versions = {}
+
+for line in raw_txt.split('\n'):
+ if line.startswith('#') or not line.strip():
+ continue
+
+ line = line.split(';')[0]
+
+ if '==' in line:
+ lib, version = line.split('==')
+ libraries_versions[lib] = '==' + version
+ elif '>=' in line:
+ lib, version = line.split('>=')
+ libraries_versions[lib] = '>=' + version
+ elif '<=' in line:
+ lib, version = line.split('<=')
+ libraries_versions[lib] = '<=' + version
+ else:
+ libraries_versions[line] = ''
+
+
+packages_to_install = []
+for lib, version in libraries_versions.items():
+
+ import_name = lib[:-1]
+
+ if lib == 'scikit-learn ':
+ import_name = 'sklearn'
+ if lib == 'umap-learn ':
+ import_name = 'umap'
+
+ packages_to_install.append(
+ PackageToInstall(
+ name=lib,
+ version=version,
+ import_name=import_name
+ )
+ )
+
+
+# packages_to_install = [
+# PackageToInstall(name='scikit-learn', version=libraries_versions['scikit-learn'], import_name='sklearn'),
+# ]
+
+if sys.platform == "linux" or sys.platform == "linux2":
+ # packages_to_install += [
+ # PackageToInstall(name='ime-gpu', version=libraries_versions['onnxruntime-gpu'], import_name='onnxruntime'),
+ # ]
+ PYTHON_EXECUTABLE_PATH = sys.executable
+elif sys.platform == "darwin": # MacOS
+ # packages_to_install += [
+ # PackageToInstall(name='onnxruntime', version=libraries_versions['onnxruntime-gpu'], import_name='onnxruntime'),
+ # ]
+ PYTHON_EXECUTABLE_PATH = str(Path(sys.prefix) / 'bin' / 'python3') # sys.executable yields QGIS in macOS
+elif sys.platform == "win32":
+ # packages_to_install += [
+ # PackageToInstall(name='onnxruntime', version=libraries_versions['onnxruntime-gpu'], import_name='onnxruntime'),
+ # ]
+ PYTHON_EXECUTABLE_PATH = 'python' # sys.executable yields QGis.exe in Windows
+else:
+ raise Exception("Unsupported operating system!")
+
+
+class PackagesInstallerDialog(QDialog, FORM_CLASS):
+ """
+ Dialog witch controls the installation process of packages.
+ UI design defined in the `packages_installer_dialog.ui` file.
+ """
+
+ signal_log_line = pyqtSignal(str) # we need to use signal because we cannot edit GUI from another thread
+
+ INSTALLATION_IN_PROGRESS = False # to make sure we will not start the installation twice
+
+ def __init__(self, iface, parent=None):
+ super(PackagesInstallerDialog, self).__init__(parent)
+ self.setupUi(self)
+ self.iface = iface
+ self.tb = self.textBrowser_log # type: QTextBrowser
+ self._create_connections()
+ self._setup_message()
+ self.aborted = False
+ self.thread = None
+
+ def move_to_top(self):
+ """ Move the window to the top.
+ Although if installed from plugin manager, the plugin manager will move itself to the top anyway.
+ """
+ self.setWindowState((self.windowState() & ~QtCore.Qt.WindowMinimized) | QtCore.Qt.WindowActive)
+
+ if sys.platform == "linux" or sys.platform == "linux2":
+ pass
+ elif sys.platform == "darwin": # MacOS
+ self.raise_() # FIXME: this does not really work, the window is still behind the plugin manager
+ elif sys.platform == "win32":
+ self.activateWindow()
+ else:
+ raise Exception("Unsupported operating system!")
+
+ def _create_connections(self):
+ self.pushButton_close.clicked.connect(self.close)
+ self.pushButton_install_packages.clicked.connect(self._run_packages_installation)
+ self.signal_log_line.connect(self._log_line)
+
+ def _log_line(self, txt):
+ txt = txt \
+ .replace(' ', ' ') \
+ .replace('\n', '<br>')
+ self.tb.append(txt)
+
+ def log(self, txt):
+ self.signal_log_line.emit(txt)
+
+ def _setup_message(self) -> None:
+
+ self.log(f'<h2><span style="color: #000080;"><strong> '
+ f'Plugin {PLUGIN_NAME} - Packages installer </strong></span></h2> \n'
+ f'\n'
+ f'<b>This plugin requires the following Python packages to be installed:</b>')
+
+ for package in packages_to_install:
+ self.log(f'\t- {package.name}{package.version}')
+
+ self.log('\n\n'
+ f'If this packages are not installed in the global environment '
+ f'(or environment in which QGIS is started) '
+ f'you can install these packages in the local directory (which is included to the Python path).\n\n'
+ f'This Dialog does it for you! (Though you can still install these packages manually instead).\n'
+ f'<b>Please click "Install packages" button below to install them automatically, </b>'
+ f'or "Test and Close" if you installed them manually...\n')
+
+ def _run_packages_installation(self):
+ if self.INSTALLATION_IN_PROGRESS:
+ self.log(f'Error! Installation already in progress, cannot start again!')
+ return
+ self.aborted = False
+ self.INSTALLATION_IN_PROGRESS = True
+ self.thread = Thread(target=self._install_packages)
+ self.thread.start()
+
+ def _install_packages(self) -> None:
+ self.log('\n\n')
+ self.log('=' * 60)
+ self.log(f'<h3><b>Attempting to install required packages...</b></h3>')
+ os.makedirs(PACKAGES_INSTALL_DIR, exist_ok=True)
+
+ self._install_pip_if_necessary()
+
+ self.log(f'<h3><b>Attempting to install required packages...</b></h3>\n')
+ try:
+ self._pip_install_packages(packages_to_install)
+ except Exception as e:
+ msg = (f'\n <span style="color: {_ERROR_COLOR};"><b> '
+ f'Packages installation failed with exception: {e}!\n'
+ f'Please try to install the packages again. </b></span>'
+ f'\nCheck if there is no error related to system packages, '
+ f'which may be required to be installed by your system package manager, e.g. "apt". '
+ f'Copy errors from the stack above and google for possible solutions. '
+ f'Please report these as an issue on the plugin repository tracker!')
+ self.log(msg)
+
+ # finally, validate the installation, if there was no error so far...
+ self.log('\n\n <b>Installation of required packages finished. Validating installation...</b>')
+ self._check_packages_installation_and_log()
+ self.INSTALLATION_IN_PROGRESS = False
+
+ def reject(self) -> None:
+ self.close()
+
+ def closeEvent(self, event: QCloseEvent):
+ self.aborted = True
+ if self._check_packages_installation_and_log():
+ event.accept()
+ return
+
+ res = QMessageBox.question(self.iface.mainWindow(),
+ f'{PLUGIN_NAME} - skip installation?',
+ 'Are you sure you want to abort the installation of the required python packages? '
+ 'The plugin may not function correctly without them!',
+ QMessageBox.No, QMessageBox.Yes)
+ log_msg = 'User requested to close the dialog, but the packages are not installed correctly!\n'
+ if res == QMessageBox.Yes:
+ log_msg += 'And the user confirmed to close the dialog, knowing the risk!'
+ event.accept()
+ else:
+ log_msg += 'The user reconsidered their decision, and will try to install the packages again!'
+ event.ignore()
+ log_msg += '\n'
+ self.log(log_msg)
+
+ def _install_pip_if_necessary(self):
+ """
+ Install pip if not present.
+ It happens e.g. in flatpak applications.
+
+ TODO - investigate whether we can also install pip in local directory
+ """
+
+ self.log(f'<h4><b>Making sure pip is installed...</b></h4>')
+ if check_pip_installed():
+ self.log(f'<em>Pip is installed, skipping installation...</em>\n')
+ return
+
+ install_pip_command = [PYTHON_EXECUTABLE_PATH, '-m', 'ensurepip']
+ self.log(f'<em>Running command to install pip: \n $ {" ".join(install_pip_command)} </em>')
+ with subprocess.Popen(install_pip_command,
+ stdout=subprocess.PIPE,
+ universal_newlines=True,
+ stderr=subprocess.STDOUT,
+ env={'SETUPTOOLS_USE_DISTUTILS': 'stdlib'}) as process:
+ try:
+ self._do_process_output_logging(process)
+ except InterruptedError as e:
+ self.log(str(e))
+ return False
+
+ if process.returncode != 0:
+ msg = (f'<span style="color: {_ERROR_COLOR};"><b>'
+ f'pip installation failed! Consider installing it manually.'
+ f'<b></span>')
+ self.log(msg)
+ self.log('\n')
+
+ def _pip_install_packages(self, packages: List[PackageToInstall]) -> None:
+ cmd = [PYTHON_EXECUTABLE_PATH, '-m', 'pip', 'install', '-U', f'--target={PACKAGES_INSTALL_DIR}']
+ cmd_string = ' '.join(cmd)
+
+ for pck in packages:
+ cmd.append(f"{pck}")
+ cmd_string += f"{pck}"
+
+ self.log(f'<em>Running command: \n $ {cmd_string} </em>')
+ with subprocess.Popen(cmd,
+ stdout=subprocess.PIPE,
+ universal_newlines=True,
+ stderr=subprocess.STDOUT) as process:
+ self._do_process_output_logging(process)
+
+ if process.returncode != 0:
+ raise RuntimeError('Installation with pip failed')
+
+ msg = (f'\n<b>'
+ f'Packages installed correctly!'
+ f'<b>\n\n')
+ self.log(msg)
+
+ def _do_process_output_logging(self, process: subprocess.Popen) -> None:
+ """
+ :param process: instance of 'subprocess.Popen'
+ """
+ for stdout_line in iter(process.stdout.readline, ""):
+ if stdout_line.isspace():
+ continue
+ txt = f'<span style="color: #999999;">{stdout_line.rstrip(os.linesep)}</span>'
+ self.log(txt)
+ if self.aborted:
+ raise InterruptedError('Installation aborted by user')
+
+ def _check_packages_installation_and_log(self) -> bool:
+ packages_ok = are_packages_importable()
+ self.pushButton_install_packages.setEnabled(not packages_ok)
+
+ if packages_ok:
+ msg1 = f'All required packages are importable! You can close this window now!'
+ self.log(msg1)
+ return True
+
+ try:
+ import_packages()
+ raise Exception("Unexpected successful import of packages?!? It failed a moment ago, we shouldn't be here!")
+ except Exception:
+ msg_base = '<b>Python packages required by the plugin could not be loaded due to the following error:</b>'
+ logging.exception(msg_base)
+ tb = traceback.format_exc()
+ msg1 = (f'<span style="color: {_ERROR_COLOR};">'
+ f'{msg_base} \n '
+ f'{tb}\n\n'
+ f'<b>Please try installing the packages again.<b>'
+ f'</span>')
+ self.log(msg1)
+
+ return False
+
+
+dialog = None
+
+
+def import_package(package: PackageToInstall):
+ # print(package.import_name)
+ importlib.import_module(package.import_name)
+
+
+def import_packages():
+ for package in packages_to_install:
+ import_package(package)
+
+
+def are_packages_importable() -> bool:
+ try:
+ import_packages()
+ except Exception:
+ logging.exception(f'Python packages required by the plugin could not be loaded due to the following error:')
+ return False
+
+ return True
+
+
+def check_pip_installed() -> bool:
+ try:
+ subprocess.check_output([PYTHON_EXECUTABLE_PATH, '-m', 'pip', '--version'])
+ return True
+ except subprocess.CalledProcessError:
+ return False
+
+
+def check_required_packages_and_install_if_necessary(iface):
+ os.makedirs(PACKAGES_INSTALL_DIR, exist_ok=True)
+ if PACKAGES_INSTALL_DIR not in sys.path:
+ sys.path.append(PACKAGES_INSTALL_DIR) # TODO: check for a less intrusive way to do this
+
+ if are_packages_importable():
+ # if packages are importable we are fine, nothing more to do then
+ return
+
+ global dialog
+ dialog = PackagesInstallerDialog(iface)
+ dialog.setWindowModality(QtCore.Qt.WindowModal)
+ dialog.show()
+ dialog.move_to_top()
diff --git a/src/iamap/dialogs/packages_installer/packages_installer_dialog.ui b/src/iamap/dialogs/packages_installer/packages_installer_dialog.ui
new file mode 100644
index 0000000..de4c975
--- /dev/null
+++ b/src/iamap/dialogs/packages_installer/packages_installer_dialog.ui
@@ -0,0 +1,65 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>PackagesInstallerDialog</class>
+ <widget class="QDialog" name="PackagesInstallerDialog">
+ <property name="geometry">
+ <rect>
+ <x>0</x>
+ <y>0</y>
+ <width>693</width>
+ <height>494</height>
+ </rect>
+ </property>
+ <property name="windowTitle">
+ <string>Featmap - Packages Installer Dialog</string>
+ </property>
+ <layout class="QGridLayout" name="gridLayout_2">
+ <item row="0" column="0">
+ <layout class="QGridLayout" name="gridLayout">
+ <item row="0" column="0">
+ <widget class="QTextBrowser" name="textBrowser_log"/>
+ </item>
+ </layout>
+ </item>
+ <item row="1" column="0">
+ <layout class="QHBoxLayout" name="horizontalLayout">
+ <item>
+ <spacer name="horizontalSpacer">
+ <property name="orientation">
+ <enum>Qt::Horizontal</enum>
+ </property>
+ <property name="sizeHint" stdset="0">
+ <size>
+ <width>40</width>
+ <height>20</height>
+ </size>
+ </property>
+ </spacer>
+ </item>
+ <item>
+ <widget class="QPushButton" name="pushButton_install_packages">
+ <property name="font">
+ <font>
+ <weight>75</weight>
+ <bold>true</bold>
+ </font>
+ </property>
+ <property name="text">
+ <string>Install packages</string>
+ </property>
+ </widget>
+ </item>
+ <item>
+ <widget class="QPushButton" name="pushButton_close">
+ <property name="text">
+ <string>Test and Close</string>
+ </property>
+ </widget>
+ </item>
+ </layout>
+ </item>
+ </layout>
+ </widget>
+ <resources/>
+ <connections/>
+</ui>
diff --git a/src/iamap/dialogs/resizable_message_box.py b/src/iamap/dialogs/resizable_message_box.py
new file mode 100644
index 0000000..b2a6794
--- /dev/null
+++ b/src/iamap/dialogs/resizable_message_box.py
@@ -0,0 +1,20 @@
+from qgis.PyQt.QtWidgets import QMessageBox, QTextEdit
+
+
+class ResizableMessageBox(QMessageBox):
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+ self.setSizeGripEnabled(True)
+
+ def event(self, event):
+ if event.type() in (event.LayoutRequest, event.Resize):
+ if event.type() == event.Resize:
+ res = super().event(event)
+ else:
+ res = False
+ details = self.findChild(QTextEdit)
+ if details:
+ details.setMaximumSize(16777215, 16777215)
+ self.setMaximumSize(16777215, 16777215)
+ return res
+ return super().event(event)
diff --git a/src/iamap/encoder.py b/src/iamap/encoder.py
new file mode 100644
index 0000000..7f4448a
--- /dev/null
+++ b/src/iamap/encoder.py
@@ -0,0 +1,817 @@
+import os
+import time
+import re
+import hashlib
+import numpy as np
+from pathlib import Path
+from typing import Dict, Any, List
+from processing import defaultOutputFolder
+
+import rasterio
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (QgsProcessing, Qgis,
+ QgsGeometry,
+ QgsRasterLayer,
+ QgsRectangle,
+ QgsCoordinateReferenceSystem,
+ QgsUnitTypes,
+ QgsRasterBandStats,
+ QgsCoordinateTransform,
+ QgsFeatureSink,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProcessingAlgorithm,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterBand,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterString,
+ QgsProcessingParameterEnum,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterCrs,
+ QgsProcessingParameterScale,
+ QgsProcessingParameterExpression,
+ QgsProcessingParameterRange,
+ QgsProcessingParameterFeatureSource,
+ QgsProcessingParameterDefinition,
+ QgsProcessingParameterFeatureSink,
+ )
+
+import torch
+import torch.nn as nn
+from torch import Tensor
+import torch.quantization
+from torch.utils.data import DataLoader
+import torchvision.transforms as T
+import kornia.augmentation as K
+import timm
+
+from torchgeo.datasets import RasterDataset, BoundingBox,stack_samples
+from torchgeo.samplers import GridGeoSampler, Units
+from torchgeo.transforms import AugmentationSequential
+
+from .utils.geo import get_mean_sd_by_band
+from .utils.geo import merge_tiles
+from .utils.torchgeo import NoBordersGridGeoSampler
+
+def get_model_size(model):
+ torch.save(model.state_dict(), "temp.p")
+ size = os.path.getsize("temp.p")/1e6
+ os.remove('temp.p')
+ return size
+
+class EncoderAlgorithm(QgsProcessingAlgorithm):
+ """
+ """
+
+ FEAT_OPTION= 'FEAT_OPTION'
+ INPUT = 'INPUT'
+ CKPT = 'CKPT'
+ BANDS = 'BANDS'
+ STRIDE = 'STRIDE'
+ SIZE = 'SIZE'
+ EXTENT = 'EXTENT'
+ QUANT = 'QUANT'
+ OUTPUT = 'OUTPUT'
+ RESOLUTION = 'RESOLUTION'
+ CRS = 'CRS'
+ CUDA = 'CUDA'
+ BATCH_SIZE = 'BATCH_SIZE'
+ CUDA_ID = 'CUDA_ID'
+ BACKBONE_CHOICE = 'BACKBONE_CHOICE'
+ MERGE_METHOD = 'MERGE_METHOD'
+ WORKERS = 'WORKERS'
+ PAUSES = 'PAUSES'
+
+
+ def initAlgorithm(self, config=None):
+ """
+ Here we define the inputs and output of the algorithm, along
+ with some other properties.
+ """
+ cwd = Path(__file__).parent.absolute()
+
+ self.addParameter(
+ QgsProcessingParameterRasterLayer(
+ name=self.INPUT,
+ description=self.tr(
+ 'Input raster layer or image file path'),
+ defaultValue=os.path.join(cwd,'rasters','test.tif'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBand(
+ name=self.BANDS,
+ description=self.tr('Selected Bands (defaults to all bands selected)'),
+ defaultValue = None,
+ parentLayerParameterName=self.INPUT,
+ optional=True,
+ allowMultiple=True,
+ )
+ )
+
+ crs_param = QgsProcessingParameterCrs(
+ name=self.CRS,
+ description=self.tr('Target CRS (default to original CRS)'),
+ optional=True,
+ )
+
+ res_param = QgsProcessingParameterNumber(
+ name=self.RESOLUTION,
+ description=self.tr(
+ 'Target resolution in meters (default to native resolution)'),
+ type=QgsProcessingParameterNumber.Double,
+ optional=True,
+ minValue=0,
+ maxValue=100000
+ )
+
+ cuda_id_param = QgsProcessingParameterNumber(
+ name=self.CUDA_ID,
+ description=self.tr(
+ 'CUDA Device ID (choose which GPU to use, default to device 0)'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=0,
+ minValue=0,
+ maxValue=9
+ )
+ nworkers_param = QgsProcessingParameterNumber(
+ name=self.WORKERS,
+ description=self.tr(
+ 'Number of CPU workers for dataloader (0 selects all)'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=0,
+ minValue=0,
+ maxValue=10
+ )
+ pauses_param = QgsProcessingParameterNumber(
+ name=self.PAUSES,
+ description=self.tr(
+ 'Schedule pauses between batches to ease CPU usage (in seconds).'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=0,
+ minValue=0,
+ maxValue=10000
+ )
+
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ name=self.EXTENT,
+ description=self.tr(
+ 'Processing extent (default to the entire image)'),
+ optional=True
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name=self.SIZE,
+ description=self.tr(
+ 'Sampling size (the raster will be sampled in a square with a side of that many pixel)'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue = 224,
+ minValue=1,
+ maxValue=1024
+ )
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name=self.STRIDE,
+ description=self.tr(
+ 'Stride (If smaller than the sampling size, tiles will overlap. If larger, it may cause errors.)'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue = 224,
+ minValue=1,
+ maxValue=1024
+ )
+ )
+
+ chkpt_param = QgsProcessingParameterFile(
+ name=self.CKPT,
+ description=self.tr(
+ 'Pretrained checkpoint'),
+ extension='pth',
+ optional=True,
+ defaultValue=None
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT,
+ self.tr(
+ "Output directory (choose the location that the image features will be saved)"),
+ defaultValue=os.path.join(cwd,'features'),
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ self.CUDA,
+ self.tr("Use GPU if CUDA is available."),
+ defaultValue=True
+ )
+ )
+ self.addParameter (
+ QgsProcessingParameterString(
+ name = self.BACKBONE_CHOICE,
+ description = self.tr(
+ 'Backbone choice (see huggingface.co/timm/)'),
+ defaultValue = 'vit_base_patch16_224.dino',
+ # defaultValue = 'vit_small_patch16_224.dino',
+ )
+ )
+
+
+
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ self.FEAT_OPTION,
+ self.tr("Display features map"),
+ defaultValue=True
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name=self.BATCH_SIZE,
+ # large images will be sampled into patches in a grid-like fashion
+ description=self.tr(
+ 'Batch size (take effect if choose to use GPU and CUDA is available)'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=1,
+ minValue=1,
+ maxValue=1024
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ self.QUANT,
+ self.tr("Quantization of the model to reduce space"),
+ defaultValue=True
+ )
+ )
+
+ self.merge_options = ['first', 'min', 'max','average','sum', 'count', 'last']
+ merge_param = QgsProcessingParameterEnum(
+ name=self.MERGE_METHOD,
+ description=self.tr(
+ 'Merge method at the end of inference.'),
+ options=self.merge_options,
+ defaultValue=0,
+ )
+
+ for param in (
+ crs_param,
+ res_param,
+ chkpt_param,
+ cuda_id_param,
+ merge_param,
+ nworkers_param,
+ pauses_param
+ ):
+ param.setFlags(
+ param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(param)
+
+
+
+ @torch.no_grad()
+ def processAlgorithm(self, parameters, context, feedback):
+ """
+ Here is where the processing itself takes place.
+ """
+ self.process_options(parameters, context, feedback)
+
+ RasterDataset.filename_glob = self.rlayer_name
+ RasterDataset.all_bands = [
+ self.rlayer.bandName(i_band) for i_band in range(1, self.rlayer.bandCount()+1)
+ ]
+ # currently only support rgb bands
+ input_bands = [self.rlayer.bandName(i_band)
+ for i_band in self.selected_bands]
+
+ feedback.pushInfo(f'create dataset')
+ if self.crs == self.rlayer.crs():
+ dataset = RasterDataset(
+ paths=self.rlayer_dir, crs=None, res=self.res, bands=input_bands, cache=False)
+ else:
+ dataset = RasterDataset(
+ paths=self.rlayer_dir, crs=self.crs.toWkt(), res=self.res, bands=input_bands, cache=False)
+ extent_bbox = BoundingBox(minx=self.extent.xMinimum(), maxx=self.extent.xMaximum(), miny=self.extent.yMinimum(), maxy=self.extent.yMaximum(),
+ mint=dataset.index.bounds[4], maxt=dataset.index.bounds[5])
+
+
+ feedback.pushInfo(f'create model')
+ model = timm.create_model(
+ self.backbone_name,
+ pretrained=True,
+ in_chans=len(input_bands),
+ num_classes=0
+ )
+
+ feedback.pushInfo(f'model done')
+ data_config = timm.data.resolve_model_data_config(model)
+ _, h, w, = data_config['input_size']
+
+ if self.quantization:
+ feedback.pushInfo(f'before quantization : {get_model_size(model)}')
+
+ model = torch.quantization.quantize_dynamic(
+ model, {nn.Linear}, dtype=torch.qint8
+ )
+ feedback.pushInfo(f'after quantization : {get_model_size(model)}')
+
+ transform = AugmentationSequential(
+ T.ConvertImageDtype(torch.float32), # change dtype for normalize to be possible
+ K.Normalize(self.means,self.sds), # normalize occurs only on raster, not mask
+ K.Resize((h, w)), # resize to 224*224 pixels, regardless of sampling size
+ data_keys=["image"],
+ )
+ dataset.transforms = transform
+
+
+ # sampler = GridGeoSampler(
+ # dataset,
+ # size=self.size,
+ # stride=self.stride,
+ # roi=extent_bbox,
+ # units=Units.PIXELS
+ # ) # Units.CRS or Units.PIXELS
+ sampler = NoBordersGridGeoSampler(
+ dataset,
+ size=self.size,
+ stride=self.stride,
+ roi=extent_bbox,
+ units=Units.PIXELS
+ ) # Units.CRS or Units.PIXELS
+
+ if len(sampler) == 0:
+ self.load_feature = False
+ feedback.pushWarning(f'\n !!!No available patch sample inside the chosen extent!!! \n')
+
+ if torch.cuda.is_available() and self.use_gpu:
+ if self.cuda_id + 1 > torch.cuda.device_count():
+ self.cuda_id = torch.cuda.device_count() - 1
+ cuda_device = f'cuda:{self.cuda_id}'
+ device = f'cuda:{self.cuda_id}'
+ else:
+ self.batch_size = 1
+ device = 'cpu'
+
+ feedback.pushInfo(f'Device id: {device}')
+
+ feedback.pushInfo(f'model to dedvice')
+ model.to(device=device)
+
+ feedback.pushInfo(f'Batch size: {self.batch_size}')
+ dataloader = DataLoader(
+ dataset,
+ batch_size=self.batch_size,
+ sampler=sampler,
+ collate_fn=stack_samples,
+ num_workers=self.nworkers,
+ )
+
+ feedback.pushInfo(f'Patch sample num: {len(sampler)}')
+ feedback.pushInfo(f'Total batch num: {len(dataloader)}')
+ feedback.pushInfo(f'\n\n{"-"*16}\nBegining inference \n{"-"*16}\n\n')
+
+ ## compute parameters hash to have a unique identifier for the run
+ ## some parameters do not change the encoding part of the algorithm
+ keys_to_remove = ['MERGE_METHOD', 'WORKERS', 'PAUSES']
+ param_encoder = {key: parameters[key] for key in parameters if key not in keys_to_remove}
+
+ param_hash = hashlib.md5(str(param_encoder).encode("utf-8")).hexdigest()
+ output_subdir = os.path.join(self.output_dir,param_hash)
+ output_subdir = Path(output_subdir)
+ output_subdir.mkdir(parents=True, exist_ok=True)
+ self.output_subdir = output_subdir
+ feedback.pushInfo(f'output_subdir: {output_subdir}')
+
+
+ last_batch_done = self.get_last_batch_done()
+ if last_batch_done >= 0:
+ feedback.pushInfo(f"\n\n {'-'*8} \n Resuming at batch number {last_batch_done}\n {'-'*8} \n\n")
+
+ bboxes = [] # keep track of bboxes to have coordinates at the end
+ elapsed_time_list = []
+ total = 100 / len(dataloader) if len(dataloader) else 0
+
+ for current, sample in enumerate(dataloader):
+
+ if current <= last_batch_done:
+ continue
+
+ start_time = time.time()
+ # Stop the algorithm if cancel button has been clicked
+ if feedback.isCanceled():
+ self.load_feature = False
+ feedback.pushWarning(
+ self.tr("\n !!!Processing is canceled by user!!! \n"))
+ break
+ feedback.pushInfo(f'\n{"-"*8}\nBatch no. {current} loaded')
+
+ images = sample['image'].to(device)
+ if len(images.shape) > 4:
+ images = images.squeeze(1)
+ feedback.pushInfo(f'Batch shape {images.shape}')
+
+ features = model.forward_features(images)
+ features = features[:,1:,:] # take only patch tokens
+ if current <= last_batch_done + 1:
+ n_patches = int(np.sqrt(features.shape[1]))
+ features = features.view(features.shape[0],n_patches,n_patches,features.shape[-1])
+ features = features.detach().cpu().numpy()
+ feedback.pushInfo(f'Features shape {features.shape}')
+ self.save_features(features,sample['bbox'], current)
+ feedback.pushInfo(f'Features saved')
+
+ bboxes.extend(sample['bbox'])
+
+ if self.pauses != 0:
+ time.sleep(self.pauses)
+
+ end_time = time.time()
+ # get the execution time of encoder, ms
+ elapsed_time = (end_time - start_time)
+ elapsed_time_list.append(elapsed_time)
+ time_spent = sum(elapsed_time_list)
+ time_remain = (time_spent / (current + 1)) * \
+ (len(dataloader) - current - 1)
+
+ # TODO: show gpu usage info
+ # if torch.cuda.is_available() and self.use_gpu:
+ # gpu_mem_used = torch.cuda.max_memory_reserved(self.sam_model.device) / (1024 ** 3)
+ # # gpu_mem_free = torch.cuda.mem_get_info(self.sam_model.device)[0] / (1024 ** 3)
+ # gpu_mem_total = torch.cuda.mem_get_info(self.sam_model.device)[1] / (1024 ** 3)
+ # feedback.pushInfo(
+ # f'GPU memory usage: {gpu_mem_used:.2f}GB / {gpu_mem_total:.2f}GB')
+ # feedback.pushInfo(str(torch.cuda.memory_summary(self.sam_model.device)))
+
+ feedback.pushInfo(f"Encoder executed with {elapsed_time:.3f}s")
+ feedback.pushInfo(f"Time spent: {time_spent:.3f}s")
+
+ if time_remain <= 60:
+ feedback.pushInfo(f"Estimated time remaining: {time_remain:.3f}s \n {'-'*8}")
+ else:
+ time_remain_m, time_remain_s = divmod(int(time_remain), 60)
+ time_remain_h, time_remain_m = divmod(time_remain_m, 60)
+ feedback.pushInfo(f"Estimated time remaining: {time_remain_h:d}h:{time_remain_m:02d}m:{time_remain_s:02d}s \n" )
+
+ # Update the progress bar
+ feedback.setProgress(int((current+1) * total))
+
+
+ all_tiles = [os.path.join(self.output_subdir,f) for f in os.listdir(self.output_subdir) if f.endswith('.tif')]
+ dst_path = Path(os.path.join(self.output_subdir,'merged.tiff'))
+ feedback.pushInfo(f"\n\n{'-'*8}\n Merging tiles \n{'-'*8}\n" )
+
+ merge_tiles(
+ tiles = all_tiles,
+ dst_path = dst_path,
+ method = self.merge_method,
+ )
+
+ parameters['OUTPUT_RASTER']=dst_path
+
+ return {"Output feature path": self.output_subdir, 'Patch samples saved': self.iPatch, 'OUTPUT_RASTER':dst_path}
+
+ def get_last_batch_done(self):
+
+ ## get largest batch_number achieved
+ ## files are saved with the pattern '{batch_number}_{image_id_within_batch}.tif'
+ # Regular expression pattern to extract numbers
+ pattern = re.compile(r'^(\d+)_\d+\.tif$')
+
+ # Initialize a set to store unique first numbers
+ batch_numbers = set()
+
+ # Iterate over all files in the directory
+ for filename in os.listdir(self.output_subdir):
+ # Match the filename pattern
+ match = pattern.match(filename)
+ if match:
+ # Extract the batch number
+ batch_number = int(match.group(1))
+ # Add to the set of batch numbers
+ batch_numbers.add(batch_number)
+
+ # Find the maximum value of the batch numbers
+ if batch_numbers:
+ return max(batch_numbers)
+ else:
+ return -1
+
+
+ def save_features(
+ self,
+ feature: np.ndarray,
+ bboxes: BoundingBox,
+ nbatch: int,
+ ):
+
+ # iterate over batch_size dimension
+ for idx in range(feature.shape[0]):
+ _, height, width, channels = feature.shape
+ bbox = bboxes[idx]
+ rio_transform = rasterio.transform.from_bounds(bbox.minx, bbox.miny, bbox.maxx, bbox.maxy, width, height) # west, south, east, north, width, height
+ feature_path = os.path.join(self.output_subdir, f"{nbatch}_{idx}.tif")
+ with rasterio.open(
+ feature_path,
+ mode="w",
+ driver="GTiff",
+ height=height,
+ width=width,
+ count=channels,
+ dtype='float32',
+ crs=self.crs.toWkt(),
+ transform=rio_transform
+ ) as ds:
+ ds.write(np.transpose(feature[idx, ...], (2, 0, 1)))
+ tags = {
+ "model_type": self.backbone_name,
+ }
+ ds.update_tags(**tags)
+
+ self.iPatch += 1
+
+ return
+
+ def process_options(self,parameters, context, feedback):
+ self.iPatch = 0
+
+ self.feature_dir = ""
+
+ self.FEAT_OPTION = self.parameterAsBoolean(
+ parameters, self.FEAT_OPTION, context)
+
+ feedback.pushInfo(
+ f'PARAMETERS :\n{parameters}')
+
+ feedback.pushInfo(
+ f'CONTEXT :\n{context}')
+
+ feedback.pushInfo(
+ f'FEEDBACK :\n{feedback}')
+
+ rlayer = self.parameterAsRasterLayer(
+ parameters, self.INPUT, context)
+
+ if rlayer is None:
+ raise QgsProcessingException(
+ self.invalidRasterError(parameters, self.INPUT))
+
+ self.selected_bands = self.parameterAsInts(
+ parameters, self.BANDS, context)
+
+ if len(self.selected_bands) == 0:
+ self.selected_bands = list(range(1, rlayer.bandCount()+1))
+
+ if max(self.selected_bands) > rlayer.bandCount():
+ raise QgsProcessingException(
+ self.tr("The chosen bands exceed the largest band number!")
+ )
+
+ ckpt_path = self.parameterAsFile(
+ parameters, self.CKPT, context)
+
+ self.backbone_name = self.parameterAsString(
+ parameters, self.BACKBONE_CHOICE, context)
+
+ self.stride = self.parameterAsInt(
+ parameters, self.STRIDE, context)
+ self.size = self.parameterAsInt(
+ parameters, self.SIZE, context)
+ res = self.parameterAsDouble(
+ parameters, self.RESOLUTION, context)
+ crs = self.parameterAsCrs(
+ parameters, self.CRS, context)
+ extent = self.parameterAsExtent(
+ parameters, self.EXTENT, context)
+ self.quantization = self.parameterAsBoolean(
+ parameters, self.QUANT, context)
+ self.use_gpu = self.parameterAsBoolean(
+ parameters, self.CUDA, context)
+ self.batch_size = self.parameterAsInt(
+ parameters, self.BATCH_SIZE, context)
+ self.output_dir = self.parameterAsString(
+ parameters, self.OUTPUT, context)
+ self.cuda_id = self.parameterAsInt(
+ parameters, self.CUDA_ID, context)
+ self.pauses = self.parameterAsInt(
+ parameters, self.PAUSES, context)
+ self.nworkers = self.parameterAsInt(
+ parameters, self.WORKERS, context)
+ merge_method_idx = self.parameterAsEnum(
+ parameters, self.MERGE_METHOD, context)
+ self.merge_method = self.merge_options[merge_method_idx]
+
+ rlayer_data_provider = rlayer.dataProvider()
+
+ # handle crs
+ if crs is None or not crs.isValid():
+ crs = rlayer.crs()
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}') # 0 for meters, 6 for degrees, 9 for unknown
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ if crs.mapUnits() == Qgis.DistanceUnit.Degrees:
+ crs = self.estimate_utm_crs(rlayer.extent())
+
+ # target crs should use meters as units
+ if crs.mapUnits() != Qgis.DistanceUnit.Meters:
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}')
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ raise QgsProcessingException(
+ self.tr("Only support CRS with the units as meters")
+ )
+
+ # 0 for meters, 6 for degrees, 9 for unknown
+ UNIT_METERS = 0
+ UNIT_DEGREES = 6
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ layer_units = 'degrees'
+ else:
+ layer_units = 'meters'
+ # if res is not provided, get res info from rlayer
+ if np.isnan(res) or res == 0:
+ res = rlayer.rasterUnitsPerPixelX() # rasterUnitsPerPixelY() is negative
+ target_units = layer_units
+ else:
+ # when given res in meters by users, convert crs to utm if the original crs unit is degree
+ if crs.mapUnits() != UNIT_METERS: # Qgis.DistanceUnit.Meters:
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ # estimate utm crs based on layer extent
+ crs = self.estimate_utm_crs(rlayer.extent())
+ else:
+ raise QgsProcessingException(
+ f"Resampling of image with the CRS of {crs.authid()} in meters is not supported.")
+ target_units = 'meters'
+ # else:
+ # res = (rlayer_extent.xMaximum() -
+ # rlayer_extent.xMinimum()) / rlayer.width()
+ self.res = res
+
+ # handle extent
+ if extent.isNull():
+ extent = rlayer.extent() # QgsProcessingUtils.combineLayerExtents(layers, crs, context)
+ extent_crs = rlayer.crs()
+ else:
+ if extent.isEmpty():
+ raise QgsProcessingException(
+ self.tr("The extent for processing can not be empty!"))
+ extent_crs = self.parameterAsExtentCrs(
+ parameters, self.EXTENT, context)
+ # if extent crs != target crs, convert it to target crs
+ if extent_crs != crs:
+ transform = QgsCoordinateTransform(
+ extent_crs, crs, context.transformContext())
+ # extent = transform.transformBoundingBox(extent)
+ # to ensure coverage of the transformed extent
+ # convert extent to polygon, transform polygon, then get boundingBox of the new polygon
+ extent_polygon = QgsGeometry.fromRect(extent)
+ extent_polygon.transform(transform)
+ extent = extent_polygon.boundingBox()
+ extent_crs = crs
+
+ # check intersects between extent and rlayer_extent
+ if rlayer.crs() != crs:
+ transform = QgsCoordinateTransform(
+ rlayer.crs(), crs, context.transformContext())
+ rlayer_extent = transform.transformBoundingBox(
+ rlayer.extent())
+ else:
+ rlayer_extent = rlayer.extent()
+ if not rlayer_extent.intersects(extent):
+ raise QgsProcessingException(
+ self.tr("The extent for processing is not intersected with the input image!"))
+
+ feedback.pushInfo(f'backbne type : {self.backbone_name}')
+
+ img_width_in_extent = round(
+ (extent.xMaximum() - extent.xMinimum())/self.res)
+ img_height_in_extent = round(
+ (extent.yMaximum() - extent.yMinimum())/self.res)
+
+ # Send some information to the user
+ feedback.pushInfo(
+ f'Layer path: {rlayer_data_provider.dataSourceUri()}')
+ # feedback.pushInfo(
+ # f'Layer band scale: {rlayer_data_provider.bandScale(self.selected_bands[0])}')
+ feedback.pushInfo(f'Layer name: {rlayer.name()}')
+ if rlayer.crs().authid():
+ feedback.pushInfo(f'Layer CRS: {rlayer.crs().authid()}')
+ else:
+ feedback.pushInfo(
+ f'Layer CRS in WKT format: {rlayer.crs().toWkt()}')
+ feedback.pushInfo(
+ f'Layer pixel size: {rlayer.rasterUnitsPerPixelX()}, {rlayer.rasterUnitsPerPixelY()} {layer_units}')
+
+ feedback.pushInfo(f'Bands selected: {self.selected_bands}')
+
+ if crs.authid():
+ feedback.pushInfo(f'Target CRS: {crs.authid()}')
+ else:
+ feedback.pushInfo(f'Target CRS in WKT format: {crs.toWkt()}')
+ # feedback.pushInfo('Band number is {}'.format(rlayer.bandCount()))
+ # feedback.pushInfo('Band name is {}'.format(rlayer.bandName(1)))
+ feedback.pushInfo(f'Target resolution: {self.res} {target_units}')
+ # feedback.pushInfo('Layer display band name is {}'.format(
+ # rlayer.dataProvider().displayBandName(1)))
+ feedback.pushInfo(
+ (f'Processing extent: minx:{extent.xMinimum():.6f}, maxx:{extent.xMaximum():.6f},'
+ f'miny:{extent.yMinimum():.6f}, maxy:{extent.yMaximum():.6f}'))
+ feedback.pushInfo(
+ (f'Processing image size: (width {img_width_in_extent}, '
+ f'height {img_height_in_extent})'))
+
+ # feedback.pushInfo(
+ # f'SAM Image Size: {self.sam_model.image_encoder.img_size}')
+
+ self.rlayer_path = rlayer.dataProvider().dataSourceUri()
+ self.rlayer_dir = os.path.dirname(self.rlayer_path)
+ self.rlayer_name = os.path.basename(self.rlayer_path)
+
+ # get mean and sd of dataset from raster metadata
+ means, sds = get_mean_sd_by_band(self.rlayer_path)
+ # subset with selected_bands
+ feedback.pushInfo(f'Selected bands: {self.selected_bands}')
+ self.means = [means[i-1] for i in self.selected_bands]
+ self.sds = [sds[i-1] for i in self.selected_bands]
+ feedback.pushInfo(f'Means for normalization: {self.means}')
+ feedback.pushInfo(f'Std. dev. for normalization: {self.sds}')
+
+ ## passing parameters to self once everything has been processed
+ self.extent = extent
+ self.rlayer = rlayer
+ self.crs = crs
+
+
+ # used to handle any thread-sensitive cleanup which is required by the algorithm.
+ def postProcessAlgorithm(self, context, feedback) -> Dict[str, Any]:
+ return {}
+
+
+ def tr(self, string):
+ """
+ Returns a translatable string with the self.tr() function.
+ """
+ return QCoreApplication.translate('Processing', string)
+
+ def createInstance(self):
+ return EncoderAlgorithm()
+
+ def name(self):
+ """
+ Returns the algorithm name, used for identifying the algorithm. This
+ string should be fixed for the algorithm, and must not be localised.
+ The name should be unique within each provider. Names should contain
+ lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return 'encoder'
+
+ def displayName(self):
+ """
+ Returns the translated algorithm name, which should be used for any
+ user-visible display of the algorithm name.
+ """
+ return self.tr('Image Encoder')
+
+ def group(self):
+ """
+ Returns the name of the group this algorithm belongs to. This string
+ should be localised.
+ """
+ return self.tr('')
+
+ def groupId(self):
+ """
+ Returns the unique ID of the group this algorithm belongs to. This
+ string should be fixed for the algorithm, and must not be localised.
+ The group id should be unique within each provider. Group id should
+ contain lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return ''
+
+ def shortHelpString(self):
+ """
+ Returns a localised short helper string for the algorithm. This string
+ should provide a basic description about what the algorithm does and the
+ parameters and outputs associated with it..
+ """
+ return self.tr("Generate image features using a deep learning backbone.")
+
+ def icon(self):
+ return 'E'
+
diff --git a/src/iamap/iamap.py b/src/iamap/iamap.py
new file mode 100644
index 0000000..a0a36a0
--- /dev/null
+++ b/src/iamap/iamap.py
@@ -0,0 +1,234 @@
+import processing
+from PyQt5.QtWidgets import (
+ QAction,
+ QToolBar,
+ QApplication,
+ QDialog
+)
+from PyQt5.QtCore import pyqtSignal, QObject
+from qgis.core import QgsApplication
+from qgis.gui import QgisInterface
+from .provider import IAMapProvider
+from .icons import QIcon_EncoderTool, QIcon_ReductionTool, QIcon_ClusterTool, QIcon_SimilarityTool, QIcon_RandomforestTool
+
+
+class IAMap(QObject):
+ execute_iamap = pyqtSignal()
+
+ def __init__(self, iface: QgisInterface, cwd: str):
+ super().__init__()
+ self.iface = iface
+ self.cwd = cwd
+
+ def initProcessing(self):
+ self.provider = IAMapProvider()
+ QgsApplication.processingRegistry().addProvider(self.provider)
+
+ def initGui(self):
+ self.initProcessing()
+
+ self.toolbar: QToolBar = self.iface.addToolBar('IAMap Toolbar')
+ self.toolbar.setObjectName('IAMapToolbar')
+ self.toolbar.setToolTip('IAMap Toolbar')
+
+ self.actionEncoder = QAction(
+ QIcon_EncoderTool,
+ "Deep Learning Image Encoder",
+ self.iface.mainWindow()
+ )
+ self.actionReducer = QAction(
+ QIcon_ReductionTool,
+ "Reduce dimensions",
+ self.iface.mainWindow()
+ )
+ self.actionCluster = QAction(
+ QIcon_ClusterTool,
+ "Cluster raster",
+ self.iface.mainWindow()
+ )
+ self.actionSimilarity = QAction(
+ QIcon_SimilarityTool,
+ "Compute similarity",
+ self.iface.mainWindow()
+ )
+ self.actionRF = QAction(
+ QIcon_RandomforestTool,
+ "Use Random Forest algorithm",
+ self.iface.mainWindow()
+ )
+ self.actionEncoder.setObjectName("mActionEncoder")
+ self.actionReducer.setObjectName("mActionReducer")
+ self.actionCluster.setObjectName("mActionCluster")
+ self.actionSimilarity.setObjectName("mactionSimilarity")
+ self.actionRF.setObjectName("mactionRF")
+
+ self.actionEncoder.setToolTip(
+ "Encode a raster with a deep learning backbone")
+ self.actionReducer.setToolTip(
+ "Reduce raster dimensions")
+ self.actionCluster.setToolTip(
+ "Cluster raster")
+ self.actionSimilarity.setToolTip(
+ "Compute similarity")
+ self.actionRF.setToolTip(
+ "Use Random Forest ")
+
+ self.actionEncoder.triggered.connect(self.encodeImage)
+ self.actionReducer.triggered.connect(self.reduceImage)
+ self.actionCluster.triggered.connect(self.clusterImage)
+ self.actionSimilarity.triggered.connect(self.similarityImage)
+ self.actionRF.triggered.connect(self.rfImage)
+
+ self.toolbar.addAction(self.actionEncoder)
+ self.toolbar.addAction(self.actionReducer)
+ self.toolbar.addAction(self.actionCluster)
+ self.toolbar.addAction(self.actionSimilarity)
+ self.toolbar.addAction(self.actionRF)
+
+ def unload(self):
+ # self.wdg_select.setVisible(False)
+ self.iface.removeToolBarIcon(self.actionEncoder)
+ self.iface.removeToolBarIcon(self.actionReducer)
+ self.iface.removeToolBarIcon(self.actionCluster)
+ self.iface.removeToolBarIcon(self.actionSimilarity)
+ self.iface.removeToolBarIcon(self.actionRF)
+
+ del self.actionEncoder
+ del self.actionReducer
+ del self.actionCluster
+ del self.actionSimilarity
+ del self.actionRF
+ del self.toolbar
+ QgsApplication.processingRegistry().removeProvider(self.provider)
+
+ def encodeImage(self):
+ '''
+ '''
+ result = processing.execAlgorithmDialog('iamap:encoder', {})
+ print(result)
+ # Check if algorithm execution was successful
+ if result:
+ # Retrieve output parameters from the result dictionary
+ if 'OUTPUT_RASTER' in result:
+ output_raster_path = result['OUTPUT_RASTER']
+
+ # Add the output raster layer to the map canvas
+ self.iface.addRasterLayer(str(output_raster_path), 'merged features')
+ else:
+ # Handle missing or unexpected output
+ print('Output raster not found in algorithm result.')
+ else:
+ # Handle algorithm execution failure or cancellation
+ print('Algorithm execution was not successful.')
+ # processing.execAlgorithmDialog('', {})
+ # self.close_all_dialogs()
+
+
+ def reduceImage(self):
+ '''
+ '''
+ result = processing.execAlgorithmDialog('iamap:reduction', {})
+ print(result)
+ # Check if algorithm execution was successful
+ if result:
+ # Retrieve output parameters from the result dictionary
+ if 'OUTPUT_RASTER' in result:
+ output_raster_path = result['OUTPUT_RASTER']
+
+ # Add the output raster layer to the map canvas
+ self.iface.addRasterLayer(str(output_raster_path), 'reduced features')
+ else:
+ # Handle missing or unexpected output
+ print('Output raster not found in algorithm result.')
+ else:
+ # Handle algorithm execution failure or cancellation
+ print('Algorithm execution was not successful.')
+ # processing.execAlgorithmDialog('', {})
+
+
+ def clusterImage(self):
+ '''
+ '''
+ result = processing.execAlgorithmDialog('iamap:cluster', {})
+ print(result)
+ # Check if algorithm execution was successful
+ if result:
+ # Retrieve output parameters from the result dictionary
+ if 'OUTPUT_RASTER' in result:
+ output_raster_path = result['OUTPUT_RASTER']
+
+ # Add the output raster layer to the map canvas
+ self.iface.addRasterLayer(str(output_raster_path), 'clustering')
+ else:
+ # Handle missing or unexpected output
+ print('Output raster not found in algorithm result.')
+ else:
+ # Handle algorithm execution failure or cancellation
+ print('Algorithm execution was not successful.')
+ # processing.execAlgorithmDialog('', {})
+
+
+ def similarityImage(self):
+ '''
+ '''
+ result = processing.execAlgorithmDialog('iamap:similarity', {})
+ print(result)
+ # Check if algorithm execution was successful
+ if result:
+ # Retrieve output parameters from the result dictionary
+ if 'OUTPUT_RASTER' in result:
+ output_raster_path = result['OUTPUT_RASTER']
+
+ # Add the output raster layer to the map canvas
+ self.iface.addRasterLayer(str(output_raster_path), 'similarity map')
+ else:
+ # Handle missing or unexpected output
+ print('Output raster not found in algorithm result.')
+ else:
+ # Handle algorithm execution failure or cancellation
+ print('Algorithm execution was not successful.')
+ # processing.execAlgorithmDialog('', {})
+
+ def rfImage(self):
+ '''
+ '''
+ result = processing.execAlgorithmDialog('iamap:Random_forest', {})
+ print(result)
+ # Check if algorithm execution was successful
+ if result:
+ # Retrieve output parameters from the result dictionary
+ if 'OUTPUT_RASTER' in result:
+ output_raster_path = result['OUTPUT_RASTER']
+
+ # Add the output raster layer to the map canvas
+ self.iface.addRasterLayer(str(output_raster_path), 'random forest map')
+ else:
+ # Handle missing or unexpected output
+ print('Output raster not found in algorithm result.')
+ else:
+ # Handle algorithm execution failure or cancellation
+ print('Algorithm execution was not successful.')
+ # processing.execAlgorithmDialog('', {})
+
+
+
+ def close_all_dialogs(self):
+ # Get the main QGIS window (QgisInterface)
+ qgis_main_window = self.iface.mainWindow()
+
+ # Get all open dialogs associated with the main window
+ open_dialogs = qgis_main_window.findChildren(QDialog)
+
+ # Iterate through the open dialogs and close them
+ for dialog in open_dialogs:
+ # Check if the dialog is visible (to avoid closing hidden dialogs)
+ if dialog.isVisible():
+ # Close the dialog
+ dialog.close()
+
+
+
+
+
+
+
diff --git a/src/iamap/icons/__init__.py b/src/iamap/icons/__init__.py
new file mode 100644
index 0000000..32f705a
--- /dev/null
+++ b/src/iamap/icons/__init__.py
@@ -0,0 +1,16 @@
+import os
+from PyQt5.QtGui import QIcon
+
+cwd = os.path.abspath(os.path.dirname(__file__))
+# encoder_tool_path = os.path.join(cwd, 'encoder_tool.svg')
+encoder_tool_path = os.path.join(cwd, 'network-base-solid-svgrepo-com.svg')
+reduction_tool_path = os.path.join(cwd, 'chart-scatterplot-solid-svgrepo-com.svg')
+cluster_tool_path = os.path.join(cwd, 'chart-scatter-3d-svgrepo-com.svg')
+similarity_tool_path = os.path.join(cwd, 'scatter-plot-svgrepo-com.svg')
+random_forest_tool_path = os.path.join(cwd, 'forest-svgrepo-com.svg')
+
+QIcon_EncoderTool = QIcon(encoder_tool_path)
+QIcon_ReductionTool = QIcon(reduction_tool_path)
+QIcon_ClusterTool = QIcon(cluster_tool_path)
+QIcon_SimilarityTool = QIcon(similarity_tool_path)
+QIcon_RandomforestTool = QIcon(random_forest_tool_path)
diff --git a/src/iamap/icons/chart-scatter-3d-svgrepo-com.svg b/src/iamap/icons/chart-scatter-3d-svgrepo-com.svg
new file mode 100644
index 0000000..45be800
--- /dev/null
+++ b/src/iamap/icons/chart-scatter-3d-svgrepo-com.svg
@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="utf-8"?><!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg width="800px" height="800px" viewBox="0 0 24 24" fill="none" xmlns="http://www.w3.org/2000/svg">
+<path d="M12 4V14M12 14L4 20M12 14L20 20M12 20H12.01M4 13H4.01M17 13H17.01M8 9H8.01M20 9H20.01M5 5H5.01M18 5H18.01" stroke="#000000" stroke-width="2" stroke-linecap="round" stroke-linejoin="round"/>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/chart-scatterplot-solid-svgrepo-com.svg b/src/iamap/icons/chart-scatterplot-solid-svgrepo-com.svg
new file mode 100644
index 0000000..756b3b4
--- /dev/null
+++ b/src/iamap/icons/chart-scatterplot-solid-svgrepo-com.svg
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+
<!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg width="800px" height="800px" viewBox="0 0 48 48" xmlns="http://www.w3.org/2000/svg">
+ <title>chart-data-sets-solid</title>
+ <g id="Layer_2" data-name="Layer 2">
+ <g id="invisible_box" data-name="invisible box">
+ <rect width="48" height="48" fill="none"/>
+ </g>
+ <g id="icons_Q2" data-name="icons Q2">
+ <circle cx="14" cy="12" r="2"/>
+ <circle cx="14" cy="22" r="2"/>
+ <circle cx="19" cy="17" r="2"/>
+ <circle cx="36" cy="34" r="2"/>
+ <circle cx="36" cy="24" r="2"/>
+ <circle cx="31" cy="29" r="2"/>
+ <path d="M41.7,20A2,2,0,0,0,44,18V6a2,2,0,0,0-2-2H30a2,2,0,0,0-2,2.3A2.1,2.1,0,0,0,30.1,8h7.1L8,37.2V6A2,2,0,0,0,4,6V44H42a2,2,0,0,0,0-4H10.8L40,10.8v7.1A2.1,2.1,0,0,0,41.7,20Z"/>
+ </g>
+ </g>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/encoder_tool.svg b/src/iamap/icons/encoder_tool.svg
new file mode 100644
index 0000000..4cb16a1
--- /dev/null
+++ b/src/iamap/icons/encoder_tool.svg
@@ -0,0 +1,92 @@
+<?xml version="1.0" encoding="utf-8"?>
+<!-- Generator: Adobe Illustrator 25.4.1, SVG Export Plug-In . SVG Version: 6.00 Build 0) -->
+<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
+ viewBox="0 0 64 64" style="enable-background:new 0 0 64 64;" xml:space="preserve">
+<style type="text/css">
+ .st0{fill:#FF5252;stroke:#000000;stroke-width:6.390949e-04;stroke-miterlimit:10;}
+ .st1{fill:#1AA1DB;stroke:#000000;stroke-width:6.390949e-04;stroke-miterlimit:10;}
+ .st2{fill:#91C965;}
+ .st3{fill:#88BA5D;}
+ .st4{fill:#9BDC6E;}
+ .st5{fill:#3C8763;}
+</style>
+<g id="net" xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" xmlns:svg="http://www.w3.org/2000/svg">
+ <g id="g16">
+ <path id="path4" class="st0" d="M5.5,24.7c-0.6,0-1.1-0.5-1.1-1.1v-4.3c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v4.3
+ C6.6,24.2,6.1,24.7,5.5,24.7 M5.5,34.3c-0.6,0-1.1-0.5-1.1-1.1V30c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v3.2
+ C6.6,33.8,6.1,34.3,5.5,34.3 M5.5,44.9c-0.6,0-1.1-0.5-1.1-1.1v-6.4c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v6.4
+ C6.6,44.4,6.1,44.9,5.5,44.9"/>
+ <path id="path6" class="st0" d="M13.6,53.1c-0.2,0-0.3,0-0.5-0.1l-3.9-1.8c-0.5-0.3-0.8-0.9-0.5-1.4s0.9-0.8,1.4-0.5l3.9,1.8
+ c0.5,0.3,0.8,0.9,0.5,1.4C14.4,52.9,14,53.1,13.6,53.1 M21.3,56.8c-0.2,0-0.3,0-0.5-0.1L18,55.3c-0.5-0.3-0.8-0.9-0.5-1.4
+ s0.9-0.8,1.4-0.5l2.9,1.4c0.5,0.3,0.8,0.9,0.5,1.4C22.1,56.5,21.7,56.8,21.3,56.8 M27.9,59.8c-0.2,0-0.3,0-0.5-0.1l-2.7-1.3
+ c-0.5-0.3-0.8-0.9-0.5-1.4s0.9-0.8,1.4-0.5l2.7,1.3c0.5,0.3,0.8,0.9,0.5,1.4C28.7,59.6,28.3,59.8,27.9,59.8"/>
+ <path id="path8" class="st0" d="M58.8,24.7c-0.6,0-1.1-0.5-1.1-1.1v-4.3c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v4.3
+ C59.8,24.2,59.4,24.7,58.8,24.7 M58.8,31.1c-0.6,0-1.1-0.5-1.1-1.1v-3.2c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1V30
+ C59.8,30.6,59.4,31.1,58.8,31.1 M58.8,38.5c-0.6,0-1.1-0.5-1.1-1.1v-3.2c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v3.2
+ C59.8,38.1,59.4,38.5,58.8,38.5 M58.8,44.9c-0.6,0-1.1-0.5-1.1-1.1v-3.2c0-0.6,0.5-1.1,1.1-1.1s1.1,0.5,1.1,1.1v3.2
+ C59.8,44.4,59.4,44.9,58.8,44.9"/>
+ <path id="path10" class="st0" d="M50.6,53.1c-0.4,0-0.8-0.2-1-0.6c-0.3-0.5,0-1.2,0.5-1.4l3.9-1.8c0.5-0.3,1.2,0,1.4,0.5
+ c0.3,0.5,0,1.2-0.5,1.4L51.1,53C51,53.1,50.8,53.1,50.6,53.1 M42.9,56.8c-0.4,0-0.8-0.2-1-0.6c-0.3-0.5,0-1.2,0.5-1.4l2.9-1.4
+ c0.5-0.2,1.2,0,1.4,0.5c0.3,0.5,0,1.2-0.5,1.4l-2.9,1.4C43.2,56.7,43.1,56.8,42.9,56.8 M36.4,59.8c-0.4,0-0.8-0.2-1-0.6
+ c-0.3-0.5,0-1.2,0.5-1.4l2.7-1.3c0.5-0.2,1.2,0,1.4,0.5s0,1.2-0.5,1.4l-2.7,1.3C36.7,59.8,36.5,59.8,36.4,59.8"/>
+ <path id="path12" class="st0" d="M38.7,6.5c-0.1,0-0.3,0-0.4-0.1l-2.3-1c-0.5-0.2-0.8-0.8-0.6-1.4c0.2-0.5,0.8-0.8,1.4-0.6l2.3,1
+ c0.5,0.2,0.8,0.8,0.6,1.4C39.6,6.2,39.2,6.5,38.7,6.5 M45.6,9.3c-0.1,0-0.3,0-0.4-0.1l-3-1.2c-0.5-0.2-0.8-0.8-0.6-1.4
+ C41.8,6.1,42.4,5.8,43,6l3,1.2c0.5,0.2,0.8,0.8,0.6,1.4C46.5,9.1,46.1,9.3,45.6,9.3 M54.5,13c-0.1,0-0.3,0-0.4-0.1l-3.9-1.6
+ c-0.5-0.2-0.8-0.8-0.6-1.4c0.2-0.5,0.8-0.8,1.4-0.6l3.9,1.6c0.5,0.2,0.8,0.8,0.6,1.4C55.3,12.7,54.9,13,54.5,13"/>
+ <path id="path14" class="st0" d="M25.5,6.5c-0.4,0-0.8-0.2-1-0.7s0-1.2,0.6-1.4l2.3-1c0.5-0.2,1.2,0,1.4,0.6
+ c0.2,0.5,0,1.2-0.6,1.4l-2.3,1C25.8,6.4,25.7,6.5,25.5,6.5 M18.6,9.3c-0.4,0-0.8-0.2-1-0.7s0-1.2,0.6-1.4l3-1.2
+ c0.5-0.2,1.2,0,1.4,0.6c0.2,0.5,0,1.2-0.6,1.4l-3,1.2C18.9,9.3,18.8,9.3,18.6,9.3 M9.8,13c-0.4,0-0.8-0.2-1-0.7s0-1.2,0.6-1.4
+ l3.9-1.6c0.5-0.2,1.2,0,1.4,0.6c0.2,0.5,0,1.2-0.6,1.4l-3.9,1.6C10,12.9,9.9,13,9.8,13"/>
+ </g>
+ <g id="g32">
+ <path id="path24" class="st0" d="M16.2,23.6c-0.2,0-0.4-0.1-0.6-0.2l-6.4-4.3c-0.5-0.3-0.6-1-0.3-1.5s1-0.6,1.5-0.3l6.4,4.3
+ c0.5,0.3,0.6,1,0.3,1.5C16.8,23.4,16.5,23.6,16.2,23.6"/>
+ <path id="path26" class="st0" d="M48.1,23.6c-0.3,0-0.7-0.2-0.9-0.5c-0.3-0.5-0.2-1.2,0.3-1.5l6.4-4.3c0.5-0.3,1.2-0.2,1.5,0.3
+ s0.2,1.2-0.3,1.5l-6.4,4.3C48.5,23.6,48.3,23.6,48.1,23.6"/>
+ <path id="path28" class="st0" d="M54.5,46c-0.2,0-0.4-0.1-0.6-0.2l-6.4-4.3c-0.5-0.3-0.6-1-0.3-1.5s1-0.6,1.5-0.3l6.4,4.3
+ c0.5,0.3,0.6,1,0.3,1.5C55.2,45.8,54.8,46,54.5,46"/>
+ <path id="path30" class="st0" d="M9.8,46c-0.3,0-0.7-0.2-0.9-0.5C8.6,45,8.7,44.4,9.2,44l6.4-4.3c0.5-0.3,1.2-0.2,1.5,0.3
+ s0.2,1.2-0.3,1.5l-6.4,4.3C10.2,45.9,10,46,9.8,46"/>
+ </g>
+ <g id="g46">
+ <path id="path34" class="st1" d="M5.5,20.4c-2.9,0-5.3-2.4-5.3-5.3s2.4-5.3,5.3-5.3s5.3,2.4,5.3,5.3S8.4,20.4,5.5,20.4"/>
+ <path id="path36" class="st1" d="M32.1,10.8c-2.9,0-5.3-2.4-5.3-5.3s2.4-5.3,5.3-5.3s5.3,2.4,5.3,5.3S35.1,10.8,32.1,10.8"/>
+ <path id="path38" class="st1" d="M32.1,64.1c-2.9,0-5.3-2.4-5.3-5.3c0-2.9,2.4-5.3,5.3-5.3s5.3,2.4,5.3,5.3S35.1,64.1,32.1,64.1"
+ />
+ <path id="path40" class="st1" d="M58.8,20.4c-2.9,0-5.3-2.4-5.3-5.3s2.4-5.3,5.3-5.3c2.9,0,5.3,2.4,5.3,5.3S61.7,20.4,58.8,20.4"
+ />
+ <path id="path42" class="st1" d="M58.8,53.4c-2.9,0-5.3-2.4-5.3-5.3c0-2.9,2.4-5.3,5.3-5.3c2.9,0,5.3,2.4,5.3,5.3
+ S61.7,53.4,58.8,53.4"/>
+ <path id="path44" class="st1" d="M5.5,53.4c-2.9,0-5.3-2.4-5.3-5.3s2.4-5.3,5.3-5.3s5.3,2.4,5.3,5.3S8.4,53.4,5.5,53.4"/>
+ </g>
+</g>
+<g id="cube">
+ <g>
+ <polygon id="polygon18" class="st2" points="16.7,41.6 16.7,22.6 32.2,30.7 32.2,49.5 "/>
+ <polygon id="polygon20" class="st3" points="47.8,41.6 47.8,22.6 32.2,30.7 32.2,49.6 "/>
+ <polygon id="polygon22" class="st4" points="16.7,22.6 32.5,14.6 47.8,22.6 32.2,30.7 "/>
+ <path id="path4_00000163069267588332021340000012919285466319501952_" class="st5" d="M48.4,22.7C48.4,22.6,48.4,22.6,48.4,22.7
+ C48.4,22.6,48.4,22.6,48.4,22.7c0-0.1,0-0.1,0-0.2c0,0,0,0,0,0c0,0,0,0,0-0.1c0,0,0,0,0,0c0,0,0,0,0-0.1c0,0,0,0,0,0
+ c0,0,0,0,0-0.1c0,0,0,0,0,0c0,0,0,0-0.1-0.1c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0-0.1,0c0,0,0,0,0,0l-15.5-8
+ c-0.2-0.1-0.4-0.1-0.6,0l-15.5,8c0,0,0,0,0,0c0,0,0,0-0.1,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0.1c0,0,0,0,0,0
+ c0,0,0,0,0,0.1c0,0,0,0,0,0c0,0,0,0,0,0.1c0,0,0,0,0,0c0,0,0,0,0,0.1c0,0,0,0,0,0c0,0,0,0,0,0.1c0,0,0,0,0,0c0,0,0,0.1,0,0.1v18.8
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0l5.2,2.6c0,0,0,0,0.1,0l5.1,2.6c0,0,0.1,0.1,0.2,0.1L32,50c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0l5.1-2.6c0.1,0,0.1,0,0.2-0.1l5.1-2.6
+ c0,0,0,0,0.1,0l5.2-2.6c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0
+ c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0c0,0,0,0,0,0L48.4,22.7L48.4,22.7z M27.1,18l3.9,2l-3.9,2l-3.9-2
+ L27.1,18z M32.2,20.7l3.9,2l-3.9,2l-3.9-2L32.2,20.7z M41.3,20l-3.9,2l-3.9-2l3.9-2L41.3,20z M17.3,30.2l4,2v4.4l-4-2V30.2z
+ M22.5,32.9l4,2v4.4l-4-2V32.9z M31.6,36.2l-4-2V29l4,2V36.2z M32.8,31l4-2v5.2l-4,2V31z M38,28.3l4-2v5.2l-4,2V28.3z M32.2,29.9
+ l-3.8-2l3.8-2l3.8,2L32.2,29.9z M26.5,33.5l-4-2v-5.2l4,2V33.5z M27.7,35.5l4,2V42l-4-2V35.5z M32.8,37.5l4-2v4.4l-4,2V37.5z
+ M38,34.9l4-2v4.4l-4,2V34.9z M43.2,32.2l4-2v4.4l-4,2V32.2z M47.2,28.8l-4,2v-5.2l4-2V28.8z M37.4,27.3l-3.8-2l3.9-2l3.8,2
+ L37.4,27.3z M30.9,25.3l-3.8,2l-3.9-2l3.8-2L30.9,25.3z M21.3,30.9l-4-2v-5.2l4,2V30.9z M17.3,36l4,2v5.1l-4-2V36z M22.5,38.6l4,2
+ v5.1l-4-2V38.6z M27.7,41.3l4,2v5.1l-4-2V41.3z M32.8,43.3l4-2v5.1l-4,2V43.3z M38,40.7l4-2v5.1l-4,2V40.7z M43.2,38l4-2v5.1l-4,2
+ V38z M42.6,24.6l-3.8-2l3.9-2l3.8,2L42.6,24.6z M32.2,15.4l3.9,2l-3.9,2l-3.9-2L32.2,15.4z M21.9,20.7l3.9,2l-3.8,2l-3.9-2
+ L21.9,20.7z"/>
+ </g>
+</g>
+</svg>
diff --git a/src/iamap/icons/forest-svgrepo-com.svg b/src/iamap/icons/forest-svgrepo-com.svg
new file mode 100644
index 0000000..eeefe73
--- /dev/null
+++ b/src/iamap/icons/forest-svgrepo-com.svg
@@ -0,0 +1,26 @@
+<?xml version="1.0" encoding="utf-8"?>
+
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg height="800px" width="800px" version="1.1" id="_x32_" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink"
+ viewBox="0 0 512 512" xml:space="preserve">
+<style type="text/css">
+ .st0{fill:#000000;}
+</style>
+<g>
+ <path class="st0" d="M346.483,226.653c-58.176-75.765-90.498-181.813-90.498-181.813s-32.318,106.048-90.505,181.813
+ c0,0,26.66,16.09,41.21,7.569c0,0-14.55,65.341-79.995,151.514c58.176,18.923,101.81-12.328,101.81-12.328v93.75h21.025h12.916
+ h21.021v-93.75c0,0,43.642,31.25,101.817,12.328c-65.457-86.174-79.995-151.514-79.995-151.514
+ C319.826,242.743,346.483,226.653,346.483,226.653z"/>
+ <path class="st0" d="M160.886,307.087c-19.185-35.761-24.363-59.015-24.363-59.015c8.768,5.141,23.33-1.454,29.058-4.376
+ c1.522-0.84,2.417-1.379,2.417-1.379c-5.313-6.985-10.353-14.276-15.186-21.718c-34.855-54.482-53.972-117.26-53.972-117.26
+ s-24.711,81.041-69.23,138.977c0,0,20.361,12.283,31.542,5.756c0,0-11.181,49.956-61.151,115.88
+ c44.451,14.426,77.788-9.443,77.788-9.443v71.674h42.034v-71.674c0,0,3.035,2.151,8.415,4.759
+ C141.633,340.391,152.332,322.817,160.886,307.087z"/>
+ <path class="st0" d="M450.849,248.071c11.121,6.527,31.474-5.756,31.474-5.756c-44.454-57.936-69.155-138.977-69.155-138.977
+ s-19.125,62.778-54.05,117.26c-4.766,7.441-9.803,14.733-15.123,21.718c0,0,0.906,0.54,2.428,1.379
+ c5.725,2.922,20.29,9.517,29.058,4.376c0,0-5.178,23.328-24.442,59.09c8.566,15.655,19.331,33.303,32.723,52.106
+ c5.381-2.608,8.423-4.759,8.423-4.759v71.674h41.967v-71.674c0,0,33.394,23.869,77.848,9.443
+ C461.97,298.027,450.849,248.071,450.849,248.071z"/>
+</g>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/network-base-solid-svgrepo-com.svg b/src/iamap/icons/network-base-solid-svgrepo-com.svg
new file mode 100644
index 0000000..a93ee40
--- /dev/null
+++ b/src/iamap/icons/network-base-solid-svgrepo-com.svg
@@ -0,0 +1,22 @@
+<?xml version="1.0" encoding="utf-8"?>
+
<!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg width="800px" height="800px" viewBox="0 0 48 48" xmlns="http://www.w3.org/2000/svg">
+ <g id="Layer_2" data-name="Layer 2">
+ <g id="invisible_box" data-name="invisible box">
+ <rect width="48" height="48" fill="none"/>
+ </g>
+ <g id="Layer_7" data-name="Layer 7">
+ <g>
+ <path d="M25.5,25.7l-1.5,1-1.5-1a2.2,2.2,0,0,0-2.8.7,2.1,2.1,0,0,0,.7,2.8L23,30.7a1.8,1.8,0,0,0,2,0l2.6-1.5a2.1,2.1,0,0,0,.7-2.8A2.2,2.2,0,0,0,25.5,25.7Z"/>
+ <path d="M3,18.7l2.5,1.5a2.3,2.3,0,0,0,1.1.3,2.1,2.1,0,0,0,1.7-.9A2,2,0,0,0,7.8,17a2,2,0,0,0,.5-2.6,2.2,2.2,0,0,0-2.8-.7L3,15.3A1.9,1.9,0,0,0,2,17,2.1,2.1,0,0,0,3,18.7Z"/>
+ <path d="M12,14.2a1.9,1.9,0,0,0,1-.3l4.1-2.4a2,2,0,0,0,.6-2.7A1.9,1.9,0,0,0,15,8.1l-4.1,2.4a2,2,0,0,0-.6,2.7A1.8,1.8,0,0,0,12,14.2Z"/>
+ <path d="M22.5,8.3l1.5-1,1.5,1a2.2,2.2,0,0,0,1.1.2,2,2,0,0,0,1.7-.9,2.1,2.1,0,0,0-.7-2.8L25,3.3a1.8,1.8,0,0,0-2,0L20.4,4.8a2.1,2.1,0,0,0-.7,2.8A2.2,2.2,0,0,0,22.5,8.3Z"/>
+ <path d="M30.9,11.5,35,13.9a1.9,1.9,0,0,0,1,.3,1.8,1.8,0,0,0,1.7-1,2,2,0,0,0-.6-2.7L33,8.1a1.9,1.9,0,0,0-2.7.7A2,2,0,0,0,30.9,11.5Z"/>
+ <path d="M40.2,17a2,2,0,0,0-.5,2.6,2.1,2.1,0,0,0,1.7.9,2.3,2.3,0,0,0,1.1-.3L45,18.7A2.1,2.1,0,0,0,46,17a1.9,1.9,0,0,0-1-1.7l-2.5-1.6a2.2,2.2,0,0,0-2.8.7A2,2,0,0,0,40.2,17Z"/>
+ <path d="M10.9,22.5h0L15,24.9a1.9,1.9,0,0,0,1,.3,1.8,1.8,0,0,0,1.7-1,2,2,0,0,0-.6-2.7L13,19.1a1.9,1.9,0,0,0-2.7.6A2.1,2.1,0,0,0,10.9,22.5Z"/>
+ <path d="M32,25.2a1.9,1.9,0,0,0,1-.3l4.1-2.4a2.1,2.1,0,0,0,.6-2.8,1.9,1.9,0,0,0-2.7-.6l-4.1,2.4a2,2,0,0,0-.6,2.7A1.8,1.8,0,0,0,32,25.2Z"/>
+ <path d="M45,28.3l-6.2-3.8L24,33.2,9.2,24.5,3,28.3A1.9,1.9,0,0,0,2,30a2.1,2.1,0,0,0,1,1.7l20,12a1.8,1.8,0,0,0,2,0l20-12A2.1,2.1,0,0,0,46,30,1.9,1.9,0,0,0,45,28.3Z"/>
+ </g>
+ </g>
+ </g>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/network-platform-svgrepo-com.svg b/src/iamap/icons/network-platform-svgrepo-com.svg
new file mode 100644
index 0000000..2d1a1de
--- /dev/null
+++ b/src/iamap/icons/network-platform-svgrepo-com.svg
@@ -0,0 +1,15 @@
+<?xml version="1.0" encoding="utf-8"?>
+
<!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg width="800px" height="800px" viewBox="0 0 48 48" xmlns="http://www.w3.org/2000/svg">
+ <title>network-platform</title>
+ <g id="Layer_2" data-name="Layer 2">
+ <g id="invisible_box" data-name="invisible box">
+ <rect width="48" height="48" fill="none"/>
+ <rect width="48" height="48" fill="none"/>
+ <rect width="48" height="48" fill="none"/>
+ </g>
+ <g id="icons_Q2" data-name="icons Q2">
+ <path d="M40,18a6.2,6.2,0,0,0-5.7,4H25.8a9.9,9.9,0,0,0-1.5-3.3l3.5-4.1A7.8,7.8,0,0,0,30,15a6,6,0,1,0-6-6,6.1,6.1,0,0,0,.8,3l-3.6,4.1A8.5,8.5,0,0,0,17,15l-1.7.2-1.2-2.8A6,6,0,1,0,10,14h.4l1.3,2.8a8.9,8.9,0,0,0,0,14.4L10.4,34H10a6,6,0,1,0,4.1,1.6l1.2-2.8L17,33a8.5,8.5,0,0,0,4.2-1.1L24.8,36a6.1,6.1,0,0,0-.8,3,6,6,0,1,0,6-6,7.8,7.8,0,0,0-2.2.4l-3.5-4.1A9.9,9.9,0,0,0,25.8,26h8.5A6.2,6.2,0,0,0,40,30a6,6,0,0,0,0-12ZM8,8a2,2,0,1,1,2,2A2,2,0,0,1,8,8Zm2,34a2,2,0,1,1,2-2A2,2,0,0,1,10,42ZM30,7a2,2,0,1,1-2,2A2,2,0,0,1,30,7ZM12,24a5,5,0,1,1,5,5A5,5,0,0,1,12,24ZM32,39a2,2,0,1,1-2-2A2,2,0,0,1,32,39Zm8-13a2,2,0,1,1,2-2A2,2,0,0,1,40,26Z"/>
+ </g>
+ </g>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/scatter-plot-outline-alerted-svgrepo-com.svg b/src/iamap/icons/scatter-plot-outline-alerted-svgrepo-com.svg
new file mode 100644
index 0000000..22ffd24
--- /dev/null
+++ b/src/iamap/icons/scatter-plot-outline-alerted-svgrepo-com.svg
@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="utf-8"?><!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg fill="#000000" width="800px" height="800px" viewBox="0 0 36 36" version="1.1" preserveAspectRatio="xMidYMid meet" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ <title>scatter-plot-outline-alerted</title>
+ <path class="clr-i-outline--alerted clr-i-outline-path-1--alerted" d="M 34 29 C 34 30.105 33.105 31 32 31 L 4 31 C 2.895 31 2 30.105 2 29 L 2 7 C 2 5.895 2.895 5 4 5 L 21.958 5 L 20.786 7 L 4 7 L 4 29 L 32 29 L 32 15.357 L 34 15.357 Z"></path><path class="clr-i-outline--alerted clr-i-outline-path-2--alerted" d="M 9.101 15.8 C 9.413 16.111 9.919 16.111 10.231 15.8 L 11.391 14.64 L 12.551 15.8 C 12.964 16.256 13.717 16.094 13.905 15.507 C 14.002 15.208 13.914 14.881 13.681 14.67 L 12.531 13.54 L 13.691 12.38 C 14.147 11.966 13.985 11.214 13.399 11.025 C 13.1 10.929 12.772 11.017 12.561 11.25 L 11.401 12.41 L 10.231 11.22 C 9.817 10.763 9.065 10.926 8.877 11.512 C 8.78 11.811 8.868 12.139 9.101 12.35 L 10.261 13.54 L 9.101 14.67 C 8.789 14.982 8.789 15.487 9.101 15.8 Z"></path><path class="clr-i-outline--alerted clr-i-outline-path-3--alerted" d="M 15.176 25.536 C 15.488 25.847 15.994 25.847 16.306 25.536 L 17.466 24.376 L 18.626 25.536 C 19.039 25.992 19.792 25.83 19.98 25.243 C 20.077 24.944 19.989 24.617 19.756 24.406 L 18.606 23.276 L 19.766 22.116 C 20.222 21.702 20.06 20.95 19.474 20.761 C 19.175 20.665 18.847 20.753 18.636 20.986 L 17.476 22.146 L 16.306 20.956 C 15.892 20.499 15.14 20.662 14.952 21.248 C 14.855 21.547 14.943 21.875 15.176 22.086 L 16.336 23.276 L 15.176 24.406 C 14.864 24.718 14.864 25.223 15.176 25.536 Z"></path><path class="clr-i-outline--alerted clr-i-outline-path-4--alerted" d="M 22.912 20.343 C 23.224 20.654 23.73 20.654 24.042 20.343 L 25.202 19.183 L 26.362 20.343 C 26.775 20.799 27.528 20.637 27.716 20.05 C 27.813 19.751 27.725 19.424 27.492 19.213 L 26.342 18.083 L 27.502 16.923 C 27.958 16.509 27.796 15.757 27.21 15.568 C 26.911 15.472 26.583 15.56 26.372 15.793 L 25.212 16.953 L 24.042 15.763 C 23.628 15.306 22.876 15.469 22.688 16.055 C 22.591 16.354 22.679 16.682 22.912 16.893 L 24.072 18.083 L 22.912 19.213 C 22.6 19.525 22.6 20.03 22.912 20.343 Z"></path><path class="clr-i-outline--alerted clr-i-outline-path-5--alerted clr-i-alert" d="M 26.854 1.144 L 21.134 11.004 C 20.579 11.818 21.114 12.928 22.097 13.001 C 22.142 13.005 22.188 13.006 22.234 13.004 L 33.684 13.004 C 34.669 13.036 35.319 11.991 34.855 11.122 C 34.834 11.081 34.81 11.042 34.784 11.004 L 29.064 1.144 C 28.57 0.299 27.348 0.299 26.854 1.144 Z"></path>
+ <rect x="0" y="0" width="36" height="36" fill-opacity="0"/>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/icons/scatter-plot-svgrepo-com.svg b/src/iamap/icons/scatter-plot-svgrepo-com.svg
new file mode 100644
index 0000000..543d6cd
--- /dev/null
+++ b/src/iamap/icons/scatter-plot-svgrepo-com.svg
@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="utf-8"?><!-- Uploaded to: SVG Repo, www.svgrepo.com, Generator: SVG Repo Mixer Tools -->
+<svg fill="#000000" width="800px" height="800px" viewBox="0 0 36 36" version="1.1" preserveAspectRatio="xMidYMid meet" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ <title>scatter-plot-line</title>
+ <path class="clr-i-outline clr-i-outline-path-1" d="M 32 5 L 4 5 C 2.895 5 2 5.895 2 7 L 2 29 C 2 30.105 2.895 31 4 31 L 32 31 C 33.105 31 34 30.105 34 29 L 34 7 C 34 5.895 33.105 5 32 5 Z M 4 29 L 4 7 L 32 7 L 32 29 Z"></path><path class="clr-i-outline clr-i-outline-path-2" d="M 9.101 15.8 C 9.413 16.111 9.919 16.111 10.231 15.8 L 11.391 14.64 L 12.551 15.8 C 12.964 16.256 13.717 16.094 13.905 15.507 C 14.002 15.208 13.914 14.881 13.681 14.67 L 12.531 13.54 L 13.691 12.38 C 14.147 11.966 13.985 11.214 13.399 11.025 C 13.1 10.929 12.772 11.017 12.561 11.25 L 11.401 12.41 L 10.231 11.22 C 9.817 10.763 9.065 10.926 8.877 11.512 C 8.78 11.811 8.868 12.139 9.101 12.35 L 10.261 13.54 L 9.101 14.67 C 8.789 14.982 8.789 15.487 9.101 15.8 Z"></path><path class="clr-i-outline clr-i-outline-path-3" d="M 15.176 25.536 C 15.488 25.847 15.994 25.847 16.306 25.536 L 17.466 24.376 L 18.626 25.536 C 19.039 25.992 19.792 25.83 19.98 25.243 C 20.077 24.944 19.989 24.617 19.756 24.406 L 18.606 23.276 L 19.766 22.116 C 20.222 21.702 20.06 20.95 19.474 20.761 C 19.175 20.665 18.847 20.753 18.636 20.986 L 17.476 22.146 L 16.306 20.956 C 15.892 20.499 15.14 20.662 14.952 21.248 C 14.855 21.547 14.943 21.875 15.176 22.086 L 16.336 23.276 L 15.176 24.406 C 14.864 24.718 14.864 25.223 15.176 25.536 Z"></path><path class="clr-i-outline clr-i-outline-path-4" d="M 22.912 20.343 C 23.224 20.654 23.73 20.654 24.042 20.343 L 25.202 19.183 L 26.362 20.343 C 26.775 20.799 27.528 20.637 27.716 20.05 C 27.813 19.751 27.725 19.424 27.492 19.213 L 26.342 18.083 L 27.502 16.923 C 27.958 16.509 27.796 15.757 27.21 15.568 C 26.911 15.472 26.583 15.56 26.372 15.793 L 25.212 16.953 L 24.042 15.763 C 23.628 15.306 22.876 15.469 22.688 16.055 C 22.591 16.354 22.679 16.682 22.912 16.893 L 24.072 18.083 L 22.912 19.213 C 22.6 19.525 22.6 20.03 22.912 20.343 Z"></path>
+ <rect x="0" y="0" width="36" height="36" fill-opacity="0"/>
+</svg>
\ No newline at end of file
diff --git a/src/iamap/provider.py b/src/iamap/provider.py
new file mode 100644
index 0000000..5d7797b
--- /dev/null
+++ b/src/iamap/provider.py
@@ -0,0 +1,46 @@
+from qgis.core import QgsProcessingProvider
+
+from .encoder import EncoderAlgorithm
+from .reduction import ReductionAlgorithm
+from .clustering import ClusterAlgorithm
+from .similarity import SimilarityAlgorithm
+from .random_forest import RFAlgorithm
+from .icons import QIcon_EncoderTool
+
+
+class IAMapProvider(QgsProcessingProvider):
+
+ def loadAlgorithms(self, *args, **kwargs):
+ self.addAlgorithm(EncoderAlgorithm())
+ self.addAlgorithm(ReductionAlgorithm())
+ self.addAlgorithm(ClusterAlgorithm())
+ self.addAlgorithm(SimilarityAlgorithm())
+ self.addAlgorithm(RFAlgorithm())
+ # add additional algorithms here
+ # self.addAlgorithm(MyOtherAlgorithm())
+
+ def id(self, *args, **kwargs):
+ """The ID of your plugin, used for identifying the provider.
+
+ This string should be a unique, short, character only string,
+ eg "qgis" or "gdal". This string should not be localised.
+ """
+ return 'iamap'
+
+ def name(self, *args, **kwargs):
+ """The human friendly name of your plugin in Processing.
+
+ This string should be as short as possible (e.g. "Lastools", not
+ "Lastools version 1.0.1 64-bit") and localised.
+ """
+ return self.tr('IAMap')
+
+ def icon(self):
+ """Should return a QIcon which is used for your provider inside
+ the Processing toolbox.
+ """
+ return QIcon_EncoderTool
+
+ def longName(self) -> str:
+ return self.name()
+
diff --git a/src/iamap/random_forest.py b/src/iamap/random_forest.py
new file mode 100644
index 0000000..c71d334
--- /dev/null
+++ b/src/iamap/random_forest.py
@@ -0,0 +1,599 @@
+import os
+import numpy as np
+from pathlib import Path
+from typing import Dict, Any
+import joblib
+
+import rasterio
+from rasterio import windows
+import geopandas as gpd
+import pandas as pd
+from shapely.geometry import box
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (Qgis,
+ QgsGeometry,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterEnum,
+ QgsCoordinateTransform,
+ QgsProcessingException,
+ QgsProcessingAlgorithm,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterString,
+ QgsProcessingParameterBand,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterCrs,
+ QgsProcessingParameterDefinition,
+ )
+import torch
+import torch.nn as nn
+from sklearn.ensemble import RandomForestClassifier
+from sklearn.metrics import accuracy_score
+from sklearn.preprocessing import LabelEncoder
+from sklearn.model_selection import train_test_split
+
+import json
+
+class RFAlgorithm(QgsProcessingAlgorithm):
+ """
+ """
+
+ INPUT = 'INPUT'
+ BANDS = 'BANDS'
+ EXTENT = 'EXTENT'
+ LOAD = 'LOAD'
+ OUTPUT = 'OUTPUT'
+ RESOLUTION = 'RESOLUTION'
+ CRS = 'CRS'
+ TEMPLATE = 'TEMPLATE'
+ COLONNE_RF = 'COLONNE_RF'
+ TEMPLATE_TEST = 'TEMPLATE_TEST'
+
+ def initAlgorithm(self, config=None):
+ """
+ Here we define the inputs and output of the algorithm, along
+ with some other properties.
+ """
+ cwd = Path(__file__).parent.absolute()
+
+ self.addParameter(
+ QgsProcessingParameterRasterLayer(
+ name=self.INPUT,
+ description=self.tr(
+ 'Input raster layer or image file path'),
+ defaultValue=os.path.join(cwd,'rasters','test.tif'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBand(
+ name=self.BANDS,
+ description=self.tr(
+ 'Selected Bands (defaults to all bands selected)'),
+ defaultValue=None,
+ parentLayerParameterName=self.INPUT,
+ optional=True,
+ allowMultiple=True,
+ )
+ )
+
+ crs_param = QgsProcessingParameterCrs(
+ name=self.CRS,
+ description=self.tr('Target CRS (default to original CRS)'),
+ optional=True,
+ )
+
+ res_param = QgsProcessingParameterNumber(
+ name=self.RESOLUTION,
+ description=self.tr(
+ 'Target resolution in meters (default to native resolution)'),
+ type=QgsProcessingParameterNumber.Double,
+ optional=True,
+ minValue=0,
+ maxValue=100000
+ )
+
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ name=self.EXTENT,
+ description=self.tr(
+ 'Processing extent (default to the entire image)'),
+ optional=True
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterFile(
+ name=self.TEMPLATE,
+ description=self.tr(
+ 'Input shapefile path for training data set for random forest (if no test data_set, will be devised in train and test)'),
+ defaultValue=os.path.join(cwd,'rasters','RF.gpkg'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterFile(
+ name=self.TEMPLATE_TEST,
+ description=self.tr(
+ 'Input shapefile path for test data set for random forest (optional)'),
+ optional = True
+ ),
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT,
+ self.tr(
+ "Output directory (choose the location that the image features will be saved)"),
+ defaultValue=os.path.join(cwd,'features'),
+ )
+ )
+
+
+ self.addParameter (
+ QgsProcessingParameterString(
+ name = self.COLONNE_RF,
+ description = self.tr(
+ 'Name of the column you want random forest to work on'),
+ defaultValue = 'Type',
+ # defaultValue = 'vit_small_patch16_224.dino',
+ )
+ )
+
+
+ for param in (crs_param, res_param):
+ param.setFlags(
+ param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(param)
+
+
+
+ def processAlgorithm(self, parameters, context, feedback):
+ """
+ Here is where the processing itself takes place.
+ """
+ self.process_options(parameters, context, feedback)
+
+ gdf = gpd.read_file(self.template)
+ feedback.pushInfo(f"self_template_test : {self.template_test}")
+ DATA_SET_TEST=False
+ if(self.template_test != ''):
+ gdf_test = gpd.read_file(self.template_test)
+ gdf_test = gdf_test.to_crs(self.crs.toWkt())
+ DATA_SET_TEST = True
+ feedback.pushInfo(f"In good loop !")
+
+ gdf = gdf.to_crs(self.crs.toWkt())
+
+ feedback.pushInfo(f"DATA_SET_TEST : {DATA_SET_TEST}")
+
+
+ feedback.pushInfo(f'before extent: {len(gdf)}')
+ bounds = box(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ )
+ feedback.pushInfo(f'xmin: {self.extent.xMinimum()},ymin: {self.extent.yMinimum()}, xmax: {self.extent.xMaximum()}, ymax: {self.extent.yMaximum()} ')
+ gdf = gdf[gdf.within(bounds)]
+ feedback.pushInfo(f'after extent: {len(gdf)}')
+
+ if len(gdf) == 0:
+ feedback.pushWarning("No template points within extent !")
+ return False
+
+
+
+
+ input_bands = [i_band -1 for i_band in self.selected_bands]
+
+
+ with rasterio.open(self.rlayer_path) as ds:
+
+
+
+
+ gdf = gdf.to_crs(ds.crs)
+
+ pixel_values_test = []
+
+ pixel_values = []
+
+ transform = ds.transform
+ crs = ds.crs
+ win = windows.from_bounds(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ transform=transform
+ )
+ raster = ds.read(window=win)
+ transform = ds.window_transform(win)
+ raster = raster[input_bands,:,:]
+
+ if (DATA_SET_TEST == True):
+ gdf_test = gdf_test.to_crs(ds.crs)
+ for index, data in gdf_test.iterrows():
+ # Get the coordinates of the point in the raster's pixel space
+ x, y = data.geometry.x, data.geometry.y
+ feedback.pushInfo (f"x : {x}, y : {y}")
+ feedback.pushInfo (f"gdf geometry : {gdf['geometry']}")
+
+ # Convert point coordinates to pixel coordinates within the window
+ col, row = ~transform * (x, y) # Convert from map coordinates to pixel coordinates
+ col, row = int(col), int(row)
+ feedback.pushInfo(f'after extent: {row, col}')
+ pixel_values_test.append(list(raster[:,row, col]))
+ template_npy_test = np.asarray(pixel_values_test)
+ template_test = torch.from_numpy(template_npy_test).to(torch.float32)
+
+ for index, data in gdf.iterrows():
+ # Get the coordinates of the point in the raster's pixel space
+ x, y = data.geometry.x, data.geometry.y
+ feedback.pushInfo (f"x : {x}, y : {y}")
+ feedback.pushInfo (f"gdf geometry : {gdf['geometry']}")
+
+ # Convert point coordinates to pixel coordinates within the window
+ col, row = ~transform * (x, y) # Convert from map coordinates to pixel coordinates
+ col, row = int(col), int(row)
+ feedback.pushInfo(f'after extent: {row, col}')
+ pixel_values.append(list(raster[:,row, col]))
+
+ raster = np.transpose(raster, (1,2,0))
+
+ feedback.pushInfo(f'{raster.shape}')
+
+ template_npy = np.asarray(pixel_values)
+
+ feedback.pushInfo(f'points : {template_npy}')
+ feedback.pushInfo(f'dim points : {template_npy.shape}')
+ template = torch.from_numpy(template_npy).to(torch.float32)
+
+
+
+
+ feat_img = torch.from_numpy(raster)
+
+ #template contient les valeurs
+ #y = gdf['Type']
+ #y=gdf ['Desc_']
+
+ if (DATA_SET_TEST == False):
+
+ if self.colonne_rf in gdf.columns :
+
+ y = gdf[self.colonne_rf]
+ else :
+ feedback.pushWarning (f'{self.colonne_rf} is not a valid column name of the dataset !!')
+
+
+
+ X_train, X_test, y_train, y_test = train_test_split(template, y, test_size=0.4, random_state=55)
+
+ rf_classifier = RandomForestClassifier(n_estimators=100, min_samples_split=4, random_state=42)
+ rf_classifier.fit(X_train, y_train)
+
+ if (DATA_SET_TEST == True):
+
+ y_train=gdf[self.colonne_rf]
+ y_test=gdf_test[self.colonne_rf]
+
+ X_train = template
+
+
+ X_test = template_test
+
+
+ rf_classifier = RandomForestClassifier(n_estimators=100, min_samples_split=4, random_state=42)
+
+ rf_classifier.fit(X_train, y_train)
+
+ params = rf_classifier.get_params()
+
+ #joblib.dump(rf_classifier, model_file)
+
+ y_pred = rf_classifier.predict(X_test)
+ accuracy = accuracy_score(y_test, y_pred)
+ feedback.pushInfo(f"Accuracy ; {accuracy}")
+
+ predicted_types = rf_classifier.predict(raster.reshape(-1, raster.shape[-1]))
+ feedback.pushInfo(f"predicted types ; {predicted_types.shape}")
+ predicted_types_image = predicted_types.reshape(raster.shape[:-1])
+ feedback.pushInfo(f"predicted_types_image ; {predicted_types_image.shape}")
+
+ label_encoder = LabelEncoder()
+ predicted_types_numeric = label_encoder.fit_transform(predicted_types_image.flatten())
+ feedback.pushInfo(f'carte avant transfo : {predicted_types_numeric.shape}')
+ predicted_types_numeric = predicted_types_numeric.reshape(predicted_types_image.shape)
+ feedback.pushInfo(f'carte après transfo : {predicted_types_numeric.shape}')
+
+
+ height, width = predicted_types_numeric.shape
+ channels = 1
+
+
+
+
+
+ dst_path = os.path.join(self.output_dir,'random_forest.tif')
+ params_file = os.path.join(self.output_dir, 'random_forest_parameters.json')
+
+
+
+ if os.path.exists(dst_path):
+ i = 1
+ while True:
+ modified_output_file = os.path.join(self.output_dir, f"random_forest_{i}.tif")
+ if not os.path.exists(modified_output_file):
+ dst_path = modified_output_file
+ break
+ i += 1
+
+ if os.path.exists(params_file):
+ i = 1
+ while True:
+ modified_output_file_params = os.path.join(self.output_dir, f"random_forest_parameters_{i}.json")
+ if not os.path.exists(modified_output_file_params):
+ params_file = modified_output_file_params
+ break
+ i += 1
+
+ with rasterio.open(dst_path, 'w', driver='GTiff',
+ height=height, width=width, count=channels, dtype='float32',
+ crs=crs, transform=transform) as dst_ds:
+ dst_ds.write(predicted_types_numeric, 1)
+
+ with open(params_file, 'w') as f:
+ json.dump(params, f, indent=4)
+ feedback.pushInfo(f"Parameters saved to {params_file}")
+
+ parameters['OUTPUT_RASTER']=dst_path
+
+ return {'OUTPUT_RASTER':dst_path}
+
+ def process_options(self,parameters, context, feedback):
+ self.iPatch = 0
+
+ self.feature_dir = ""
+
+ feedback.pushInfo(
+ f'PARAMETERS :\n{parameters}')
+
+ feedback.pushInfo(
+ f'CONTEXT :\n{context}')
+
+ feedback.pushInfo(
+ f'FEEDBACK :\n{feedback}')
+
+ rlayer = self.parameterAsRasterLayer(
+ parameters, self.INPUT, context)
+
+ if rlayer is None:
+ raise QgsProcessingException(
+ self.invalidRasterError(parameters, self.INPUT))
+
+ self.selected_bands = self.parameterAsInts(
+ parameters, self.BANDS, context)
+
+ if len(self.selected_bands) == 0:
+ self.selected_bands = list(range(1, rlayer.bandCount()+1))
+
+ if max(self.selected_bands) > rlayer.bandCount():
+ raise QgsProcessingException(
+ self.tr("The chosen bands exceed the largest band number!")
+ )
+
+ self.template = self.parameterAsFile(
+ parameters, self.TEMPLATE, context)
+
+ self.template_test = self.parameterAsFile(
+ parameters, self.TEMPLATE_TEST, context)
+
+ self.colonne_rf = self.parameterAsString(
+ parameters, self.COLONNE_RF, context)
+
+ res = self.parameterAsDouble(
+ parameters, self.RESOLUTION, context)
+ crs = self.parameterAsCrs(
+ parameters, self.CRS, context)
+ extent = self.parameterAsExtent(
+ parameters, self.EXTENT, context)
+ self.output_dir = self.parameterAsString(
+ parameters, self.OUTPUT, context)
+
+ rlayer_data_provider = rlayer.dataProvider()
+
+ # handle crs
+ if crs is None or not crs.isValid():
+ crs = rlayer.crs()
+ feedback.pushInfo(f'crs : {crs}')
+
+ """
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}') # 0 for meters, 6 for degrees, 9 for unknown
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ if crs.mapUnits() == Qgis.DistanceUnit.Degrees:
+ crs = self.estimate_utm_crs(rlayer.extent())
+
+ # target crs should use meters as units
+ if crs.mapUnits() != Qgis.DistanceUnit.Meters:
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}')
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ raise QgsProcessingException(
+ self.tr("Only support CRS with the units as meters")
+ )
+ """
+
+ # 0 for meters, 6 for degrees, 9 for unknown
+ UNIT_METERS = 0
+ UNIT_DEGREES = 6
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ layer_units = 'degrees'
+ else:
+ layer_units = 'meters'
+ # if res is not provided, get res info from rlayer
+ if np.isnan(res) or res == 0:
+ res = rlayer.rasterUnitsPerPixelX() # rasterUnitsPerPixelY() is negative
+ target_units = layer_units
+ else:
+ # when given res in meters by users, convert crs to utm if the original crs unit is degree
+ if crs.mapUnits() != UNIT_METERS: # Qgis.DistanceUnit.Meters:
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ # estimate utm crs based on layer extent
+ crs = self.estimate_utm_crs(rlayer.extent())
+ else:
+ raise QgsProcessingException(
+ f"Resampling of image with the CRS of {crs.authid()} in meters is not supported.")
+ target_units = 'meters'
+ # else:
+ # res = (rlayer_extent.xMaximum() -
+ # rlayer_extent.xMinimum()) / rlayer.width()
+ self.res = res
+
+ # handle extent
+ if extent.isNull():
+ extent = rlayer.extent() # QgsProcessingUtils.combineLayerExtents(layers, crs, context)
+ extent_crs = rlayer.crs()
+ else:
+ if extent.isEmpty():
+ raise QgsProcessingException(
+ self.tr("The extent for processing can not be empty!"))
+ extent_crs = self.parameterAsExtentCrs(
+ parameters, self.EXTENT, context)
+ # if extent crs != target crs, convert it to target crs
+ if extent_crs != crs:
+ transform = QgsCoordinateTransform(
+ extent_crs, crs, context.transformContext())
+ # extent = transform.transformBoundingBox(extent)
+ # to ensure coverage of the transformed extent
+ # convert extent to polygon, transform polygon, then get boundingBox of the new polygon
+ extent_polygon = QgsGeometry.fromRect(extent)
+ extent_polygon.transform(transform)
+ extent = extent_polygon.boundingBox()
+ extent_crs = crs
+
+ # check intersects between extent and rlayer_extent
+ if rlayer.crs() != crs:
+ transform = QgsCoordinateTransform(
+ rlayer.crs(), crs, context.transformContext())
+ rlayer_extent = transform.transformBoundingBox(
+ rlayer.extent())
+ else:
+ rlayer_extent = rlayer.extent()
+ if not rlayer_extent.intersects(extent):
+ raise QgsProcessingException(
+ self.tr("The extent for processing is not intersected with the input image!"))
+
+ img_width_in_extent = round(
+ (extent.xMaximum() - extent.xMinimum())/self.res)
+ img_height_in_extent = round(
+ (extent.yMaximum() - extent.yMinimum())/self.res)
+
+ # Send some information to the user
+ feedback.pushInfo(
+ f'Layer path: {rlayer_data_provider.dataSourceUri()}')
+ # feedback.pushInfo(
+ # f'Layer band scale: {rlayer_data_provider.bandScale(self.selected_bands[0])}')
+ feedback.pushInfo(f'Layer name: {rlayer.name()}')
+ if rlayer.crs().authid():
+ feedback.pushInfo(f'Layer CRS: {rlayer.crs().authid()}')
+ else:
+ feedback.pushInfo(
+ f'Layer CRS in WKT format: {rlayer.crs().toWkt()}')
+ feedback.pushInfo(
+ f'Layer pixel size: {rlayer.rasterUnitsPerPixelX()}, {rlayer.rasterUnitsPerPixelY()} {layer_units}')
+
+ feedback.pushInfo(f'Bands selected: {self.selected_bands}')
+
+ if crs.authid():
+ feedback.pushInfo(f'Target CRS: {crs.authid()}')
+ else:
+ feedback.pushInfo(f'Target CRS in WKT format: {crs.toWkt()}')
+ feedback.pushInfo(f'Target resolution: {self.res} {target_units}')
+ feedback.pushInfo(
+ (f'Processing extent: minx:{extent.xMinimum():.6f}, maxx:{extent.xMaximum():.6f},'
+ f'miny:{extent.yMinimum():.6f}, maxy:{extent.yMaximum():.6f}'))
+ feedback.pushInfo(
+ (f'Processing image size: (width {img_width_in_extent}, '
+ f'height {img_height_in_extent})'))
+
+ self.rlayer_path = rlayer.dataProvider().dataSourceUri()
+
+ feedback.pushInfo(f'Selected bands: {self.selected_bands}')
+
+ ## passing parameters to self once everything has been processed
+ self.extent = extent
+ self.rlayer = rlayer
+ self.crs = crs
+
+
+ # used to handle any thread-sensitive cleanup which is required by the algorithm.
+ def postProcessAlgorithm(self, context, feedback) -> Dict[str, Any]:
+ return {}
+
+
+ def tr(self, string):
+ """
+ Returns a translatable string with the self.tr() function.
+ """
+ return QCoreApplication.translate('Processing', string)
+
+ def createInstance(self):
+ return RFAlgorithm()
+
+ def name(self):
+ """
+ Returns the algorithm name, used for identifying the algorithm. This
+ string should be fixed for the algorithm, and must not be localised.
+ The name should be unique within each provider. Names should contain
+ lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return 'Random_forest'
+
+ def displayName(self):
+ """
+ Returns the translated algorithm name, which should be used for any
+ user-visible display of the algorithm name.
+ """
+ return self.tr('Random_forest')
+
+ def group(self):
+ """
+ Returns the name of the group this algorithm belongs to. This string
+ should be localised.
+ """
+ return self.tr('')
+
+ def groupId(self):
+ """
+ Returns the unique ID of the group this algorithm belongs to. This
+ string should be fixed for the algorithm, and must not be localised.
+ The group id should be unique within each provider. Group id should
+ contain lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return ''
+
+ def shortHelpString(self):
+ """
+ Returns a localised short helper string for the algorithm. This string
+ should provide a basic description about what the algorithm does and the
+ parameters and outputs associated with it..
+ """
+ return self.tr("Compute Random forest using input template")
+
+ def icon(self):
+ return 'E'
+
+
+
+
diff --git a/src/iamap/reduction.py b/src/iamap/reduction.py
new file mode 100644
index 0000000..0d73303
--- /dev/null
+++ b/src/iamap/reduction.py
@@ -0,0 +1,568 @@
+import os
+import numpy as np
+from pathlib import Path
+from typing import Dict, Any
+import joblib
+import pandas as pd
+from sklearn.preprocessing import StandardScaler
+
+import rasterio
+from rasterio import windows
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (Qgis,
+ QgsGeometry,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterEnum,
+ QgsCoordinateTransform,
+ QgsProcessingException,
+ QgsProcessingAlgorithm,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterBand,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterCrs,
+ QgsProcessingParameterDefinition,
+ )
+
+
+from sklearn.decomposition import PCA
+#from umap.umap_ import UMAP
+from sklearn.cluster import KMeans
+
+import json
+
+
+
+class ReductionAlgorithm(QgsProcessingAlgorithm):
+ """
+ """
+
+ INPUT = 'INPUT'
+ BANDS = 'BANDS'
+ EXTENT = 'EXTENT'
+ LOAD = 'LOAD'
+ OUTPUT = 'OUTPUT'
+ RESOLUTION = 'RESOLUTION'
+ CRS = 'CRS'
+ COMPONENTS = 'COMPONENTS'
+ SUBSET = 'SUBSET'
+ METHOD = 'METHOD'
+ SAVE_MODEL = 'SAVE_MODEL'
+ THRESOLD_PCA = 'THRESOLD_PCA'
+
+
+ def initAlgorithm(self, config=None):
+ """
+ Here we define the inputs and output of the algorithm, along
+ with some other properties.
+ """
+ cwd = Path(__file__).parent.absolute()
+
+ self.addParameter(
+ QgsProcessingParameterRasterLayer(
+ name=self.INPUT,
+ description=self.tr(
+ 'Input raster layer or image file path'),
+ defaultValue=os.path.join(cwd,'rasters','test.tif'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBand(
+ name=self.BANDS,
+ description=self.tr(
+ 'Selected Bands (defaults to all bands selected)'),
+ defaultValue=None,
+ parentLayerParameterName=self.INPUT,
+ optional=True,
+ allowMultiple=True,
+ )
+ )
+
+ crs_param = QgsProcessingParameterCrs(
+ name=self.CRS,
+ description=self.tr('Target CRS (default to original CRS)'),
+ optional=True,
+ )
+
+ res_param = QgsProcessingParameterNumber(
+ name=self.RESOLUTION,
+ description=self.tr(
+ 'Target resolution in meters (default to native resolution)'),
+ type=QgsProcessingParameterNumber.Double,
+ optional=True,
+ minValue=0,
+ maxValue=100000
+ )
+
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ name=self.EXTENT,
+ description=self.tr(
+ 'Processing extent (default to the entire image)'),
+ optional=True
+ )
+ )
+ self.method_opt = ['PCA', 'UMAP', 'K-means', '--Empty--']
+ self.addParameter (
+ QgsProcessingParameterEnum(
+ name = self.METHOD,
+ description = self.tr(
+ 'Method for the dimension reduction'),
+ defaultValue = 0,
+ options = self.method_opt,
+
+ )
+ )
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name = self.THRESOLD_PCA,
+ description = self.tr (
+ 'Thresold for displaying contribution of each variables when using PCA (between 0 and 1)'
+ ),
+ type = QgsProcessingParameterNumber.Double,
+ minValue = 0,
+ defaultValue = 0.5,
+ maxValue = 1
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ name=self.COMPONENTS,
+ description=self.tr(
+ 'Number of target components'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue = 3,
+ minValue=1,
+ maxValue=1024
+ )
+ )
+ subset_param = QgsProcessingParameterNumber(
+ name=self.SUBSET,
+ description=self.tr(
+ 'Select a subset of random pixels of the image to fit transform'),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=None,
+ minValue=1,
+ maxValue=10_000,
+ optional=True,
+ )
+
+ save_param = QgsProcessingParameterBoolean(
+ self.SAVE_MODEL,
+ self.tr("Save projection model after fit."),
+ defaultValue=True
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT,
+ self.tr(
+ "Output directory (choose the location that the image features will be saved)"),
+ defaultValue=os.path.join(cwd,'features'),
+ )
+ )
+
+
+ for param in (crs_param, res_param, subset_param, save_param):
+ param.setFlags(
+ param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(param)
+
+
+
+ def processAlgorithm(self, parameters, context, feedback):
+ """
+ Here is where the processing itself takes place.
+ """
+ self.process_options(parameters, context, feedback)
+
+ input_bands = [i_band -1 for i_band in self.selected_bands]
+
+ if self.method == 'PCA':
+ proj = PCA(int(self.ncomponents))
+ save_file = 'pca.pkl'
+ if self.method == 'UMAP':
+ proj = UMAP(int(self.ncomponents))
+ save_file = 'umap.pkl'
+ if self.method == 'K-means':
+ proj = KMeans(int(self.ncomponents))
+ save_file = 'kmeans_transform.pkl'
+
+ out_path = os.path.join(self.output_dir, save_file)
+
+ with rasterio.open(self.rlayer_path) as ds:
+ transform = ds.transform
+ crs = ds.crs
+ win = windows.from_bounds(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ transform=transform
+ )
+ raster = ds.read(window=win)
+ transform = ds.window_transform(win)
+ raster = np.transpose(raster, (1,2,0))
+ raster = raster[:,:,input_bands]
+
+ feedback.pushInfo(f'{raster.shape}')
+ feedback.pushInfo(f'{raster.reshape(-1, raster.shape[0]).shape}')
+
+ if self.subset:
+ feedback.pushInfo(f'Using a random subset of {self.subset} pixels')
+ fit_raster = raster.reshape(-1, raster.shape[-1])
+ nsamples = fit_raster.shape[0]
+
+ # Generate random indices to select subset_size number of samples
+ np.random.seed(42)
+ random_indices = np.random.choice(nsamples, size=self.subset, replace=False)
+ fit_raster = fit_raster[random_indices,:]
+ feedback.pushInfo(f'Starting fit\n')
+ proj.fit(fit_raster)
+
+ if self.method=='PCA' :
+ eig = proj.components_
+
+ absolute_eig = np.abs(eig)
+ contributions = pd.DataFrame(eig.T, columns=[f'PC{i+1}' for i in range(eig.shape[0])])
+ #feedback.pushInfo(f"Contributions des variables aux composantes principales : {contributions}")
+
+
+ # Filtrer les contributions supérieures à 0.8
+ #threshold = 0.01
+ significant_contributions = contributions[contributions > self.thresold_pca].dropna(how='all')
+ feedback.pushInfo(f"Contributions des variables aux composantes principales : {significant_contributions}")
+ if self.save_model:
+ joblib.dump(proj, out_path)
+
+ feedback.pushInfo(f'starting inference\n')
+ proj_img = proj.transform(raster.reshape(-1, raster.shape[-1]))
+
+
+ else:
+ scaler = StandardScaler()
+ #raster_normalized=scaler.fit_transform(raster)
+ feedback.pushInfo(f"Mean raster : {np.mean(raster)}")
+ feedback.pushInfo(f"Standart dev : {np.std(raster)}")
+ raster = (raster-np.mean(raster))/np.std(raster)
+ feedback.pushInfo(f"Mean raster normalized : {np.mean(raster)}")
+ feedback.pushInfo(f"Standart dev normalized : {np.std(raster)}")
+ #feedback.pushInfo(f"Mean raster normalized: {raster_normalized.mean}")
+ proj_img = proj.fit_transform(raster.reshape(-1, raster.shape[-1]))
+
+ if self.method=='PCA' :
+ eig = proj.components_.T
+ feedback.pushInfo(f"Thresold: {self.thresold_pca}")
+
+ cos2 = eig ** 2
+ cos2_filtered = cos2 > self.thresold_pca
+ feedback.pushInfo(f"Contributions des variables aux composantes principales : {cos2_filtered}")
+ num_values_above_threshold = np.sum(cos2_filtered)
+
+ feedback.pushInfo(f"Number of values in cos2 > threshold: {num_values_above_threshold}")
+ #contributions = pd.DataFrame(eig, columns=[f'PC{i+1}' for i in range(eig.T.shape[0])])
+ #feedback.pushInfo(f"Contributions des variables aux composantes principales : {contributions}")
+
+
+ # Filtrer les contributions supérieures à 0.8
+ #threshold = 0.01
+ #significant_contributions = contributions[contributions > self.thresold_pca].dropna(how='all')
+ #feedback.pushInfo(f"Contributions des variables aux composantes principales : {significant_contributions}")
+
+
+ """
+ for i in range(self.ncomponents):
+ component_loadings = absolute_eig[:, i]
+ top_indices = np.argsort(component_loadings)[-top_n:][::-1]
+ #top_indices = np.argsort(component_loadings)
+ feedback.pushInfo(f'eigenvectors significatifs : {top_indices}')
+ """
+ if self.save_model:
+ joblib.dump(proj, out_path)
+ params = proj.get_params()
+
+ proj_img = proj_img.reshape((raster.shape[0], raster.shape[1],-1))
+ height, width, channels = proj_img.shape
+
+
+ dst_path = os.path.join(self.output_dir,'proj')
+ params_file = os.path.join(self.output_dir, 'reduction_parameters.json')
+
+ if os.path.exists(params_file):
+ i = 1
+ while True:
+ modified_output_file_params = os.path.join(self.output_dir, f"reductions_parameters_{i}.json")
+ if not os.path.exists(modified_output_file_params):
+ params_file = modified_output_file_params
+ break
+ i += 1
+
+
+
+
+ if os.path.exists(dst_path):
+ i = 1
+ while True:
+ modified_output_file = os.path.join(self.output_dir, f"proj_{i}.tif")
+ if not os.path.exists(modified_output_file):
+ dst_path = modified_output_file
+ break
+ i += 1
+
+ with rasterio.open(dst_path, 'w', driver='GTiff',
+ height=height, width=width, count=channels, dtype='float32',
+ crs=crs, transform=transform) as dst_ds:
+ dst_ds.write(np.transpose(proj_img, (2, 0, 1)))
+
+ with open(params_file, 'w') as f:
+ json.dump(params, f, indent=4)
+ feedback.pushInfo(f"Parameters saved to {params_file}")
+
+ parameters['OUTPUT_RASTER']=dst_path
+
+ return {'OUTPUT_RASTER':dst_path}
+
+ def process_options(self,parameters, context, feedback):
+ self.iPatch = 0
+
+ self.feature_dir = ""
+
+ feedback.pushInfo(
+ f'PARAMETERS :\n{parameters}')
+
+ feedback.pushInfo(
+ f'CONTEXT :\n{context}')
+
+ feedback.pushInfo(
+ f'FEEDBACK :\n{feedback}')
+
+ rlayer = self.parameterAsRasterLayer(
+ parameters, self.INPUT, context)
+
+ if rlayer is None:
+ raise QgsProcessingException(
+ self.invalidRasterError(parameters, self.INPUT))
+
+ self.selected_bands = self.parameterAsInts(
+ parameters, self.BANDS, context)
+
+ if len(self.selected_bands) == 0:
+ self.selected_bands = list(range(1, rlayer.bandCount()+1))
+
+ if max(self.selected_bands) > rlayer.bandCount():
+ raise QgsProcessingException(
+ self.tr("The chosen bands exceed the largest band number!")
+ )
+
+ self.ncomponents = self.parameterAsInt(
+ parameters, self.COMPONENTS, context)
+ self.subset = self.parameterAsInt(
+ parameters, self.SUBSET, context)
+ method_idx = self.parameterAsEnum(
+ parameters, self.METHOD, context)
+ self.method = self.method_opt[method_idx]
+
+ self.thresold_pca = self.parameterAsDouble(
+ parameters, self.THRESOLD_PCA, context
+ )
+
+ res = self.parameterAsDouble(
+ parameters, self.RESOLUTION, context)
+ crs = self.parameterAsCrs(
+ parameters, self.CRS, context)
+ extent = self.parameterAsExtent(
+ parameters, self.EXTENT, context)
+ self.output_dir = self.parameterAsString(
+ parameters, self.OUTPUT, context)
+ self.save_model = self.parameterAsBoolean(
+ parameters, self.SAVE_MODEL, context)
+
+ rlayer_data_provider = rlayer.dataProvider()
+
+ # handle crs
+ if crs is None or not crs.isValid():
+ crs = rlayer.crs()
+ #feedback.pushInfo(
+ # f'Layer CRS unit is {crs.mapUnits()}') # 0 for meters, 6 for degrees, 9 for unknown
+ #feedback.pushInfo(
+ # f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ #if crs.mapUnits() == Qgis.DistanceUnit.Degrees:
+ # crs = self.estimate_utm_crs(rlayer.extent())
+
+ # target crs should use meters as units
+ #if crs.mapUnits() != Qgis.DistanceUnit.Meters:
+ # feedback.pushInfo(
+ # f'Layer CRS unit is {crs.mapUnits()}')
+ # feedback.pushInfo(
+ # f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ #raise QgsProcessingException(
+ # self.tr("Only support CRS with the units as meters")
+ #)
+
+ # 0 for meters, 6 for degrees, 9 for unknown
+ UNIT_METERS = 0
+ UNIT_DEGREES = 6
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ layer_units = 'degrees'
+ else:
+ layer_units = 'meters'
+ # if res is not provided, get res info from rlayer
+ if np.isnan(res) or res == 0:
+ res = rlayer.rasterUnitsPerPixelX() # rasterUnitsPerPixelY() is negative
+ target_units = layer_units
+ else:
+ # when given res in meters by users, convert crs to utm if the original crs unit is degree
+ if crs.mapUnits() != UNIT_METERS: # Qgis.DistanceUnit.Meters:
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ # estimate utm crs based on layer extent
+ crs = self.estimate_utm_crs(rlayer.extent())
+ else:
+ raise QgsProcessingException(
+ f"Resampling of image with the CRS of {crs.authid()} in meters is not supported.")
+ target_units = 'meters'
+ # else:
+ # res = (rlayer_extent.xMaximum() -
+ # rlayer_extent.xMinimum()) / rlayer.width()
+ self.res = res
+
+ # handle extent
+ if extent.isNull():
+ extent = rlayer.extent() # QgsProcessingUtils.combineLayerExtents(layers, crs, context)
+ extent_crs = rlayer.crs()
+ else:
+ if extent.isEmpty():
+ raise QgsProcessingException(
+ self.tr("The extent for processing can not be empty!"))
+ extent_crs = self.parameterAsExtentCrs(
+ parameters, self.EXTENT, context)
+ # if extent crs != target crs, convert it to target crs
+ if extent_crs != crs:
+ transform = QgsCoordinateTransform(
+ extent_crs, crs, context.transformContext())
+ # extent = transform.transformBoundingBox(extent)
+ # to ensure coverage of the transformed extent
+ # convert extent to polygon, transform polygon, then get boundingBox of the new polygon
+ extent_polygon = QgsGeometry.fromRect(extent)
+ extent_polygon.transform(transform)
+ extent = extent_polygon.boundingBox()
+ extent_crs = crs
+
+ # check intersects between extent and rlayer_extent
+ if rlayer.crs() != crs:
+ transform = QgsCoordinateTransform(
+ rlayer.crs(), crs, context.transformContext())
+ rlayer_extent = transform.transformBoundingBox(
+ rlayer.extent())
+ else:
+ rlayer_extent = rlayer.extent()
+ if not rlayer_extent.intersects(extent):
+ raise QgsProcessingException(
+ self.tr("The extent for processing is not intersected with the input image!"))
+
+ img_width_in_extent = round(
+ (extent.xMaximum() - extent.xMinimum())/self.res)
+ img_height_in_extent = round(
+ (extent.yMaximum() - extent.yMinimum())/self.res)
+
+ # Send some information to the user
+ feedback.pushInfo(
+ f'Layer path: {rlayer_data_provider.dataSourceUri()}')
+ # feedback.pushInfo(
+ # f'Layer band scale: {rlayer_data_provider.bandScale(self.selected_bands[0])}')
+ feedback.pushInfo(f'Layer name: {rlayer.name()}')
+ if rlayer.crs().authid():
+ feedback.pushInfo(f'Layer CRS: {rlayer.crs().authid()}')
+ else:
+ feedback.pushInfo(
+ f'Layer CRS in WKT format: {rlayer.crs().toWkt()}')
+ feedback.pushInfo(
+ f'Layer pixel size: {rlayer.rasterUnitsPerPixelX()}, {rlayer.rasterUnitsPerPixelY()} {layer_units}')
+
+ feedback.pushInfo(f'Bands selected: {self.selected_bands}')
+
+ if crs.authid():
+ feedback.pushInfo(f'Target CRS: {crs.authid()}')
+ else:
+ feedback.pushInfo(f'Target CRS in WKT format: {crs.toWkt()}')
+ feedback.pushInfo(f'Target resolution: {self.res} {target_units}')
+ feedback.pushInfo(
+ (f'Processing extent: minx:{extent.xMinimum():.6f}, maxx:{extent.xMaximum():.6f},'
+ f'miny:{extent.yMinimum():.6f}, maxy:{extent.yMaximum():.6f}'))
+ feedback.pushInfo(
+ (f'Processing image size: (width {img_width_in_extent}, '
+ f'height {img_height_in_extent})'))
+
+ self.rlayer_path = rlayer.dataProvider().dataSourceUri()
+
+ feedback.pushInfo(f'Selected bands: {self.selected_bands}')
+
+ ## passing parameters to self once everything has been processed
+ self.extent = extent
+ self.rlayer = rlayer
+ self.crs = crs
+
+
+ # used to handle any thread-sensitive cleanup which is required by the algorithm.
+ def postProcessAlgorithm(self, context, feedback) -> Dict[str, Any]:
+ return {}
+
+
+ def tr(self, string):
+ """
+ Returns a translatable string with the self.tr() function.
+ """
+ return QCoreApplication.translate('Processing', string)
+
+ def createInstance(self):
+ return ReductionAlgorithm()
+
+ def name(self):
+ """
+ Returns the algorithm name, used for identifying the algorithm. This
+ string should be fixed for the algorithm, and must not be localised.
+ The name should be unique within each provider. Names should contain
+ lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return 'reduction'
+
+ def displayName(self):
+ """
+ Returns the translated algorithm name, which should be used for any
+ user-visible display of the algorithm name.
+ """
+ return self.tr('Dimension Reduction')
+
+ def group(self):
+ """
+ Returns the name of the group this algorithm belongs to. This string
+ should be localised.
+ """
+ return self.tr('')
+
+ def groupId(self):
+ """
+ Returns the unique ID of the group this algorithm belongs to. This
+ string should be fixed for the algorithm, and must not be localised.
+ The group id should be unique within each provider. Group id should
+ contain lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return ''
+
+ def shortHelpString(self):
+ """
+ Returns a localised short helper string for the algorithm. This string
+ should provide a basic description about what the algorithm does and the
+ parameters and outputs associated with it..
+ """
+ return self.tr("Reduce the dimension of deep learning features.")
+
+ def icon(self):
+ return 'E'
+
+
diff --git a/src/iamap/requirements.txt b/src/iamap/requirements.txt
new file mode 100644
index 0000000..dfbde06
--- /dev/null
+++ b/src/iamap/requirements.txt
@@ -0,0 +1,4 @@
+torchgeo == 0.5.2
+geopandas == 0.14.4
+scikit-learn == 1.5.1
+umap-learn == 0.5.6
diff --git a/src/iamap/similarity.py b/src/iamap/similarity.py
new file mode 100644
index 0000000..c490579
--- /dev/null
+++ b/src/iamap/similarity.py
@@ -0,0 +1,473 @@
+import os
+import numpy as np
+from pathlib import Path
+from typing import Dict, Any
+import joblib
+
+import rasterio
+from rasterio import windows
+import geopandas as gpd
+from shapely.geometry import box
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (Qgis,
+ QgsGeometry,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterEnum,
+ QgsCoordinateTransform,
+ QgsProcessingException,
+ QgsProcessingAlgorithm,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterBand,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterCrs,
+ QgsProcessingParameterDefinition,
+ )
+import torch
+import torch.nn as nn
+
+import json
+
+
+class SimilarityAlgorithm(QgsProcessingAlgorithm):
+ """
+ """
+
+ INPUT = 'INPUT'
+ BANDS = 'BANDS'
+ EXTENT = 'EXTENT'
+ LOAD = 'LOAD'
+ OUTPUT = 'OUTPUT'
+ RESOLUTION = 'RESOLUTION'
+ CRS = 'CRS'
+ TEMPLATE = 'TEMPLATE'
+
+
+ def initAlgorithm(self, config=None):
+ """
+ Here we define the inputs and output of the algorithm, along
+ with some other properties.
+ """
+ cwd = Path(__file__).parent.absolute()
+
+ self.addParameter(
+ QgsProcessingParameterRasterLayer(
+ name=self.INPUT,
+ description=self.tr(
+ 'Input raster layer or image file path'),
+ defaultValue=os.path.join(cwd,'rasters','test.tif'),
+ ),
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBand(
+ name=self.BANDS,
+ description=self.tr(
+ 'Selected Bands (defaults to all bands selected)'),
+ defaultValue=None,
+ parentLayerParameterName=self.INPUT,
+ optional=True,
+ allowMultiple=True,
+ )
+ )
+
+ crs_param = QgsProcessingParameterCrs(
+ name=self.CRS,
+ description=self.tr('Target CRS (default to original CRS)'),
+ optional=True,
+ )
+
+ res_param = QgsProcessingParameterNumber(
+ name=self.RESOLUTION,
+ description=self.tr(
+ 'Target resolution in meters (default to native resolution)'),
+ type=QgsProcessingParameterNumber.Double,
+ optional=True,
+ minValue=0,
+ maxValue=100000
+ )
+
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ name=self.EXTENT,
+ description=self.tr(
+ 'Processing extent (default to the entire image)'),
+ optional=True
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterFile(
+ name=self.TEMPLATE,
+ description=self.tr(
+ 'Input shapefile path for cosine similarity'),
+ defaultValue=os.path.join(cwd,'rasters','template.shp'),
+ ),
+ )
+
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT,
+ self.tr(
+ "Output directory (choose the location that the image features will be saved)"),
+ defaultValue=os.path.join(cwd,'features'),
+ )
+ )
+
+
+ for param in (crs_param, res_param):
+ param.setFlags(
+ param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(param)
+
+
+
+ def processAlgorithm(self, parameters, context, feedback):
+ """
+ Here is where the processing itself takes place.
+ """
+ self.process_options(parameters, context, feedback)
+
+ template_path = self.template # Replace with the actual attribute holding the file path
+
+# Create a dictionary to store the template path
+ dictionnary_cosine = {
+ 'cosine_path': template_path
+ }
+
+ gdf = gpd.read_file(self.template)
+
+ gdf = gdf.to_crs(self.crs.toWkt())
+
+ feedback.pushInfo(f'before extent: {len(gdf)}')
+ bounds = box(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ )
+ gdf = gdf[gdf.within(bounds)]
+ feedback.pushInfo(f'after extent: {len(gdf)}')
+
+ if len(gdf) == 0:
+ feedback.pushWarning("No template points within extent !")
+ return False
+
+
+
+
+ input_bands = [i_band -1 for i_band in self.selected_bands]
+
+
+ with rasterio.open(self.rlayer_path) as ds:
+ gdf = gdf.to_crs(ds.crs)
+ pixel_values = []
+
+ transform = ds.transform
+ crs = ds.crs
+ win = windows.from_bounds(
+ self.extent.xMinimum(),
+ self.extent.yMinimum(),
+ self.extent.xMaximum(),
+ self.extent.yMaximum(),
+ transform=transform
+ )
+ raster = ds.read(window=win)
+ transform = ds.window_transform(win)
+ raster = raster[input_bands,:,:]
+
+
+ for index, data in gdf.iterrows():
+ # Get the coordinates of the point in the raster's pixel space
+ x, y = data.geometry.x, data.geometry.y
+
+ # Convert point coordinates to pixel coordinates within the window
+ col, row = ~transform * (x, y) # Convert from map coordinates to pixel coordinates
+ col, row = int(col), int(row)
+ feedback.pushInfo(f'after extent: {row, col}')
+ pixel_values.append(list(raster[:,row, col]))
+
+ raster = np.transpose(raster, (1,2,0))
+
+ feedback.pushInfo(f'{raster.shape}')
+
+ template_npy = np.asarray(pixel_values)
+ template = torch.from_numpy(template_npy).to(torch.float32)
+ template = torch.mean(template, dim=0)
+
+ feat_img = torch.from_numpy(raster)
+ cos = nn.CosineSimilarity(dim=-1, eps=1e-6)
+
+ sim = cos(feat_img,template)
+ feedback.pushInfo(f'{sim.shape}')
+ sim = sim.unsqueeze(-1)
+ sim = sim.numpy()
+ height, width, channels = sim.shape
+
+ dst_path = os.path.join(self.output_dir,'similarity.tif')
+ params_file = os.path.join(self.output_dir,'cosine.json')
+
+ if os.path.exists(dst_path):
+ i = 1
+ while True:
+ modified_output_file = os.path.join(self.output_dir, f"similarity_{i}.tif")
+ if not os.path.exists(modified_output_file):
+ dst_path = modified_output_file
+ break
+ i += 1
+ if os.path.exists(params_file):
+ i = 1
+ while True:
+ modified_output_file_params = os.path.join(self.output_dir, f"cosine_{i}.json")
+ if not os.path.exists(modified_output_file_params):
+ params_file = modified_output_file_params
+ break
+ i += 1
+
+ with rasterio.open(dst_path, 'w', driver='GTiff',
+ height=height, width=width, count=channels, dtype='float32',
+ crs=crs, transform=transform) as dst_ds:
+ dst_ds.write(np.transpose(sim, (2, 0, 1)))
+
+ with open(params_file, 'w') as f:
+ json.dump(dictionnary_cosine, f, indent=4)
+ feedback.pushInfo(f"Parameters saved to {params_file}")
+
+ parameters['OUTPUT_RASTER']=dst_path
+
+ return {'OUTPUT_RASTER':dst_path}
+
+ def process_options(self,parameters, context, feedback):
+ self.iPatch = 0
+
+ self.feature_dir = ""
+
+ feedback.pushInfo(
+ f'PARAMETERS :\n{parameters}')
+
+ feedback.pushInfo(
+ f'CONTEXT :\n{context}')
+
+ feedback.pushInfo(
+ f'FEEDBACK :\n{feedback}')
+
+ rlayer = self.parameterAsRasterLayer(
+ parameters, self.INPUT, context)
+
+ if rlayer is None:
+ raise QgsProcessingException(
+ self.invalidRasterError(parameters, self.INPUT))
+
+ self.selected_bands = self.parameterAsInts(
+ parameters, self.BANDS, context)
+
+ if len(self.selected_bands) == 0:
+ self.selected_bands = list(range(1, rlayer.bandCount()+1))
+
+ if max(self.selected_bands) > rlayer.bandCount():
+ raise QgsProcessingException(
+ self.tr("The chosen bands exceed the largest band number!")
+ )
+
+ self.template = self.parameterAsFile(
+ parameters, self.TEMPLATE, context)
+
+ res = self.parameterAsDouble(
+ parameters, self.RESOLUTION, context)
+ crs = self.parameterAsCrs(
+ parameters, self.CRS, context)
+ extent = self.parameterAsExtent(
+ parameters, self.EXTENT, context)
+ self.output_dir = self.parameterAsString(
+ parameters, self.OUTPUT, context)
+
+ rlayer_data_provider = rlayer.dataProvider()
+
+ # handle crs
+ if crs is None or not crs.isValid():
+ crs = rlayer.crs()
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}') # 0 for meters, 6 for degrees, 9 for unknown
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ if crs.mapUnits() == Qgis.DistanceUnit.Degrees:
+ crs = self.estimate_utm_crs(rlayer.extent())
+
+ # target crs should use meters as units
+ if crs.mapUnits() != Qgis.DistanceUnit.Meters:
+ feedback.pushInfo(
+ f'Layer CRS unit is {crs.mapUnits()}')
+ feedback.pushInfo(
+ f'whether the CRS is a geographic CRS (using lat/lon coordinates) {crs.isGeographic()}')
+ raise QgsProcessingException(
+ self.tr("Only support CRS with the units as meters")
+ )
+
+ # 0 for meters, 6 for degrees, 9 for unknown
+ UNIT_METERS = 0
+ UNIT_DEGREES = 6
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ layer_units = 'degrees'
+ else:
+ layer_units = 'meters'
+ # if res is not provided, get res info from rlayer
+ if np.isnan(res) or res == 0:
+ res = rlayer.rasterUnitsPerPixelX() # rasterUnitsPerPixelY() is negative
+ target_units = layer_units
+ else:
+ # when given res in meters by users, convert crs to utm if the original crs unit is degree
+ if crs.mapUnits() != UNIT_METERS: # Qgis.DistanceUnit.Meters:
+ if rlayer.crs().mapUnits() == UNIT_DEGREES: # Qgis.DistanceUnit.Degrees:
+ # estimate utm crs based on layer extent
+ crs = self.estimate_utm_crs(rlayer.extent())
+ else:
+ raise QgsProcessingException(
+ f"Resampling of image with the CRS of {crs.authid()} in meters is not supported.")
+ target_units = 'meters'
+ # else:
+ # res = (rlayer_extent.xMaximum() -
+ # rlayer_extent.xMinimum()) / rlayer.width()
+ self.res = res
+
+ # handle extent
+ if extent.isNull():
+ extent = rlayer.extent() # QgsProcessingUtils.combineLayerExtents(layers, crs, context)
+ extent_crs = rlayer.crs()
+ else:
+ if extent.isEmpty():
+ raise QgsProcessingException(
+ self.tr("The extent for processing can not be empty!"))
+ extent_crs = self.parameterAsExtentCrs(
+ parameters, self.EXTENT, context)
+ # if extent crs != target crs, convert it to target crs
+ if extent_crs != crs:
+ transform = QgsCoordinateTransform(
+ extent_crs, crs, context.transformContext())
+ # extent = transform.transformBoundingBox(extent)
+ # to ensure coverage of the transformed extent
+ # convert extent to polygon, transform polygon, then get boundingBox of the new polygon
+ extent_polygon = QgsGeometry.fromRect(extent)
+ extent_polygon.transform(transform)
+ extent = extent_polygon.boundingBox()
+ extent_crs = crs
+
+ # check intersects between extent and rlayer_extent
+ if rlayer.crs() != crs:
+ transform = QgsCoordinateTransform(
+ rlayer.crs(), crs, context.transformContext())
+ rlayer_extent = transform.transformBoundingBox(
+ rlayer.extent())
+ else:
+ rlayer_extent = rlayer.extent()
+ if not rlayer_extent.intersects(extent):
+ raise QgsProcessingException(
+ self.tr("The extent for processing is not intersected with the input image!"))
+
+ img_width_in_extent = round(
+ (extent.xMaximum() - extent.xMinimum())/self.res)
+ img_height_in_extent = round(
+ (extent.yMaximum() - extent.yMinimum())/self.res)
+
+ # Send some information to the user
+ feedback.pushInfo(
+ f'Layer path: {rlayer_data_provider.dataSourceUri()}')
+ # feedback.pushInfo(
+ # f'Layer band scale: {rlayer_data_provider.bandScale(self.selected_bands[0])}')
+ feedback.pushInfo(f'Layer name: {rlayer.name()}')
+ if rlayer.crs().authid():
+ feedback.pushInfo(f'Layer CRS: {rlayer.crs().authid()}')
+ else:
+ feedback.pushInfo(
+ f'Layer CRS in WKT format: {rlayer.crs().toWkt()}')
+ feedback.pushInfo(
+ f'Layer pixel size: {rlayer.rasterUnitsPerPixelX()}, {rlayer.rasterUnitsPerPixelY()} {layer_units}')
+
+ feedback.pushInfo(f'Bands selected: {self.selected_bands}')
+
+ if crs.authid():
+ feedback.pushInfo(f'Target CRS: {crs.authid()}')
+ else:
+ feedback.pushInfo(f'Target CRS in WKT format: {crs.toWkt()}')
+ feedback.pushInfo(f'Target resolution: {self.res} {target_units}')
+ feedback.pushInfo(
+ (f'Processing extent: minx:{extent.xMinimum():.6f}, maxx:{extent.xMaximum():.6f},'
+ f'miny:{extent.yMinimum():.6f}, maxy:{extent.yMaximum():.6f}'))
+ feedback.pushInfo(
+ (f'Processing image size: (width {img_width_in_extent}, '
+ f'height {img_height_in_extent})'))
+
+ self.rlayer_path = rlayer.dataProvider().dataSourceUri()
+
+ feedback.pushInfo(f'Selected bands: {self.selected_bands}')
+
+ ## passing parameters to self once everything has been processed
+ self.extent = extent
+ self.rlayer = rlayer
+ self.crs = crs
+
+
+ # used to handle any thread-sensitive cleanup which is required by the algorithm.
+ def postProcessAlgorithm(self, context, feedback) -> Dict[str, Any]:
+ return {}
+
+
+ def tr(self, string):
+ """
+ Returns a translatable string with the self.tr() function.
+ """
+ return QCoreApplication.translate('Processing', string)
+
+ def createInstance(self):
+ return SimilarityAlgorithm()
+
+ def name(self):
+ """
+ Returns the algorithm name, used for identifying the algorithm. This
+ string should be fixed for the algorithm, and must not be localised.
+ The name should be unique within each provider. Names should contain
+ lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return 'similarity'
+
+ def displayName(self):
+ """
+ Returns the translated algorithm name, which should be used for any
+ user-visible display of the algorithm name.
+ """
+ return self.tr('Similarity')
+
+ def group(self):
+ """
+ Returns the name of the group this algorithm belongs to. This string
+ should be localised.
+ """
+ return self.tr('')
+
+ def groupId(self):
+ """
+ Returns the unique ID of the group this algorithm belongs to. This
+ string should be fixed for the algorithm, and must not be localised.
+ The group id should be unique within each provider. Group id should
+ contain lowercase alphanumeric characters only and no spaces or other
+ formatting characters.
+ """
+ return ''
+
+ def shortHelpString(self):
+ """
+ Returns a localised short helper string for the algorithm. This string
+ should provide a basic description about what the algorithm does and the
+ parameters and outputs associated with it..
+ """
+ return self.tr("Compute Cosine similarity between features")
+
+ def icon(self):
+ return 'E'
+
+
+
+
diff --git a/src/iamap/utils/geo.py b/src/iamap/utils/geo.py
new file mode 100644
index 0000000..4d5d6f8
--- /dev/null
+++ b/src/iamap/utils/geo.py
@@ -0,0 +1,118 @@
+from typing import Callable, Union
+import rasterio
+import geopandas as gpd
+import numpy as np
+import warnings
+from rasterio.io import MemoryFile
+from rasterio.merge import merge
+
+def replace_nan_with_zero(array):
+ array[array != array] = 0 # Replace NaN values with zero
+ return array
+
+def custom_method_avg(merged_data, new_data, merged_mask, new_mask, **kwargs):
+ """Returns the average value pixel.
+ cf. https://amanbagrecha.github.io/posts/2022-07-31-merge-rasters-the-modern-way-using-python/index.html
+ """
+ mask = np.empty_like(merged_mask, dtype="bool")
+ np.logical_or(merged_mask, new_mask, out=mask)
+ np.logical_not(mask, out=mask)
+ np.nanmean([merged_data, new_data], axis=0, out=merged_data, where=mask)
+ np.logical_not(new_mask, out=mask)
+ np.logical_and(merged_mask, mask, out=mask)
+ np.copyto(merged_data, new_data, where=mask, casting="unsafe")
+
+def merge_tiles(
+ tiles:list,
+ dst_path,
+ dtype:str = 'float32',
+ nodata=None,
+ #method:str | Callable ='first',
+ method: Union[str, Callable] = 'first',
+ ):
+ """
+ cf. https://amanbagrecha.github.io/posts/2022-07-31-merge-rasters-the-modern-way-using-python/index.html
+ """
+
+ file_handler = [rasterio.open(ds) for ds in tiles]
+ extents = [ds.bounds for ds in file_handler]
+ # Extract individual bounds
+ lefts, bottoms, rights, tops = zip(*extents)
+ union_extent = (
+ min(lefts), # Left
+ min(bottoms), # Bottom
+ max(rights), # Right
+ max(tops) # Top
+ )
+
+ if method == 'average':
+ method = custom_method_avg
+
+ # memfile = MemoryFile()
+ merge(datasets=file_handler, # list of dataset objects opened in 'r' mode
+ bounds=union_extent, # tuple
+ nodata=nodata, # float
+ dtype=dtype, # dtype
+ # resampling=Resampling.nearest,
+ method=method, # strategy to combine overlapping rasters
+ # dst_path=memfile.name, # str or PathLike to save raster
+ dst_path=dst_path,
+ # dst_kwds={'blockysize':512, 'blockxsize':512} # Dictionary
+ )
+
+def get_mean_sd_by_band(path, force_compute=True, ignore_zeros=True):
+ '''
+ Reads metadata or computes mean and sd of each band of a geotiff.
+ If the metadata is not available, mean and standard deviation can be computed via numpy.
+
+ Parameters
+ ----------
+ path : str
+ path to a geotiff file
+ ignore_zeros : boolean
+ ignore zeros when computing mean and sd via numpy
+
+ Returns
+ -------
+ means : list
+ list of mean values per band
+ sds : list
+ list of standard deviation values per band
+ '''
+
+ src = rasterio.open(path)
+ means = []
+ sds = []
+ if force_compute:
+ for band in range(1, src.count+1):
+ arr = src.read(band)
+ arr = replace_nan_with_zero(arr)
+ if ignore_zeros:
+ mean = np.ma.masked_equal(arr, 0).mean()
+ sd = np.ma.masked_equal(arr, 0).std()
+ else:
+ mean = np.mean(arr)
+ sd = np.std(arr)
+ means.append(float(mean))
+ sds.append(float(sd))
+
+ else:
+ for band in range(1, src.count+1):
+ try:
+ tags = src.tags(band)
+ if 'STATISTICS_MEAN' in tags and 'STATISTICS_STDDEV' in tags:
+ mean = float(tags['STATISTICS_MEAN'])
+ sd = float(tags['STATISTICS_STDDEV'])
+ means.append(mean)
+ sds.append(sd)
+ else:
+ raise KeyError("Statistics metadata not found.")
+
+ except KeyError:
+ warnings.warn("Statistics metadata not found and computation not enabled.", UserWarning)
+ except Exception as e:
+ print(f"Error processing band {band}: {e}")
+
+ src.close()
+ return means, sds
+
diff --git a/src/iamap/utils/torchgeo.py b/src/iamap/utils/torchgeo.py
new file mode 100644
index 0000000..122fdfc
--- /dev/null
+++ b/src/iamap/utils/torchgeo.py
@@ -0,0 +1,44 @@
+from torchgeo.samplers import GridGeoSampler
+from collections.abc import Iterator
+
+from torchgeo.datasets import BoundingBox
+from torchgeo.samplers.utils import tile_to_chips
+
+class NoBordersGridGeoSampler(GridGeoSampler):
+
+ def __iter__(self) -> Iterator[BoundingBox]:
+ """
+ Modification of original Torchgeo sampler to avoid overlapping borders of a dataset.
+ Return the index of a dataset.
+
+ Returns:
+ (minx, maxx, miny, maxy, mint, maxt) coordinates to index a dataset
+ """
+ # For each tile...
+ for hit in self.hits:
+ bounds = BoundingBox(*hit.bounds)
+ rows, cols = tile_to_chips(bounds, self.size, self.stride)
+ mint = bounds.mint
+ maxt = bounds.maxt
+
+ # For each row...
+ for i in range(rows):
+ miny = bounds.miny + i * self.stride[0]
+ maxy = miny + self.size[0]
+ if maxy > bounds.maxy:
+ maxy = bounds.maxy
+ miny = bounds.maxy - self.size[0]
+ if miny < bounds.miny:
+ miny = bounds.miny
+
+ # For each column...
+ for j in range(cols):
+ minx = bounds.minx + j * self.stride[1]
+ maxx = minx + self.size[1]
+ if maxx > bounds.maxx:
+ maxx = bounds.maxx
+ minx = bounds.maxx - self.size[1]
+ if minx < bounds.minx:
+ minx = bounds.minx
+
+ yield BoundingBox(minx, maxx, miny, maxy, mint, maxt)
--
GitLab