#!/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
import re
import requests
import configparser
import tempfile
from requests.exceptions import HTTPError
from xdg import xdg_config_home
from pathlib import Path
from osgeo import gdal, ogr, osr
gdal.UseExceptions()

from qgis.core import (
    QgsApplication,
    QgsCoordinateReferenceSystem,
    QgsLayerDefinition,
    QgsLayerTreeModel,
    QgsMapLayer,
    QgsPointXY,
    QgsProject,
    QgsRectangle,
    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')

programName = 'gis-observation-map'

config = configparser.ConfigParser()
config.read(Path(xdg_config_home()).joinpath(programName).joinpath('config'))
if 'QGIS' not in config.keys() or config['QGIS'] is None:
    config['QGIS'] = {}

def ePath(v):
    return Path(v).expanduser()

parser = argparse.ArgumentParser(
    description='Create a QGIS project with observations from Artdatabanken.',
    prog = programName,
    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=config['QGIS']['topo-basedir'], type=ePath,
    help='Base directory for "Topografi 10", "Topografi 50", etc.')
parser.add_argument('--project-name',
    help='Project filename (and title)')
parser.add_argument('--project-home', type=ePath,
    help='Project home directory')
parser.add_argument('--geometry', nargs='*', default=[], type=ePath,
    help='Geometry filename(s)')
parser.add_argument('--point', nargs='*', metavar='X,Y', default=[],
    help=f'Coordinates of interest (in {target_srs.GetName()})')
parser.add_argument('--observation-format', default='GPKG',
    help='Format for the observation file (default: %(default)s)')
parser.add_argument('--style', default=config['QGIS']['style'], type=ePath,
    help='QGIS Layer Style File (*.qml) to apply to the observation layer (default: %(default)s)')

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

    args.project_home.mkdir(mode=0o755, exist_ok=True)
    projectPath = args.project_home.joinpath(args.project_name)
    if projectPath.suffix not 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.print_usage()
    exit(1)

if len(args.geometry) == 0 and len(args.point) == 0:
    parser.print_usage()
    exit(1)
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):
        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() )

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()
    source_srs = layer.GetSpatialRef()
    transform_srs = osr.CoordinateTransformation(source_srs, target_srs)
    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:
    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()

if len(args.geometry) == 0:
    r = geometries[0].GetGeometryRef(0)
    pt0X, pt0Y = r.GetPoint_2D(0)
    pt1X, pt1Y = r.GetPoint_2D(2)
    canvas = QgsMapCanvas()
    extent = QgsRectangle(QgsPointXY(pt0X, pt0Y), QgsPointXY(pt1X, pt1Y))
    canvas.setExtent(extent)

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):
            lon, lat = ring.GetPoint_2D(i)
            ring2.append([lon, lat])
        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
    wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    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

artDataBankenURL = 'https://api.artdatabanken.se'
artDataBankenHeaders = {
    'X-Api-Version': '1.5',
    'X-Requesting-System': '',
    # https://github.com/biodiversitydata-se/SOS/blob/master/Docs/Authentication.md
    'Ocp-Apim-Subscription-Key': config['ArtDataBanken']['Subscription-Key']
}

