aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xgis-observation-map112
1 files changed, 55 insertions, 57 deletions
diff --git a/gis-observation-map b/gis-observation-map
index 8fc48de..a1ce006 100755
--- a/gis-observation-map
+++ b/gis-observation-map
@@ -208,15 +208,8 @@ if projectInstance is not None:
else:
projectInstance.addMapLayer(layer, False)
sourceGroup.addLayer(layer)
- if i == 0:
- #layerTreeModel = QgsLayerTreeModel(layerTreeRoot)
- #layerTreeView = QgsLayerTreeView()
- #layerTreeView.setModel(layerTreeModel)
- #layerTreeView.setCurrentLayer(layer)
- canvas = QgsMapCanvas()
- canvas.setExtent( layer.extent() )
-
-def gpkg_intersects(path, geometries):
+
+def gpkg_intersects(path, source_geometry):
ds = gpkg_drv.Open(path, update=0)
for i in range(ds.GetLayerCount()):
layer = ds.GetLayerByIndex(i)
@@ -227,14 +220,13 @@ def gpkg_intersects(path, geometries):
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
+ if geometry.Intersects(source_geometry):
+ ds = None
+ return True
feature = layer.GetNextFeature()
return False
-def idx_intersects(idx_path, geometries):
+def idx_intersects(idx_path, source_geometry):
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)
@@ -250,23 +242,16 @@ def idx_intersects(idx_path, geometries):
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
+ if geometry.Intersects(source_geometry):
+ loc = feature.GetField(locFieldIdx)
+ locations.append(loc)
if union is None:
union = geometry.Clone()
else:
union = union.Union(geometry)
feature = layer.GetNextFeature()
- contains = True
- for source_geometry in geometries:
- if not union.Contains(source_geometry):
- contains = False
- break
-
+ contains = union.Contains(source_geometry)
ds = None
return locations, contains
@@ -279,8 +264,9 @@ def find_qlr(basedir):
return None
return qlr[0]
-geometries = []
+clipGeometry = ogr.Geometry(ogr.wkbMultiPolygon)
for pt in args.point:
+ # northing,easting in target SRS
y, x = pt.split(',', 1)
x = int(x)
y = int(y)
@@ -300,7 +286,7 @@ for pt in args.point:
bbox = ogr.Geometry(ogr.wkbPolygon)
bbox.AddGeometry(ring)
- geometries.append(bbox)
+ clipGeometry = clipGeometry.Union(bbox)
for geom in args.geometry:
ds = ogr.Open(geom['path'].as_posix(), update=0)
@@ -314,12 +300,12 @@ for geom in args.geometry:
geometry = feature.GetGeometryRef()
if not geometry.IsValid():
raise Exception('Invalid geometry')
+ geometryType = geometry.GetGeometryType()
+ geometry = geometry.Clone()
+ geometry.FlattenTo2D()
geometry.Transform(transform_srs)
- if args.margin is None:
- geometry = geometry.Clone()
- geometry.FlattenTo2D()
- else:
+ if args.margin is not None:
minX, maxX, minY, maxY = geometry.GetEnvelope()
minX -= args.margin
maxX += args.margin
@@ -335,15 +321,27 @@ for geom in args.geometry:
bbox = ogr.Geometry(ogr.wkbPolygon)
bbox.AddGeometry(ring)
- geometry = bbox
+ clipGeometry = clipGeometry.Union(bbox)
+ elif 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:
+ clipGeometry = clipGeometry.Union(geometry)
+ elif geometryType == ogr.wkbMultiPolygon:
+ for i in range(geometry.GetGeometryCount()):
+ geom2 = geometry.GetGeometryRef(i)
+ clipGeometry = clipGeometry.Union(geom2)
+ else:
+ raise Exception('Unsuported geometry (try --margin=0)')
- geometries.append(geometry)
feature = layer.GetNextFeature()
-if len(args.geometry) == 0 and projectInstance is not None:
- r = geometries[0].GetGeometryRef(0)
- pt0X, pt0Y = r.GetPoint_2D(0)
- pt1X, pt1Y = r.GetPoint_2D(2)
+if not clipGeometry.IsValid():
+ raise Exception('Invalid geometry')
+
+if projectInstance is not None:
+ pt0X, pt1X, pt0Y, pt1Y = clipGeometry.GetEnvelope()
canvas = QgsMapCanvas()
extent = QgsRectangle(QgsPointXY(pt0X, pt0Y), QgsPointXY(pt1X, pt1Y))
canvas.setExtent(extent)
@@ -372,31 +370,31 @@ def geometricFilterPolygon(polygon):
}
# https://github.com/biodiversitydata-se/SOS/blob/master/Docs/SearchFilter.md#geographics-filter
-def geometricFilter(geometries):
+def geometricFilter(geometry):
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326) # WGS84
wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
towgs84 = osr.CoordinateTransformation(target_srs, wgs84)
+ geometry = geometry.Clone()
+ geometry.Transform(towgs84)
+ geometryType = geometry.GetGeometryType()
+
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:
+ 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(geometry)
+ myGeometries.append(gfp)
+ elif geometryType == ogr.wkbMultiPolygon:
+ for i in range(geometry.GetGeometryCount()):
+ geom = geometry.GetGeometryRef(i)
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)')
+ else:
+ raise Exception('Unsuported geometry (try --margin=0)')
return myGeometries
@@ -686,7 +684,7 @@ def getTaxonRedlistCategories(taxonLists, i):
if not args.no_observations:
geographicsFilter = {
- 'geometries': geometricFilter(geometries),
+ 'geometries': geometricFilter(clipGeometry),
'maxAccuracy': 2500,
'considerObservationAccuracy': True,
'considerDisturbanceRadius': False
@@ -834,7 +832,7 @@ if args.topo_basedir is not None and projectInstance is not None:
idx = d.joinpath('index.geojson')
if idx.is_file():
- tiles, idxContains[topo] = idx_intersects(idx.as_posix(), geometries)
+ tiles, idxContains[topo] = idx_intersects(idx.as_posix(), clipGeometry)
if len(tiles) == 0:
continue
@@ -887,7 +885,7 @@ if args.topo_basedir is not None and projectInstance is not None:
print(f'WARN: {t}: lacks {gpkg_name}, skipping', file=sys.stderr)
continue
- if not gpkg_intersects(gpkg_path.as_posix(), geometries):
+ if not gpkg_intersects(gpkg_path.as_posix(), clipGeometry):
continue
qlr_paths.append(qlr_path)