datacheck/core/gis_.py

'''

@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 "<Field: nom={}, type={}, longueur={}>".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 "<Feature: type={}, bounding_box={}>".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)