diff options
| -rwxr-xr-x | gis-observation-map | 391 |
1 files changed, 391 insertions, 0 deletions
diff --git a/gis-observation-map b/gis-observation-map new file mode 100755 index 0000000..7e97fb9 --- /dev/null +++ b/gis-observation-map @@ -0,0 +1,391 @@ +#!/usr/bin/python3 + +#---------------------------------------------------------------------- +# Create a QGIS project with observations from Artdatabanken +# Copyright © 2023 Guilhem Moulin <guilhem@fripost.org> +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <https://www.gnu.org/licenses/>. +#---------------------------------------------------------------------- + +import argparse +import sys +import json +from pathlib import Path +from osgeo import gdal, ogr, osr +gdal.UseExceptions() + +from qgis.core import ( + QgsApplication, + QgsCoordinateReferenceSystem, + QgsLayerDefinition, + QgsLayerTreeModel, + QgsMapLayer, + QgsProject, + QgsUnitTypes, + QgsVectorLayer, + QgsWkbTypes +) +from qgis.gui import QgsLayerTreeView, QgsMapCanvas +from PyQt5.QtGui import QColor +from PyQt5.QtCore import Qt + +target_srs = osr.SpatialReference() +target_srs.ImportFromEPSG(3006) +target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) +gpkg_drv = ogr.GetDriverByName('GPKG') +geojson_drv = ogr.GetDriverByName('GeoJSON') + +parser = argparse.ArgumentParser( + description='Create a QGIS project with observations from Artdatabanken.', + usage='''%(prog)s --project-home=DIR --project-name=NAME {--geometry=FILE|--point=X,Y} ... + %(prog)s {--geometry=FILE|--point=X,Y} ...''' +) +parser.add_argument('--margin', type=int, help='Margin (in meters) from geometry envelopes') +parser.add_argument('--topo-basedir', default='.', help='Base directory for "Topografi 10", "Topografi 50", etc.') +parser.add_argument('--project-name', help='Project filename (and title)') +parser.add_argument('--project-home', help='Project home directory') +parser.add_argument('--geometry', nargs='*', default=[], help='Geometry filename(s)') +parser.add_argument('--point', nargs='*', default=[], help=f'Coordinates of interest (X,Y in {target_srs.GetName()})') + +args = parser.parse_args() + +if args.project_home is not None and args.project_name is not None: + QgsApplication.setPrefixPath('/usr/bin/qgis', True) + qgs = QgsApplication([], False) + qgs.initQgis() + + projectHome = Path(args.project_home) + projectHome = projectHome.expanduser() + projectHome.mkdir(mode=0o755, exist_ok=True) + + projectPath = projectHome.joinpath(args.project_name) + if not projectPath.suffix in ['.qgs', '.qgz']: + projectPath = projectPath.with_suffix('.qgz') + if projectPath.exists(): + raise Exception(f'{projectPath}: file exists') + + projectInstance = QgsProject.instance() + projectInstance.setFileName(projectPath.as_posix()) + projectInstance.setDistanceUnits(QgsUnitTypes.DistanceMeters) + projectInstance.setAreaUnits(QgsUnitTypes.AreaSquareMeters) + projectInstance.setTitle(Path(args.project_name).stem) + layerTreeRoot = projectInstance.layerTreeRoot() + + qgis_crs = QgsCoordinateReferenceSystem(target_srs.ExportToWkt()) + if qgis_crs.isValid(): + projectInstance.setCrs(qgis_crs) + else: + raise Exception('Invalid CRS') + +elif args.project_home is None or args.project_name is None: + parser.error() + +if len(args.geometry) == 0 and len(args.point) == 0: + parser.error() +if len(args.point) > 0 and (args.margin is None or args.margin <= 0): + parser.error("--point requires positive --margin") + +if projectInstance is not None: + if len(args.geometry) > 1: + sourceGroup = layerTreeRoot.addGroup('Source geometry') + sourceGroup.setItemVisibilityChecked(False) + + for arg_idx, path in enumerate(args.geometry): + path = Path(path) + path = path.expanduser() + layer = QgsVectorLayer(path.as_posix(), path.stem, 'ogr') + if not layer.isValid(): + raise Exception(f'ERROR: {path}: failed to load in QGIS') + + layer.setReadOnly(True) + layer.setFlags(QgsMapLayer.Identifiable | QgsMapLayer.Searchable) + if layer.geometryType() == QgsWkbTypes.PointGeometry: + symbol = layer.renderer().symbol().symbolLayers()[0] + symbol.setSize(2.0) + symbol.setFillColor(QColor('#db1e2a')) + symbol.setStrokeColor(QColor('#801119')) + symbol.setStrokeStyle(Qt.SolidLine) + symbol.setStrokeWidth(0) + elif layer.geometryType() == QgsWkbTypes.LineGeometry: + symbol = layer.renderer().symbol().symbolLayers()[0] + symbol.setColor(QColor('black')) + symbol.setPenStyle(Qt.DashLine) + symbol.setWidth(0.39) + elif layer.geometryType() == QgsWkbTypes.PolygonGeometry: + symbol = layer.renderer().symbol().symbolLayers()[0] + symbol.setBrushStyle(Qt.NoBrush) + symbol.setStrokeColor(QColor('black')) + symbol.setStrokeStyle(Qt.DashLine) + symbol.setStrokeWidth(0.39) + + if len(args.geometry) == 1: + projectInstance.addMapLayer(layer) + layerTreeRoot.findLayer(layer.id()).setItemVisibilityChecked(False) + else: + projectInstance.addMapLayer(layer, False) + sourceGroup.addLayer(layer) + if arg_idx == 0: + #layerTreeModel = QgsLayerTreeModel(layerTreeRoot) + #layerTreeView = QgsLayerTreeView() + #layerTreeView.setModel(layerTreeModel) + #layerTreeView.setCurrentLayer(layer) + canvas = QgsMapCanvas() + canvas.setExtent( layer.extent() ) + canvas.refresh() + +def gpkg_intersects(path, geometries): + ds = gpkg_drv.Open(path, update=0) + for i in range(ds.GetLayerCount()): + layer = ds.GetLayerByIndex(i) + source_srs = layer.GetSpatialRef() + transform_srs = osr.CoordinateTransformation(source_srs, target_srs) + + feature = layer.GetNextFeature() + while feature is not None: + geometry = feature.GetGeometryRef() + geometry.Transform(transform_srs) + for source_geometry in geometries: + if geometry.Intersects(source_geometry): + ds = None + return True + feature = layer.GetNextFeature() + return False + +def idx_intersects(idx_path, geometries): + ds = geojson_drv.Open(idx_path, update=0) + if ds.GetLayerCount() != 1: + print(f'WARN: {idx_path}: has {ds.GetLayerCount()} != 1 layers', file=sys.stderr) + + layer = ds.GetLayer() + locFieldIdx = layer.FindFieldIndex("location", True) + locations = [] + + feature = layer.GetNextFeature() + while feature is not None: + geometry = feature.GetGeometryRef() + geometry.Transform(transform_srs) + for source_geometry in geometries: + if geometry.Intersects(source_geometry): + loc = feature.GetField(locFieldIdx) + locations.append(loc) + break + feature = layer.GetNextFeature() + ds = None + return locations + +def find_qlr(basedir): + if not basedir.is_dir(): + return None + qlr = sorted(basedir.glob('*.qlr')) + if len(qlr) == 0 or not qlr[0].is_file(): + print(f'WARN: {basedir}: does not contain any layer definition file, skipping', file=sys.stderr) + return None + return qlr[0] + +geometries = [] +for pt in args.point: + x, y = pt.split(',', 1) + x = int(x) + y = int(y) + + minX = x - args.margin + maxX = x + args.margin + minY = y - args.margin + maxY = y + args.margin + + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint_2D(minX, maxY) + ring.AddPoint_2D(maxX, maxY) + ring.AddPoint_2D(maxX, minY) + ring.AddPoint_2D(minX, minY) + ring.AddPoint_2D(minX, maxY) + + bbox = ogr.Geometry(ogr.wkbPolygon) + bbox.AddGeometry(ring) + + geometries.append(bbox) + +for path in args.geometry: + path = Path(path) + path = path.expanduser() + ds = ogr.Open(path.as_posix(), update=0) + for i in range(ds.GetLayerCount()): + layer = ds.GetLayerByIndex(i) + source_srs = layer.GetSpatialRef() + transform_srs = osr.CoordinateTransformation(source_srs, target_srs) + + feature = layer.GetNextFeature() + while feature is not None: + geometry = feature.GetGeometryRef() + geometry.Transform(transform_srs) + + if args.margin is None: + geometry = geometry.Clone() + geometry.FlattenTo2D() + else: + minX, maxX, minY, maxY = geometry.GetEnvelope() + minX -= args.margin + maxX += args.margin + minY -= args.margin + maxY += args.margin + + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint_2D(minX, maxY) + ring.AddPoint_2D(maxX, maxY) + ring.AddPoint_2D(maxX, minY) + ring.AddPoint_2D(minX, minY) + ring.AddPoint_2D(minX, maxY) + + bbox = ogr.Geometry(ogr.wkbPolygon) + bbox.AddGeometry(ring) + geometry = bbox + + geometries.append(geometry) + feature = layer.GetNextFeature() + +def geometricFilterPolygon(polygon): + polygon2 = [] + if polygon.GetGeometryType() != ogr.wkbPolygon: + raise Exception('Not a polygon') + + for i in range(polygon.GetGeometryCount()): + ring = polygon.GetGeometryRef(i) + n = ring.GetPointCount() + if ring.GetGeometryType() != ogr.wkbLineString: + raise Exception('Not a linestring') + if n > 0 and ring.GetPoint(0) != ring.GetPoint(n-1): + raise Exception('Not a linear ring') + ring2 = [] + for i in range(n): + lat, lon = ring.GetPoint_2D(i) + ring2.append([lat, lon]) + polygon2.append(ring2) + + return { + 'type': 'polygon', + 'coordinates': polygon2 + } + +# https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md#geographics-filter +def geometricFilter(geometries): + wgs84 = osr.SpatialReference() + wgs84.ImportFromEPSG(4326) # WGS84 + towgs84 = osr.CoordinateTransformation(target_srs, wgs84) + + myGeometries = [] + for geom in geometries: + geom = geom.Clone() + geom.Transform(towgs84) + geometryType = geom.GetGeometryType() + if geometryType == ogr.wkbPoint or geometryType == ogr.wkbMultiPoint: + parser.error('Point geometry requires positive --margin') + elif geometryType == ogr.wkbLineString or geometryType == ogr.wkbMultiLineString: + parser.error('Line geometry requires positive --margin') + elif geometryType == ogr.wkbPolygon: + gfp = geometricFilterPolygon(geom) + myGeometries.append(gfp) + elif geometryType == ogr.wkbMultiPolygon: + for i in range(geom.GetGeometryCount()): + geom2 = geom.GetGeometryRef(i) + gfp = geometricFilterPolygon(geom2) + myGeometries.append(gfp) + else: + raise Exception('Unsuported geometry (try --margin=0)') + + return myGeometries + +geograficsFilter = { + 'geometries': geometricFilter(geometries), + #'maxAccuracy': 5000, + 'considerObservationAccuracy': True, + 'considerDisturbanceRadius': True +} +print(json.dumps(geograficsFilter, indent=4)) +exit() + +topo_maps = ['Topografi 10', 'Topografi 50', 'Topografi 100'] +if projectInstance is not None: + currentMaxScale = 1 + basedir = Path(args.topo_basedir) + basedir = basedir.expanduser() + + for topo in topo_maps: + qlr_paths = [] + d = basedir.joinpath(topo) + if not d.is_dir(): + continue + + idx = d.joinpath('index.geojson') + if idx.is_file(): + tiles = idx_intersects(idx.as_posix(), geometries) + for tile in tiles: + t = d.joinpath(tile) + qlr_path = find_qlr(t) + if qlr_path is None: + continue + qlr_paths.append(qlr_path) + + else: + for t in sorted(d.glob('*')): + qlr_path = find_qlr(t) + if qlr_path is None: + continue + + gpkg_name = 'mark.gpkg' + gpkg_path = t.joinpath(gpkg_name) + if not gpkg_path.is_file(): + print(f'WARN: {t}: lacks {gpkg_name}, skipping', file=sys.stderr) + continue + + if not gpkg_intersects(gpkg_path.as_posix(), geometries): + continue + + qlr_paths.append(qlr_path) + + if len(qlr_paths) > 0: + rootGroup = layerTreeRoot.addGroup(topo) + for qlr_path in qlr_paths: + QgsLayerDefinition.loadLayerDefinition(qlr_path.as_posix(), projectInstance, rootGroup) + lyrSubTree = rootGroup.children() + if len(lyrSubTree) == len(qlr_paths): + for i, lt in enumerate(lyrSubTree): + name = qlr_paths[i].parent.stem + lt.setName(name) + for lt in lyrSubTree: + lt.setExpanded(False) + lt.setItemVisibilityChecked(True) + for lyr in lt.findLayers(): + lyr.setItemVisibilityChecked(True) + lyr = lyr.layer() + lyr.setReadOnly(True) + lyr.setFlags(QgsMapLayer.Removable) + + if topo == 'Topografi 10': + lyr.setScaleBasedVisibility(True) + lyr.setMaximumScale(currentMaxScale) + currentMaxScale = 50000 + lyr.setMinimumScale(currentMaxScale) + elif topo == 'Topografi 50': + lyr.setScaleBasedVisibility(True) + lyr.setMaximumScale(currentMaxScale) + currentMaxScale = 100000 + lyr.setMinimumScale(currentMaxScale) + elif topo == 'Topografi 100': + lyr.setScaleBasedVisibility(True) + lyr.setMaximumScale(currentMaxScale) + currentMaxScale = 10000000 + lyr.setMinimumScale(currentMaxScale) + rootGroup.setExpanded(False) + + projectInstance.write() |