path: root/gis-observation-map

#!/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()