2

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?

4

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
      }
    })
2
  • Thanks Thomas for this detailed explanation. It was indeed the best answer to my question. GDAL is great. I knew that "fiona" is a higher level package that helps in exploiting the underlying functions of GDAL/OGR. Until now I resisted to use it before developing a good understanding of GDAL/OGR. At some point I have to use it though :) – Sinooshka Mar 30 at 12:25
  • You can accept the answer in this case e.g stackoverflow.com/help/someone-answers – ThomasG77 Mar 30 at 12:27
4

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)
2
  • Thanks alot for this insight. You are right. The SQL dialect is also a good way. At this time however I was intending to limit myself to use the functions that I can access directly from "ogr" and see if I can do it without resorting to other dialects. – Sinooshka Mar 30 at 12:33
  • Tried echo '{"type": "FeatureCollection", "features": [{"type": "Feature","properties": {"id": 1}, "geometry": {"type": "LineString","coordinates": [[54,37], [62,35], [70,38], [74,41.5]]}}]}' | ogr2ogr out.shp /vsistdin/ also – ThomasG77 Mar 30 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.