aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xgis-observation-map391
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()