Ho trovato in rete un dataset di foto da drone con gia' classificate le fratture. Si tratta di dati con licenza permissiva
I dati derivano da questo articolo
@article{article,
author = {Ovaskainen, Nikolas and Nordbäck, Nicklas and Skyttä, Pietari and Engström, Jon},
year = {2022},
month = {01},
pages = {104528},
title = {A new subsampling methodology to optimize the characterization of two-dimensional bedrock fracture networks},
volume = {155},
journal = {Journal of Structural Geology},
doi = {10.1016/j.jsg.2022.104528}
}
e sono disponibili come ortofoto RGB da drone a questo link https://zenodo.org/records/4719627 (16 Gb) mentre le fratture sono tematizzate in un diversi file gpkg a questo link https://www.kaggle.com/datasets/nasurz/getaberget-fracture-trace-dataset/data (4Mb)
il problema e' che in questo formato le immagini non sono adatte al processamento per la rete neurale. L'idea e' quindi di dividere le ortofoto in tasselli da 15 metri e creare una immagine maschera con le sole fratturazioni. Tale compito e' stato eseguito con uno script Python in QGis
Nello script si devono impostare le coordinate geografiche della finestra ed il folder dove saranno salvate le varie tiles. In un primo passaggio si disattiva dal progetto QGis il tema lineare delle fratture per avere solo i tasselli delle ortofoto e dopo si disattiva il tematismo raster per creare le maschere
------------------------------------------------------------------------------------------
from qgis.core import (
QgsProject, QgsApplication, QgsLayout, QgsPrintLayout, QgsLayoutItemMap,
QgsUnitTypes, QgsLayoutSize, QgsLayoutPoint, QgsLayoutExporter, QgsRectangle
)
import os
# Path to your QGIS project
project_path = "C:\\Users\\l.innocenti\\Desktop\\geta.qgz"
# Output folder
output_dir = "C:\\Users\\l.innocenti\\Desktop\\geta5\\mask\\"
os.makedirs(output_dir, exist_ok=True)
# Load project
project = QgsProject.instance()
project.read(project_path)
from qgis.core import QgsRectangle
# Compute project extent from all layers
extent = None
for layer in project.mapLayers().values():
if not layer.isValid():
continue
layer_extent = layer.extent()
if extent is None:
extent = QgsRectangle(layer_extent)
else:
extent.combineExtentWith(layer_extent)
if extent is None:
raise ValueError("No valid layers found in the project.")
xmin, ymin = 108914, 6720011
xmax, ymax = 108969, 6720066
# Get the full extent in project CRS
#extent = project.boundingBox()
#xmin, ymin, xmax, ymax = extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum()
# Tile size in meters
tile_size = 15
# Iterate over the grid
x = xmin
tile_id = 0
while x < xmax:
y = ymin
while y < ymax:
# Tile extent
tile_extent = QgsRectangle(x, y, x + tile_size, y + tile_size)
# Create a print layout
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(f"tile_{tile_id}")
pc = layout.pageCollection().pages()[0]
# Map item
map_item = QgsLayoutItemMap(layout)
map_item.setRect(20, 20, 200, 200)
map_item.setExtent(tile_extent)
map_item.setBackgroundColor(QColor(255, 255, 255, 0))
map_item.attemptResize(QgsLayoutSize(200, 200, QgsUnitTypes.LayoutMillimeters))
visible_layers = [lyr for lyr in iface.mapCanvas().layers()]
map_item.setLayers(visible_layers)
layout.addLayoutItem(map_item)
layout.addLayoutItem(map_item)
map_item.attemptMove(QgsLayoutPoint(0, 0, QgsUnitTypes.LayoutMillimeters))
# Export PNG
exporter = QgsLayoutExporter(layout)
out_path = os.path.join(output_dir, f"tile_{tile_id}.png")
exporter.exportToImage(out_path, QgsLayoutExporter.ImageExportSettings())
print(f"Saved {out_path}")
tile_id += 1
y += tile_size
x += tile_size
------------------------------------------------------------------------------------------
alla fine si hanno coppie di questo tipo (le immagini sono state rifilate con imagemagick per togliere lo spazio bianco inserito da QGis)
.png)
per "aggiustare" il background con colore nero e il target come grigio come 127
import cv2
import numpy as np
import sys
# --- Check for input filename ---
if len(sys.argv) < 2:
print("Usage: python script.py <image_filename>")
sys.exit(1)
filename = sys.argv[1]
# Read image
img = cv2.imread(filename)
if img is None:
raise FileNotFoundError(f"Cannot read file: {filename}")
# Define white color in BGR
white = np.array([255, 255, 255], dtype=np.uint8)
# Create mask for white pixels
mask_white = np.all(img == white, axis=-1)
# Create output image
output = np.full_like(img, (128, 128, 128)) # Set all to gray
output[mask_white] = (0, 0, 0) # Set white pixels to black
# Save with same filename (overwrite original)
cv2.imwrite(filename, output)
print(f"Image processed and saved: {filename}")
