Create a shapefile from geometry with OGR

Question

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?

Answer 1

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
      }
    })

Answer 2

You can also create geometries with SQL by using SQLite dialect. You can execute SQL with Python but here is a one line example about how to do that with ogr2ogr. Because ogr2ogr requires two data stores, one for the source and another for the target some existing vector source must be used as a dummy input.

ogr2ogr -f "ESRI Shapefile" -a_srs epsg:4326 -dialect SQLite -sql "select ST_GeomFromText('LINESTRING (54 37,62 35,70 38,74 41.5)')" created_with_sql.shp foo.shp

Check the result with ogrinfo

ogrinfo created_with_sql.shp -al
INFO: Open of `created_with_sql.shp'
      using driver `ESRI Shapefile' successful.

Layer name: created_with_sql
Metadata:
  DBF_DATE_LAST_UPDATE=2021-03-30
Geometry: Line String
Feature Count: 1
Extent: (54.000000, 35.000000) - (74.000000, 41.500000)
Layer SRS WKT:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
FID: Integer64 (11.0)
OGRFeature(created_with_sql):0
  FID (Integer64) = 0
  LINESTRING (54 37,62 35,70 38,74.0 41.5)