Fractionner le fichier de formes par fonctionnalité en Python à l'aide de GDAL?

est-il possible de diviser un fichier de formes par fonctionnalité en python? (Le mieux serait une solution où je peux enregistrer temporairement les objets vectoriels résultants dans la mémoire à la place sur le disque).

La raison: je veux utiliser la fonction gdal rasterizeLayer avec plusieurs sous-ensembles différents de fichiers de formes. La fonction nécessite un objet osgeo.ogr.Layer.

mkay, j'ai essayé un peu et cela pourrait fonctionner comme suit. Vous pouvez obtenir la géométrie des objets de couche gdal par entité comme suit.

#  Load shape into gdal
shapefile=str(vectorPath)
layer_source = ogr.Open(shapefile)
lyr = layer_source.GetLayer(0)
for i in range(0,lyr.GetFeatureCount()):
     feat = lyr.GetFeature(i)
     ge = feat.geometry()

Maintenant, j'ai juste besoin de savoir comment créer un objet osgeo.ogr.layer basé sur cette géométrie.

Pour clarification. J'ai besoin d'une fonction en simple code ogr / gdal! Cela semble également intéresser d'autres personnes et je veux toujours une solution sans modules secondaires (bien que toute solution venant d'ici sera utilisée dans un plugin qgis disponible gratuitement).

OK donc une deuxième tentative pour répondre à votre question avec une solution GDAL pure.

Tout d'abord, GDAL (Geospatial Data Abstraction Library) n'était à l'origine qu'une bibliothèque pour travailler avec des données géospatiales raster, tandis que la bibliothèque OGR séparée était destinée à fonctionner avec des données vectorielles. Cependant, les deux bibliothèques sont maintenant partiellement fusionnées et sont généralement téléchargées et installées ensemble sous le nom combiné de GDAL. La solution relève donc vraiment de l'OGR. Vous l'avez dans votre code initial, donc je suppose que vous le saviez, mais c'est une distinction importante à retenir lors de la recherche de conseils et d'astuces.

Pour lire les données d'une couche vectorielle, votre code initial est très bien:

from osgeo import ogr
shapefile = ogr.Open(shapefile)
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    print i, name, geometry.GetGeometryName()

Nous devons créer une nouvelle fonctionnalité avant de pouvoir l'écrire dans un fichier de formes (ou tout autre ensemble de données vectorielles). Pour créer une nouvelle entité, nous avons d'abord besoin: - d'une géométrie - d'une définition d'entité, qui inclura probablement des définitions de champ. Utilisez le constructeur Geometry ogr.Geometry () pour créer un objet Geometry vide. Définissez la géométrie d'une manière différente pour chaque type (point, ligne, polygone, etc.). Ainsi, par exemple:

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20)

ou

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
line.SetPoint(0,30,30) #(10,10) -> (30,30)

Pour une définition de champ

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)

Vous pouvez maintenant créer votre couche vectorielle. Dans ce cas, un polygone carré:

#create simple square polygon shapefile:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')

datasource = driver.CreateDataSource('YOUR_PATH')
layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbPolygon)

#create polygon object:
myRing = ogr.Geometry(type=ogr.wkbLinearRing)
myRing.AddPoint(0.0, 0.0)     #LowerLeft
myRing.AddPoint(0.0, 10.0)    #UpperLeft
myRing.AddPoint(10.0, 10.0)   #UpperRight
myRing.AddPoint(10.0, 0.0)    #Lower Right
myRing.AddPoint(0.0, 0.0)     #close ring
myPoly = ogr.Geometry(type=ogr.wkbPolygon)
myPoly.AddGeometry(myRing)
print ('Polygon area =',myPoly.GetArea())  #returns correct area of 100.0

#create feature object with point geometry type from layer object:
feature = ogr.Feature( layer.GetLayerDefn())
feature.SetGeometry(myPoly)
layer.CreateFeature(feature)

#flush memory - very important
feature.Destroy()
datasource.Destroy()

J'ai eu de la chance de lire et d'écrire sur des calques. Plus précisément, j'ai du code qui lira une couche de fichiers de formes contenant des polylignes et affichera la géométrie de chaque entité dans des fichiers texte (utilisée comme entrée pour un ancien modèle).

name     = layer.name()
provider = layer.dataProvider()
feat     = QgsFeature()

# Now we can loop through all the defined features
while provider.nextFeature(feat):

    # Get layer attributes               
    attrs = feat.attributeMap()
    for (k,attr) in attrs.iteritems():
        if k == 0:
            attrOne = attr.toString()
        elif k == 1:
            attrTwo = attr.toString()
        ...

    # Gets the geometry of the feature
    geom = feat.geometry()

    # Get the coordinates of the whole line [or use asPoint()]                    
    line = geom.asPolyline()
        # all points in the line
        for point in line:
            lat = point[0]
            lon = point[1]
            # Add these to a QgsGeometry
            your_Own_QgsGeometry.add...

Cela semble utile pour obtenir chacune des fonctionnalités de vos calques.

L'écriture dans un autre calque ne devrait pas être trop complexe à partir d'ici. Quelque chose comme ça devrait fonctionner en théorie:

# New layer name
filename = "myNewLayer.shp"

# define fields for feature attributes
fields   = { 0 : QgsField("attrOne", QVariant.String),
             1 : QgsField("attrTwo", QVariant.String),
             2 : QgsField("...", QVariant.Int) }

# Create coordinate reference system as WGS84
crs    = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)

# Create the vector layer
writer = QgsVectorFileWriter(filename, "CP1250", fields, QGis.WKBLineString, crs)

# Create some features
feat = QgsFeature()
feat.addAttribute(0, QVariant(runway))
feat.addAttribute(1, QVariant(arriveDepart))
feat.addAttribute(2, QVariant(subTrack))

# Add your geometry
feat.setGeometry(your_Own_QgsGeometry)

# Add the features
writer.addFeature(feat)

# Add it to QGIS project
self.iface.addVectorLayer(filename, "MyLayerName", "ogr")

À partir de là, vous devriez pouvoir obtenir des données de chaque entité et écrire de nouvelles entités dans une nouvelle couche.

Dan