python – How to Create a Shapefile from Geometry with OGR in Python

geometryogrpythonshapefilevector-layer

As I see from the GDAL/OGR available solutions, most of the new OGR data sources (shapefile etc.) are being created from an already existing similar file. This way many steps of directly creating Layers, Fields and Attributes are avoided and thus makes life much easier to manipulate the existing structure of a vector layer.

I wonder however how to :

  1. Create a simple geometry (a line)
  2. Create a shapefile from scratch
  3. Put the geometry inside the shapefile and save it to the disc

Here is the code that I have so-far prepared.

import sys 
from osgeo import ogr
from osgeo import osr

# Creating a line geometry
linegeo = ogr.Geometry(ogr.wkbLineString)
linegeo.AddPoint(54,37)
linegeo.AddPoint(62,35)
linegeo.AddPoint(70,38)
linegeo.AddPoint(74,41.5)

  
# Set up the shapefile driver 
driver = ogr.GetDriverByName("ESRI Shapefile")

# create the data source
ds = driver.CreateDataSource("line.shp")

# create the spatial reference system, WGS84
srs =  osr.SpatialReference()
srs.ImportFromEPSG(4326)

# create one layer 
layer = ds.CreateLayer("line", srs, ogr.wkbLineString)

At this point I think I need to somehow use the "ogr.FieldDefn()" in order to place the above created line geometry inside the newly created layer.

How do I do that?

Best Answer

You need to create a feature an add the feature (a feature = a geometry + fields) to the layer

import sys 
from osgeo import ogr
from osgeo import osr

# Creating a line geometry
linegeo = ogr.Geometry(ogr.wkbLineString)
linegeo.AddPoint(54,37)
linegeo.AddPoint(62,35)
linegeo.AddPoint(70,38)
linegeo.AddPoint(74,41.5)

  
# Set up the shapefile driver 
driver = ogr.GetDriverByName("ESRI Shapefile")

# create the data source
ds = driver.CreateDataSource("line.shp")

# create the spatial reference system, WGS84
srs =  osr.SpatialReference()
srs.ImportFromEPSG(4326)

# create one layer 
layer = ds.CreateLayer("line", srs, ogr.wkbLineString)

# Add an ID field
idField = ogr.FieldDefn("id", ogr.OFTInteger)
layer.CreateField(idField)

# Create the feature and set values
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(linegeo)
feature.SetField("id", 1)
layer.CreateFeature(feature)

feature = None

# Save and close DataSource
ds = None

Code adapted from https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-layer-from-the-extent-of-an-existing-layer

You may prefer using Fiona, a library combining shapely with ogr module from GDAL as it can avoid boilerlate at GDAL level. In this case, the recipe could be

import fiona
from fiona.crs import from_epsg
from shapely.geometry import LineString, mapping

outname = 'line1.shp'
schema = {
    'geometry': 'LineString', 
    'properties': {
        'id': 'int'
    }
}

line = [
    [54,37],
    [62,35],
    [70,38],
    [74,41.5]
]

# write your shapefile, as projection of epsg: 4326
with fiona.open(outname, 'w', 'ESRI Shapefile', schema=schema, crs = from_epsg(4326)) as out:
    linestring = LineString(line)
    out.write({
      'geometry': mapping(linestring),
      'properties': {
          'id': 1
      }
    })
Related Question