I have a GeoPackage with severallayers, and a folder of shapefiles showing boundaries that define “holes.” in gpkg. Because I filled those holes in gpkg with other features, some features that should be one continuous feature have been split apart.
My goal is:
- For each polygon in each boundary shapefile, buffer it by 5 m.
- In every layer of the GeoPackage, select all features intersecting that buffer.
- Dissolve those selected features by a set of attribute fields
- Delete the originals and write the dissolved features back to the GeoPackage, so that any segments “near” the hole become one feature again.
What I’m seeing
for example, When i select a line feature intersecting a polygon, another line feature that intersects another polygon, that is very far away, also get selected. Which means those have been dissolved into multipart features, even though they are in different polygon intersections. How could this happen>? How can I restructure this workflow (in OGR/Shapely) so that each polygon’s dissolve only affects the features that genuinely intersect that polygon, without “leakage” into other areas?
def dissolve_pol(flattened_gpkg, shp_dir, unique_fields): source = ogr.Open(flattened_gpkg, 1) shp_list = [f for f in os.listdir(shp_dir) if f.endswith('.shp')] for i in range(source.GetLayerCount()): layer = source.GetLayerByIndex(i) layer_name = layer.GetName() if layer_name in layers_to_skip: continue for shp_file in shp_list: shp_path = os.path.join(shp_dir, shp_file) shp_ds = ogr.GetDriverByName('ESRI Shapefile').Open(shp_path) shp_layer = shp_ds.GetLayer() print(f" SHP: {shp_file}") logging.info(f" SHP: {shp_file}") for shp_feat in shp_layer: poly_id = shp_feat.GetFID() shp_geom = shp_feat.GetGeometryRef().Clone() buffer_geom = shp_geom.Buffer(1) # small buffer for safety intersect_feats = [] layer.ResetReading() for feat in layer: geom = feat.GetGeometryRef() if not geom.IsValid(): logging.info(f"Invalid geometry in {layer_name}, fixing...") geom = geom.MakeValid() if not geom.IsValid(): continue if geom is not None and geom.Intersects(buffer_geom): intersect_feats.append(feat.Clone()) if not intersect_feats: continue dissolved_dict = {} layer_defn = layer.GetLayerDefn() field_names = [layer_defn.GetFieldDefn(i).GetName() for i in range(layer_defn.GetFieldCount())] usable_fields = [f for f in unique_fields if f in field_names] for feat in intersect_feats: key = [] for f in usable_fields: value = feat.GetField(f) value = '' if value is None else str(value).strip().upper().replace('\n','').replace('\t','') key.append(value) key = tuple(key) geom = feat.GetGeometryRef().Clone() if key not in dissolved_dict: dissolved_dict[key] = geom else: dissolved_dict[key] = dissolved_dict[key].Union(geom) # Remove old intersecting features fid_list = [f.GetFID() for f in intersect_feats] for fid in fid_list: layer.DeleteFeature(fid) # Add new dissolved features for key, geom in dissolved_dict.items(): new_feat = ogr.Feature(layer.GetLayerDefn()) for i, field in enumerate(usable_fields): new_feat.SetField(field, key[i]) new_feat.SetGeometry(geom) layer.CreateFeature(new_feat) new_feat = None # After this polygon → reset reading to avoid stale state layer.ResetReading() shp_ds = None layer.ResetReading() source = None How can I restructure this workflow (in OGR/Shapely) so that each polygon’s dissolve only affects the features that genuinely intersect that polygon, without “leakage” into other areas?