def getObservations(taxonLists, taxonRedlistCategories, searchFilter):
    # https://api-portal.artdatabanken.se/docs/services/sos-api-v1/operations/Exports_DownloadGeoJson
    headers = artDataBankenHeaders.copy()
    headers['Content-Type'] = 'application/json'
    headers['Authorization'] = config['ArtDataBanken']['Authorization']
    params = {
        'cultureCode': 'sv-SE',
        'flat': 'true',
        'excludeNullValues': 'true',
        'gzip': 'false',
        #'sensitiveObservations': 'true',
        'validateSearchFilter': 'true',
    }
    resp = requests.post(
        artDataBankenURL + '/species-observation-system/v1/Exports/Download/GeoJson',
        headers=headers,
        params=params,
        data=json.dumps(searchFilter)
    )

    resp.raise_for_status()
    obs = resp.json()

    # https://www.rfc-editor.org/rfc/rfc7946
    if obs is None or type(obs) != dict or 'type' not in obs.keys():
        raise Exception('invalid GeoJSON output')
    if 'crs' not in obs.keys():
        print('WARN: GeoJSON output lacks CRS', file=sys.stderr)

    layername = 'Observations'
    if obs['type'] == 'FeatureCollection' and 'features' in obs.keys() and type(obs['features']) == list:
        print(f'{len(obs["features"])} observations found', file=sys.stderr)
        if 'name' not in obs.keys() or obs['name'] is None:
            obs['name'] = layername
        for feat in obs['features']:
            if (type(feat) != dict or 'type' not in feat.keys() or feat['type'] != 'Feature'
                    or 'properties' not in feat.keys() or type(feat['properties']) != dict):
                print('WARN: Invalid feature in GeoJSON output', file=sys.stderr)
                continue
            properties = feat['properties']
            if 'DyntaxaTaxonId' not in properties.keys() or properties['DyntaxaTaxonId'] is None:
                print('WARN: Feature lacks taxon ID', file=sys.stderr)
                continue

            tid = properties['DyntaxaTaxonId']
            for k, taxonList in taxonLists.items():
                v = (tid in taxonList)
                if k in properties.keys() and properties[k] != v:
                    print(f'WARN: #{tid} {k}: {properties[k]} → {v}', file=sys.stderr)
                properties[k] = v

            if tid in taxonRedlistCategories.keys():
                k = 'RedlistCategory'
                v = taxonRedlistCategories[tid]
                if k in properties.keys() and properties[k] != v:
                    print(f'WARN: #{tid} {k}: {properties[k]} → {v}', file=sys.stderr)
                properties[k] = v

    drv = ogr.GetDriverByName(args.observation_format)
    if drv.name == geojson_drv.name:
        path = args.project_home.joinpath(layername.lower() + '.geojson')
        with path.open(mode='w') as fp:
            json.dump(obs, fp, indent=2, ensure_ascii=False)
    else:
        with tempfile.NamedTemporaryFile(suffix='.geojson', mode='w') as fp:
            json.dump(obs, fp)
            fp.flush()

            suffix = drv.GetMetadataItem(gdal.DMD_EXTENSIONS).split(' ')[0]
            path = args.project_home.joinpath(layername.lower()).with_suffix('.' + suffix)
            gdal.VectorTranslate(path.as_posix(), fp.name, reproject=True, dstSRS=target_srs, format=drv.name)
    obs = None

    if projectInstance is not None:
        query = '''"IsPositiveObservation" AND "IsNaturalOccurrence" AND "CoordinateUncertaintyInMeters" <= 250'''
        layer = QgsVectorLayer(path.as_posix() + f'|subset={query}', 'Fynd', 'ogr')
        if not layer.isValid():
            raise Exception(f'ERROR: {path}: failed to load in QGIS')
        if args.style is not None:
            layer.loadNamedStyle(args.style.as_posix())
        layer.setReadOnly(True)
        layer.setFlags(QgsMapLayer.Identifiable | QgsMapLayer.Searchable)
        projectInstance.addMapLayer(layer)
        layerTreeRoot.findLayer(layer.id()).setItemVisibilityChecked(True)

