1

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

4
  • Does it need the 1st point as the 5th point to close the polygon?
    – Banger
    Jan 10, 2015 at 5:40
  • Nope. A, B, C, D (anti-clockwise) are the vertices of the bounding box of the polygon shapefile. I am trying to draw the A-D line segment and export it as a shapefile. To add some context: I will draw lines parallel to this A-D line along the entire bounding box
    – Akhil
    Jan 10, 2015 at 5:42
  • Do you need to close the shape file at the end, I'm not really familiar with Pyqgis.
    – Banger
    Jan 10, 2015 at 6:12
  • Here is some additional information that even I am looking through.. You could help out if you are familiar with Python atleast. Look through the question hyperlinked in the above post and also Vector Layer information from PyQGIS Cookbook
    – Akhil
    Jan 10, 2015 at 6:29

1 Answer 1

6

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 ] ] }
2
  • Just a clarification... If I had to add multiple lines, I would have to run a loop. Would a for loop from before Add.Point() and ending after feat.SetGeometry() suffice? Considering I have to draw equally spaced parallel lines, increasing the x coordinate by a fixed distance.
    – Akhil
    Jan 13, 2015 at 11:40
  • 1
    before linelyr = ogr.Geometry(ogr.wkbLineString)and after layr1.CreateFeature(feat) otherwise no line is created
    – gene
    Jan 13, 2015 at 20:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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