OGR – Writing Shapefile from Shapely Geometry with No Feature Added Error

ogrpyqgisqgisshapely

I am trying to write a shapefile from a Shapely geometry using PyQGIS in the QGIS 2.6.0 console.

Lines of code before the OGR method to write shapefile:

extnt = poly.extent() # gives a QtRectangle output -> extnt
x_min = extnt.xMinimum()
x_max = extnt.xMaximum()
y_min = extnt.yMinimum()
y_max = extnt.yMaximum()

A_bbox = QgsPoint(x_min, y_min)
B_bbox = QgsPoint(x_max, y_min)
C_bbox = QgsPoint(x_max, y_max)
D_bbox = QgsPoint(x_min, y_max)

where poly is a polygon shapefile.

Source for the Code below :- How to write Shapely geometries to shapefiles? and GitHub post

linelyr = ogr.Geometry(ogr.wkbLineString)
linelyr.AddPoint(x_min, y_min)
linelyr.AddPoint(x_min, y_max)
out_line = ogr.Geometry(ogr.wkbLineString)
out_line.AddGeometry(linelyr)
out_line.ExportToWkb

driver = ogr.GetDriverByName("Esri Shapefile")
ds = driver.CreateDataSource("outputlocation.shp")
layr1 = ds.CreateLayer('', None, ogr.wkbLineString)

layr1.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layr1.GetLayerDefn()

feat = ogr.Feature(defn)
feat.SetField('id', 1)

geom = ogr.CreateGeometryFromWkb(out_line.wkb) # or maybe just 'outline'
feat.SetGeometry(geom)

layr1.CreateFeature(feat)
ds = layr1 = feat = geom = None

The shapefile is created, along with all the necessary files, meaning no errors. But the shapefile, when opened in QGIS has no features. Just a constructed attribute table.

Pls assist…

P. S. Any Other solution besides OGR usage also works, as long as it can be coded in the Python console of QGIS

Best Answer

Some remarks because you mix many things ( if you want yo use your script outside the console of QGIS, install the Python module GDAL (osgeo) in your Python installation: it is installed in the Python version of QGIS if you are on Windows).

1) There are no Shapely geometries in your script, only

  • PyQGIS geometries QgsPoint(x_min, y_min), ...
  • ogr geometries linelyr = ogr.Geometry(ogr.wkbLineString), ...

2) out_line is unnecessary because linelyr is already an ogr LineString (and your creation of out_line gives invalid geometries (empty LineStrings = LineString(LineString))

linelyr = ogr.Geometry(ogr.wkbLineString)
linelyr.AddPoint(5, 47)
linelyr.AddPoint(5, 55)
print linelyr.ExportToJson()
{ "type": "LineString", "coordinates": [ [ 5.0, 47.0, 0.0 ], [ 5.0, 55.0, 0.0 ] ] }
# and
out_line = ogr.Geometry(ogr.wkbLineString)
out_line.AddGeometry(linelyr)
print out_line.ExportToWkt()
LINESTRING EMPTY

3) In the same way, you don't need geom = ogr.CreateGeometryFromWkb(out_line.wkb) because as said, out_line or linelyr are already ogr LineStrings.

3) At the end of the process, with ogr, if you not close the resulting shapefile, it remains empty.

So:

from osgeo import ogr
# create ogr geometry
linelyr = ogr.Geometry(ogr.wkbLineString)
linelyr.AddPoint(5, 47)
linelyr.AddPoint(5, 55)
# create the shapefile
driver = ogr.GetDriverByName("Esri Shapefile")
ds = driver.CreateDataSource("outputlocation.shp")
layr1 = ds.CreateLayer('',None, ogr.wkbLineString)
# create the field
layr1.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
# Create the feature and set values
defn = layr1.GetLayerDefn()
feat = ogr.Feature(defn)
feat.SetField('id', 1)
feat.SetGeometry(linelyr)
layr1.CreateFeature(feat)
# close the shapefile
ds.Destroy()

enter image description here

If you want to use Shapely:

from shapely.geometry import Point, LineString
line = LineString([Point(5,47),Point(5,55)])
# conversion to ogr geometry
linelyr = ogr.CreateGeometryFromWkt(line.wkt)
print linelyr.ExportToJson()
{ "type": "LineString", "coordinates": [ [ 5.0, 47.0, 0.0 ], [ 5.0, 55.0, 0.0 ] ] }

If you want to use PyQGIS:

line = QgsGeometry.fromPolyline([QgsPoint(5,47),QgsPoint(5,55)]))
# conversion to ogr geometry
linelyr = ogr.CreateGeometryFromWkt(line.exportToWkt())
print linelyr.ExportToJson()
{ "type": "LineString", "coordinates": [ [ 5.0, 47.0 ], [ 5.0, 55.0 ] ] }
Related Question