def getTaxonLists():
    resp = requests.get(
        artDataBankenURL + '/species-observation-system/v1/TaxonLists?cultureCode=en-US',
        headers=artDataBankenHeaders
    )
    resp.raise_for_status()
    taxonLists = resp.json()
    #print(json.dumps(taxonLists, indent=2))
    taxonLists2 = {}

    # https://github.com/biodiversitydata-se/SOS/blob/master/Src/SOS.lib/Resources/TaxonLists.json
    # https://artfakta.se/artinformation/taxa
    getTaxonList(taxonLists, taxonLists2, 'ProtectedByLaw', 'Protected by law species')
    getTaxonList(taxonLists, taxonLists2, 'TaxonIsSignalSpecie', 'Signal species')

    getTaxonList(taxonLists, taxonLists2, 'TaxonIsInvasive', 'Invasive species')
    getTaxonList(taxonLists, taxonLists2, 'TaxonIsInvasiveInSweden', 'Invasive species in Sweden')
    getTaxonList(taxonLists, taxonLists2, 'TaxonIsInvasiveEuRegulation', 'EU regulation 1143/2014')

    getTaxonList(taxonLists, taxonLists2, 'Natura2000HabitatsDirective', 'Habitats directive species')
    getTaxonList(taxonLists, taxonLists2, 'Natura2000HabitatsDirectiveArticle2', 'Habitats directive Annex 2')
    getTaxonList(taxonLists, taxonLists2, 'Natura2000HabitatsDirectiveArticle2PrioritySpecie', 'Habitats directive Annex 2, priority species')
    getTaxonList(taxonLists, taxonLists2, 'Natura2000HabitatsDirectiveArticle4', 'Habitats directive Annex 4')
    getTaxonList(taxonLists, taxonLists2, 'Natura2000HabitatsDirectiveArticle5', 'Habitats directive Annex 5')

    # sanity check
    directives = [
        'Natura2000HabitatsDirectiveArticle2',
        'Natura2000HabitatsDirectiveArticle2PrioritySpecie',
        'Natura2000HabitatsDirectiveArticle4',
        'Natura2000HabitatsDirectiveArticle5'
    ]
    d0 = 'Natura2000HabitatsDirective'
    for i in taxonLists2[d0]:
        if not any(map(lambda d: i in taxonLists2[d], directives)):
            raise Exception(f'missing taxon #{i} (in {d0}) from directives {",".join(directives)}')
    for d in directives:
        for i in taxonLists2[d]:
            if i not in taxonLists2[d0]:
                raise Exception(f'missing taxon #{i} (in {d}) from {d0}')
    for i in taxonLists2['Natura2000HabitatsDirectiveArticle2PrioritySpecie']:
        if i not in taxonLists2['Natura2000HabitatsDirectiveArticle2']:
            raise Exception(f'missing taxon #{i} (in Natura2000HabitatsDirectiveArticle2PrioritySpecie) from Natura2000HabitatsDirectiveArticle2')

    getTaxonList(taxonLists, taxonLists2, 'BirdDirective', 'Birds Directive')
    getTaxonList(taxonLists, taxonLists2, 'BirdDirectiveArticle1', 'Birds directive - Annex 1')
    getTaxonList(taxonLists, taxonLists2, 'BirdDirectiveArticle2', 'Birds directive - Annex 2')
    getTaxonList(taxonLists, taxonLists2, 'TaxonIsPriorityBird', 'Priority birds')

    directives = ['BirdDirectiveArticle1', 'BirdDirectiveArticle2']
    d0 = 'BirdDirective'
    #for i in taxonLists2[d0]:
    #    if not any(map(lambda d: i in taxonLists2[d], directives)):
    #        raise Exception(f'missing taxon #{i} (in {d0}) from directives {",".join(directives)}')
    for d in directives:
        for i in taxonLists2[d]:
            if i not in taxonLists2[d0]:
                raise Exception(f'missing taxon #{i} (in {d}) from {d0}')

    i = getTaxonList(taxonLists, taxonLists2, 'TaxonIsRedlisted', 'Redlisted species')
    taxonRedlistCategories = getTaxonRedlistCategories(taxonLists, i)

    return taxonLists2, taxonRedlistCategories

def getTaxonList(taxonLists, taxonLists2, key, name):
    i = None
    for t in taxonLists:
        if t['name'] == name:
            if i is not None:
                raise Exception(f'multiple results found for taxon list "{name}"')
            i = t['id']
    if i is None:
        raise Exception(f'no found for taxon list "{name}"')

    resp = requests.get(
        artDataBankenURL + f'/species-observation-system/v1/TaxonLists/{i}/Taxa',
        headers=artDataBankenHeaders
    )
    resp.raise_for_status()
    resp = resp.json()
    if type(resp) != list:
        raise Exception(f'expected list, got {type(resp)}')
    taxonLists2[key] = {r['id'] for r in resp}

    return i

def getTaxonRedlistCategories(taxonLists, i):
    taxonRedlistCategories = {}
    r = re.compile(r'\(([A-Z][A-Z])\)\Z')
    for t in taxonLists:
        if 'parentId' not in t.keys() or t['parentId'] != i:
            continue
        name = t['name']
        c = r.search(name)
        if not c:
            raise Exception(f'invalid redlist name "{name}"')
        c = c.group(1)

        resp = requests.get(
            artDataBankenURL + f'/species-observation-system/v1/TaxonLists/{t["id"]}/Taxa',
            headers=artDataBankenHeaders
        )
        resp.raise_for_status()
        resp = resp.json()
        if type(resp) != list:
            raise Exception(f'expected list, got {type(resp)}')
        for t in resp:
            j = t['id']
            if j in taxonRedlistCategories.keys():
                raise Exception(f'duplicate redlist classification for taxon #{j}')
            taxonRedlistCategories[j] = c
    return taxonRedlistCategories

geograficsFilter = {
    'geometries': geometricFilter(geometries),
    'maxAccuracy': 2500,
    'considerObservationAccuracy': True,
    'considerDisturbanceRadius': False
}
taxonLists, taxonRedlistCategories = getTaxonLists()

