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

Akhil
  • 1,435
  • 3
  • 19
  • 24
  • Does it need the 1st point as the 5th point to close the polygon? – Banger Jan 10 '15 at 05: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 '15 at 05:42
  • Do you need to close the shape file at the end, I'm not really familiar with Pyqgis. – Banger Jan 10 '15 at 06: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 '15 at 06:29

1 Answers1

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 ] ] }
gene
  • 54,868
  • 3
  • 110
  • 187
  • 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 '15 at 11:40
  • 1
    before linelyr = ogr.Geometry(ogr.wkbLineString)and after layr1.CreateFeature(feat) otherwise no line is created – gene Jan 13 '15 at 20:56