Comment créer un fichier de formes polygonales à partir d'une liste de coordonnées en utilisant python gdal / ogr?

J'essaie de créer un fichier de formes polygonales à partir d'une liste de coordonnées en utilisant des outils open source python. Ce qui suit est ce que j'ai jusqu'à présent qui a été piraté ensemble à partir du livre de recettes Python GDAL / OGR et ce GIS SE répond .

Il y a une question similaire Python: Comment créer un Polygon Shapefile à partir d'une liste de coordonnées X, Y? , bien que cette question se réfère à l'utilisation pyshp. Cependant, je suis intéressé par la création d'un fichier de formes de polygones en utilisant uniquement les outils Python gdal / ogr.

import ogr

def create_polygon(coords):          
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for coord in coords:
        ring.AddPoint(coord[0], coord[1])

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly.ExportToWkt()

def write_shapefile(poly, out_shp):
    """
    https://gis.stackexchange.com/a/52708/8104
    """
    # Now convert it to a shapefile with OGR    
    driver = ogr.GetDriverByName('Esri Shapefile')
    ds = driver.CreateDataSource(out_shp)
    layer = ds.CreateLayer('', None, ogr.wkbPolygon)
    # Add one attribute
    layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
    defn = layer.GetLayerDefn()

    ## If there are multiple geometries, put the "for" loop here

    # Create a new feature (attribute and geometry)
    feat = ogr.Feature(defn)
    feat.SetField('id', 123)

    # Make a geometry, from Shapely object
    geom = ogr.CreateGeometryFromWkb(poly.wkb)
    feat.SetGeometry(geom)

    layer.CreateFeature(feat)
    feat = geom = None  # destroy these

    # Save and close everything
    ds = layer = feat = geom = None

def main(coords, out_shp):
    poly = create_polygon(coords)
    write_shapefile(poly, out_shp)

if __name__ == "__main__":
    coords = [(-106.6472953, 24.0370137), (-106.4933356, 24.05293569), (-106.4941789, 24.01969175), (-106.4927777, 23.98804445), (-106.4922614, 23.95582128), (-106.4925834, 23.92302327), (-106.4924068, 23.89048039), (-106.4925695, 23.85771361), (-106.4932479, 23.82457675), (-106.4928676, 23.7922049), (-106.4925072, 23.75980241), (-106.492388, 23.72722475), (-106.4922574, 23.69464296), (-106.4921181, 23.6620529), (-106.4922734, 23.62926733), (-106.4917201, 23.59697561), (-106.4914134, 23.56449628), (-106.4912558, 23.5319045), (-106.491146, 23.49926362), (-106.4911734, 23.46653561), (-106.4910181, 23.43392476), (-106.4910156, 23.40119976), (-106.4909501, 23.3685223), (-106.4908165, 23.33586566), (-106.4907735, 23.30314904), (-106.4906954, 23.27044931), (-106.4906366, 23.23771759), (-106.4905894, 23.20499124), (-106.4905432, 23.17226022), (-106.4904748, 23.1395177), (-106.4904187, 23.10676788), (-106.4903676, 23.07401321), (-106.4903098, 23.04126832), (-106.4902512, 23.00849426), (-106.4901979, 22.97572025), (-106.490196, 22.97401001), (-106.6481193, 22.95609832), (-106.6481156, 22.95801668), (-106.6480697, 22.99082052), (-106.6480307, 23.02362441), (-106.6479937, 23.0563977), (-106.6479473, 23.0891833), (-106.647902, 23.12196713), (-106.6478733, 23.15474057), (-106.6478237, 23.18750353), (-106.6477752, 23.22026138), (-106.6477389, 23.25302505), (-106.647701, 23.28577123), (-106.6476562, 23.31851549), (-106.6476211, 23.3512557), (-106.6475745, 23.38397935), (-106.6475231, 23.41671055), (-106.6474863, 23.44942382), (-106.6474432, 23.48213255), (-106.6474017, 23.51484861), (-106.6474626, 23.54747418), (-106.647766, 23.57991134), (-106.6482374, 23.61220905), (-106.6484783, 23.64467084), (-106.6482308, 23.6775148), (-106.6479338, 23.7103854), (-106.6478395, 23.74309074), (-106.6472376, 23.77618646), (-106.6472982, 23.80876072), (-106.647127, 23.84151129), (-106.6471277, 23.8741312), (-106.6473995, 23.90656505), (-106.6473138, 23.93916488), (-106.6473408, 23.97172031), (-106.6473796, 24.00435261), (-106.6472953, 24.0370137)]
    out_shp = r'X:\temp\polygon.shp'
    main(coords, out_shp)

Voici l'erreur que j'ai reçue:

runfile('X:/temp/corner_detection.py', wdir='X:/temp')
Traceback (most recent call last):

  File "<ipython-input-40-952256a990f1>", line 1, in <module>
    runfile('X:/temp/corner_detection.py', wdir='X:/temp')

  File "C:\Program Files (x86)\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 866, in runfile
    execfile(filename, namespace)

  File "C:\Program Files (x86)\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 87, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)

  File "X:/temp/corner_detection.py", line 48, in <module>
    main(coords, out_shp)

  File "X:/temp/corner_detection.py", line 43, in main
    write_shapefile(poly, out_shp)

  File "X:/temp/corner_detection.py", line 20, in write_shapefile
    layer = ds.CreateLayer('', None, ogr.wkbPolygon)

AttributeError: 'NoneType' object has no attribute 'CreateLayer'

Je suis convaincu que le problème a quelque chose à voir avec la lecture d'une shapelygéométrie. Comment puis-je créer un fichier de formes polygonales à partir d'une liste de coordonnées en utilisant uniquement les outils python gdal / ogr?

Vous obtenez cette erreur car le fichier ShapeFile est ouvert. Il ne peut donc pas le recréer.

Mais j'ai exécuté votre script et j'ai eu une erreur différente:

Traceback (most recent call last):
  File "test.py", line 48, in <module>
    main(coords, out_shp)
  File "test.py", line 43, in main
    write_shapefile(poly, out_shp)
  File "test.py", line 32, in write_shapefile
    geom = ogr.CreateGeometryFromWkb(poly.wkb)
AttributeError: 'str' object has no attribute 'wkb'

Vous semblez essayer d'extraire le binaire bien connu (WKB) du texte bien connu (WKT).

Vous obtenez déjà la représentation WKT de:

    return poly.ExportToWkt()

Donc, pour le corriger, utilisez simplement le créateur de géométrie Wkt au lieu de Wkb.

Changement:

geom = ogr.CreateGeometryFromWkb(poly.wkb)

À:

geom = ogr.CreateGeometryFromWkt(poly)

Votre source de données ( ds) n'est pas créée, peut-être parce que vous demandez un Esri Shapefileau lieu d'un ESRI Shapefile. En général, vous devez vérifier que vous récupérez une banque de données de pilotes *

driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(out_shp)