diff options
| -rwxr-xr-x | gis-observation-map | 112 |
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) |