''' @author: olivier.massot, 2018 ''' from osgeo import ogr, osr from path import Path GEOM_UNKNOWN = 0 GEOM_POINT = 1 GEOM_LINE = 2 GEOM_POLYGON = 3 GEOM_MULTIPOINT = 4 GEOM_MULTILINE = 5 GEOM_MULTIPOLYGON = 6 GEOM_NAMES = {0: "(AUCUN)", 1: "POINT", 2: "LIGNE", 3: "POLYGONE", 4: "MULTI-POINT", 5:"MULTI-LIGNE", 6:"MULTI-POLYGONE", } def point_(x, y): point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(x, y) return point class Datasource(): DRIVER_NAME = "ESRI Shapefile" def __init__(self, filename, readonly=True): self.filename = Path(filename) driver = ogr.GetDriverByName(self.DRIVER_NAME) self._ogr_datasource = driver.Open(filename, 0 if readonly else 1) if not self._ogr_datasource: raise IOError("Unable to read the file {}".format(filename)) @property def layer(self): return Layer(self._ogr_datasource.GetLayer()) class LayerField(): def __init__(self, ogr_field): self.name = ogr_field.GetName() self.type_ = ogr_field.GetType() self.type_name = ogr_field.GetFieldTypeName(self.type_) self.width = ogr_field.GetWidth() self.precision = ogr_field.GetPrecision() def __repr__(self): return "".format(self.name, self.type_name, self.width) class Layer(): def __init__(self, ogr_layer): self._ogr_layer = ogr_layer self.name = ogr_layer.GetName() self.srid = ogr_layer.GetSpatialRef().GetAttrValue('AUTHORITY',1) self.geom_type = ogr_layer.GetGeomType() _layer_def = self._ogr_layer.GetLayerDefn() self.fields = [LayerField(_layer_def.GetFieldDefn(i)) for i in range(_layer_def.GetFieldCount())] @property def geom_name(self): return GEOM_NAMES[self._ogr_layer.GetGeomType()] def __len__(self): return self._ogr_layer.GetFeatureCount() def __iter__(self): for f in self._ogr_layer.__iter__(): yield Feature(f) self._ogr_layer.ResetReading() def __getitem__(self, i): return Feature(self._ogr_layer.GetFeature(i)) class Feature(): def __init__(self, ogr_feature): self._ogr_feature = ogr_feature self.geom = self._ogr_feature.GetGeometryRef() self.geom_type = self.geom.GetGeometryType() xmin, xmax, ymin, ymax = self.geom.GetEnvelope() self.bounding_box = (xmin, ymin, xmax, ymax) _def = ogr_feature.GetDefnRef() for i in range(_def.GetFieldCount()): field = _def.GetFieldDefn(i) value = ogr_feature[i] if value is None and field.GetType() == ogr.OFTString: value = "" self.__setattr__(field.GetName(), value) @property def geom_name(self): return GEOM_NAMES[self.geom.GetGeometryType()] def fields(self): layerDefinition = self._ogr_layer.GetLayerDefn() return [LayerField(layerDefinition.GetFieldDefn(i)) for i in range(layerDefinition.GetFieldCount())] @property def valid(self): return self.geom.IsValid() def __repr__(self): return "".format(self.geom_name, self.bounding_box) @property def points(self): if self.geom_type != GEOM_POLYGON: return [point_(*p) for p in self.geom.GetPoints()] else: return [point_(*p) for p in self.geom.GetGeometryRef(0).GetPoints()] def almost_equal(self, feature, buffer_dist=1): if feature.geom_type == self.geom_type: buffer = self.geom.Buffer(buffer_dist) if buffer.Contains(feature.geom): return True return False @classmethod def union(cls, features): g = features[0].geom.Clone() for f in features[1:]: g = g.Union(f.geom) return g @classmethod def buffered_union(cls, features, buffer_dist=1): multi = ogr.Geometry(ogr.wkbMultiPolygon) for f in features: geom = f.geom if issubclass(f.__class__, Feature) else f # au cas où ce ne seraient pas des features, mais directement des geometries multi.AddGeometry(geom.Buffer(buffer_dist)) return multi.UnionCascaded() if __name__ == "__main__": ds = Datasource(r"C:\dev\python\datacheck\work\STURNO_178AP1_REC_171121_OK\ARTERE_GEO.shp") layer = ds.layer print(layer.fields) # for f in layer: # print(f) f = layer[4] f2 = layer[5] f3 = layer[4] print(f) print(f.geom) print(f.AR_ID_INSE) gf = Feature.union(list(layer)) print(gf) # gf = Feature.buffered_union(list(layer), 0.5) # print(gf) # save to file result_file = Path(r"C:\dev\python\datacheck\work\result.shp") driver = ogr.GetDriverByName("ESRI Shapefile") if result_file.exists(): driver.DeleteDataSource(result_file) data_source = driver.CreateDataSource(result_file) srs = osr.SpatialReference() srs.ImportFromEPSG(3949) layer = data_source.CreateLayer("result", srs, gf.GetGeometryType()) feature = ogr.Feature(layer.GetLayerDefn()) feature.SetGeometry(gf) layer.CreateFeature(feature)