# https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md
searchFilter = {
    'geographics': geograficsFilter,
    'determinationFilter': 'NoFilter',
    'notRecoveredFilter': 'DontIncludeNotRecovered',
    'occurrenceStatus': 'Present',
    'verificationStatus': 'BothVerifiedAndNotVerified',
    'output': {
        'fields': [
            # https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md#fields
            # https://github.com/biodiversitydata-se/SOS/blob/master/Docs/Observation.md
            'dataProviderId',
            'datasetName',
            #'basisOfRecord',
            'rightsHolder',

            'modified',
            'sensitive',
            'measurementOrFacts',
            'projects',

            'occurrence.occurrenceId',
            'occurrence.occurrenceRemarks',
            'occurrence.recordedBy',
            'occurrence.reportedBy',
            'occurrence.reportedDate',
            'occurrence.occurrenceStatus',
            'occurrence.activity',
            'occurrence.behavior',
            'occurrence.biotope',
            'occurrence.biotopeDescription',
            'occurrence.lifeStage',
            'occurrence.reproductiveCondition',
            'occurrence.sex',
            'occurrence.associatedMedia',
            'occurrence.associatedReferences',
            'occurrence.individualCount',
            'occurrence.organismQuantityInt',
            'occurrence.organismQuantityUnit',
            'occurrence.sensitivityCategory',
            'occurrence.isNaturalOccurrence',
            'occurrence.isNeverFoundObservation',
            'occurrence.isNotRediscoveredObservation',
            'occurrence.isPositiveObservation',
            'occurrence.substrate.description',
            'occurrence.substrate.name',
            'occurrence.length',
            'occurrence.weight',
            'occurrence.url',

            'event.startDate',
            'event.endDate',
            'event.habitat',
            'event.eventRemarks',
            'event.discoveryMethod',
            'event.measurementOrFacts',
            'event.fieldNotes',

            'identification.verified',
            'identification.uncertainIdentification',
            'identification.verificationStatus',
            'identification.confirmedBy',
            'identification.confirmedDate',
            'identification.identifiedBy',
            'identification.dateIdentified',
            'identification.verifiedBy',
            'identification.determinationMethod',
            'identification.identificationRemarks',

            'location.locality',
            'location.county',
            'location.municipality',
            'location.coordinateUncertaintyInMeters',

            'taxon.id',
            'taxon.scientificName',
            'taxon.scientificNameAuthorship',
            'taxon.vernacularName',
            'taxon.genus',
            'taxon.family',
            'taxon.order',
            'taxon.class',
            'taxon.phylum',
            'taxon.kingdom',

            'taxon.attributes.taxonCategory',
            'taxon.attributes.organismGroup',

            'taxon.attributes.redlistCategory',
        ]
    }
}

# https://artfakta.se/artinformation/taxa/biota-0/detaljer
kingdomIds = {
    'Animalia': 5000001,
    'Archaea': 5000082,
    'Bacteria': 5000052,
    'Chromista': 5000055,
    'Fungi': 5000039,
    'Plantae': 5000045,
    'Protozoa': 5000060,
    'Viruses': 5000083,
    'Algae': 6001047,
}

searchFilter['taxon'] = {
    'ids': list(map(lambda k: kingdomIds[k], ['Fungi', 'Plantae'])),
    'includeUnderlyingTaxa': True,
    # https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md#taxon-lists
    'taxonListIds': [1, 2, 7, 8, 13, 18],
    'taxonListOperator': 'Filter'
}

getObservations(taxonLists, taxonRedlistCategories, searchFilter)

topo_maps = ['Topografi 10', 'Topografi 50', 'Topografi 100']
if args.topo_basedir is not None:
    for topo in topo_maps:
        qlr_paths = []
        d = args.topo_basedir.joinpath(topo)
        if not d.is_dir():
            continue

        if args.project_home is not None:
            l = args.project_home.expanduser().joinpath(topo)
            l.symlink_to(d, target_is_directory=True)
            d = l

        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 projectInstance is not None and 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)
            rootGroup.setExpanded(False)

    if projectInstance is not None:
        maxScale = 1
        for topo in topo_maps:
            minScale = None
            if topo == 'Topografi 10':
                minScale = 50000
            elif topo == 'Topografi 50':
                minScale = 100000
            elif topo == 'Topografi 100':
                minScale = 10000000

            if minScale is not None:
                for lt in QgsProject.instance().layerTreeRoot().findGroup(topo).findLayers():
                    lyr = lt.layer()
                    lyr.setScaleBasedVisibility(True)
                    lyr.setMinimumScale(minScale)
                    lyr.setMaximumScale(maxScale)
                maxScale = minScale

if projectInstance is not None:
    projectInstance.write()