PoC

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()