From 4b0aff7b2aded657b4cf107bfe543e44bbd8d747 Mon Sep 17 00:00:00 2001 From: Joaquin Gottlebe Date: Wed, 30 Jul 2025 17:17:44 +0200 Subject: [PATCH] feat: added some missing scripts. --- python-scripts/borders.py | 90 ++++++++++++++ python-scripts/center_of_points.py | 38 ++++++ python-scripts/map.py | 185 +++++++++++++++++++++++++++++ 3 files changed, 313 insertions(+) create mode 100644 python-scripts/borders.py create mode 100644 python-scripts/center_of_points.py create mode 100644 python-scripts/map.py diff --git a/python-scripts/borders.py b/python-scripts/borders.py new file mode 100644 index 0000000..41bf228 --- /dev/null +++ b/python-scripts/borders.py @@ -0,0 +1,90 @@ +from qgis.PyQt.QtCore import QVariant +from qgis.core import (QgsFeature, QgsField, QgsGeometry, QgsVectorLayer, QgsProject) +import requests +import geopandas as gpd + +# Parameters (Set these manually in the script) +country_code = 'DEU' # Country code (ISO-3) +boundary_type = 'ADM0' # Boundary type (ADM0, ADM1, ADM2, etc.) +release_type = 'gbOpen' # Release type ('gbOpen', 'gbHumanitarian', 'gbAuthorative') + +# Function to fetch geoboundary data +def fetch_geoboundary(release_type, country_code, boundary_type): + api_url = f"https://www.geoboundaries.org/api/current/{release_type}/{country_code}/{boundary_type}/" + try: + response = requests.get(api_url) + if response.status_code != 200: + print(f"Failed to fetch data: HTTP Status Code {response.status_code}") + return None + + data = response.json() + results = [] + + if not isinstance(data, list): + data = [data] + + for country_data in data: + if 'gjDownloadURL' not in country_data: + print(f"'gjDownloadURL' not found for {country_data.get('boundaryISO', 'unknown country')}") + continue + + geojson_url = country_data['gjDownloadURL'] + gdf = gpd.read_file(geojson_url) + if gdf.empty: + print(f"Loaded GeoDataFrame for {country_code} is empty.") + continue + + metadata = {key: country_data.get(key, '') for key in country_data} + results.append({'gdf': gdf, 'metadata': metadata}) + + return results + + except requests.RequestException as e: + print(f"Request error: {e}") + except Exception as e: + print(f"An unexpected error occurred: {e}") + + return None + +# Fetch geoboundary data +results = fetch_geoboundary(release_type, country_code, boundary_type) + +if not results: + raise Exception('Failed to fetch geoBoundary data') + +# Create a new memory layer to store the boundary data +layer_name = f"{country_code}_{boundary_type}_{release_type}_geoBoundaries" +vector_layer = QgsVectorLayer("Polygon?crs=epsg:4326", layer_name, "memory") +pr = vector_layer.dataProvider() + +# Define fields for the new layer if needed +pr.addAttributes([QgsField("id", QVariant.Int)]) +vector_layer.updateFields() + +# Loop over the fetched data and add features to the new layer +for result in results: + gdf = result['gdf'] + + for index, row in gdf.iterrows(): + # Check if the row has a valid geometry + if not hasattr(row['geometry'], 'wkt'): + print(f"Unexpected geometry type for row {index}") + continue + + # Create a new feature and set the geometry from the WKT + feat = QgsFeature() + try: + feat.setGeometry(QgsGeometry.fromWkt(row['geometry'].wkt)) + feat.setAttributes([index]) # Set any attributes if necessary + pr.addFeature(feat) + except Exception as e: + print(f"Error setting geometry from WKT: {e}") + continue + +# Update the layer's extents +vector_layer.updateExtents() + +# Add the new layer to the QGIS project +QgsProject.instance().addMapLayer(vector_layer) + +print(f"GeoBoundary layer '{layer_name}' added to the project.") diff --git a/python-scripts/center_of_points.py b/python-scripts/center_of_points.py new file mode 100644 index 0000000..8caad4d --- /dev/null +++ b/python-scripts/center_of_points.py @@ -0,0 +1,38 @@ +from qgis.core import QgsVectorLayer, QgsFeature, QgsGeometry, QgsPointXY, QgsFields, QgsField, QgsWkbTypes +from qgis.PyQt.QtCore import QVariant + +layer = iface.activeLayer() + +sum_x = 0 +sum_y = 0 +count = 0 + +for feature in layer.getFeatures(): + geom = feature.geometry() + if geom.isMultipart(): + points = geom.asMultiPoint() + else: + points = [geom.asPoint()] + for point in points: + sum_x += point.x() + sum_y += point.y() + count += 1 + +mean_x = sum_x / count +mean_y = sum_y / count +centroid_point = QgsGeometry.fromPointXY(QgsPointXY(mean_x, mean_y)) + +crs = layer.crs().authid() +new_layer = QgsVectorLayer(f"Point?crs={crs}", "Centroid Layer", "memory") +new_layer_data_provider = new_layer.dataProvider() +new_layer_data_provider.addAttributes([QgsField("id", QVariant.Int)]) +new_layer.updateFields() +new_feature = QgsFeature() +new_feature.setGeometry(centroid_point) +new_feature.setAttributes([1]) +new_layer_data_provider.addFeature(new_feature) +new_layer.updateExtents() + +QgsProject.instance().addMapLayer(new_layer) + +print(f"New layer 'Centroid Layer' created withcentroid at: ({mean_x}, {mean_y})") diff --git a/python-scripts/map.py b/python-scripts/map.py new file mode 100644 index 0000000..e2bb8ad --- /dev/null +++ b/python-scripts/map.py @@ -0,0 +1,185 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/Applications/QGIS.app/Contents/Resources/python" +LAYER_PATH = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/lab6/stacked/stacked.tif" +STYLE_PATH = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/styles/true_color_2.qml" +EXPORT_PATH = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/lab6/maps/true_color_water.png" +NORTH_ARROW_PATH = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +SECOND_LAYER_PATH = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/lab6/kmeans/kmeans_5.tif" +STYLE_PATH_SECOND = "/Users/huaqo/Developer/courses/2024_Remote_Sensing/styles/overlay_water.qml" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "True Color" +SECOND_LAYER_NAME = "Water" + +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "True Color Landsat and K-means Water class" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +fontfamily = "Arial" +info_font_size = 12 +info_pos_x = legend_pos_x + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +second_raster_layer = QgsRasterLayer(SECOND_LAYER_PATH, SECOND_LAYER_NAME) +if not second_raster_layer.isValid(): + raise Exception(f"{SECOND_LAYER_NAME} layer failed to load!") + + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +QgsProject.instance().addMapLayer(second_raster_layer) +if not second_raster_layer.loadNamedStyle(STYLE_PATH_SECOND): + raise Exception(f"Style {STYLE_PATH_SECOND} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis()