3

I have a PostGIS layer with linestring and I am trying to write the layer to shapefile using fiona and shapely.

Here is what I have so far:

from shapely.geometry import mapping,LineString
from shapely.wkt import loads,load
import fiona
from fiona.crs import from_epsg
import psycopg2
conn = psycopg2.connect("dbname='postgres',etc..) #connecting to DB
cur = conn.cursor()  #setting up connection cursor
crs = from_epsg(4269)
schema = {'geometry': 'LineString', 'properties': {'State': 'str', 'Date': 'str'}}
with fiona.open(r"routes.shp", "w", "ESRI Shapefile", schema, crs) as output:
    cur.execute("""select state,date,ST_AsText(geom) from trips;""")
    rows = cur.fetchall()
    conn.commit()
    for row in rows:
        print row
        g = loads(row[2])
        output.write({'properties':{'State': row[0], 'Date': row[1]},
        'geometry': mapping(g)})

It throws an error when it gets to:

('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)')
Traceback (most recent call last):
  File "C:\Users\Ralph\Documents\PythonScripts\LocationHistory.py", line 57, in <module>
    g = loads(row[2])
  File "C:\Python27\ArcGIS10.1\lib\site-packages\shapely\wkt.py", line 10, in loads
    return geom_from_wkt(data) #factory(geom)
  File "C:\Python27\ArcGIS10.1\lib\site-packages\shapely\geometry\base.py", line 55, in geom_from_wkt
    "Could not create geometry because of errors while reading input."
ReadingError: Could not create geometry because of errors while reading input.

I used the same select query in postgres on the routes table and

('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)') 

seems to be the error.

It writes the shapefile up to this record. what is strange is if I export the table to a shapefile using the PostGIS shapefile and dbf loader it exports the entire layer to a shapefile correctly.

My question is twofold:

How did postgis make a linestring from only one coordinate? How would I go about circumventing this problem?

enter image description here

shapefile from PostGIS - this record ('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)') which throws the error in shapely when trying to create a linestring (because it is only 1 coordinate) shows up in the PostGIS shapefile as a record with no geometric attribute

enter image description here

  • I haven't tested it, but your Linestring consists of only one point (one coordinate), I think it should have at least two or it may be considered invalid – chrki Nov 19 '16 at 0:56
  • yes that has to be it, but how did it create a linestring in postgis with only 1 lat and 1 long value and successfully output it to a shapefile? – ziggy Nov 19 '16 at 1:02
  • Please edit the question to identify the contents of the row was represented. There are any number of ways to write invalid output from invalid input; writing valid output is somewhat trickier, and usually involves data loss. – Vince Nov 19 '16 at 12:17
  • @Vince i have a picture of the row it throws the error at...('South Dakota', '2016-06-19', 'LINESTRING(-103.8036896 44.4258771)') – ziggy Nov 19 '16 at 18:49
  • No, the row in the shapefile corresponding to the invalid geometry. Is it also degenerate? – Vince Nov 19 '16 at 20:55
3

The problem is that a LineString with only one point is not valid with GEOS (Shapely and many others)

With Shapely:

The LineString constructor takes an ordered sequence of 2 or more (x, y[, z]) point tuples. (from Shapely LineString)

With PostGIS:

A linestring is a path between locations. It takes the form of an ordered series of two or more points (from PostGIS intro: geometries).

Therefore a LineString with one Point gives an error:

from shapely.geometry import LineString
line = LineString([(-103.8036896,44.4258771)])
Traceback (most recent call last):
....
ValueError: LineStrings must have at least 2 coordinate tuples

You can always save your original geometry as a shapefile

import fiona
line = fiona.open("line.shp")    
line.crs
{'init': u'epsg:4269'}
first =line.next()
print first
{'geometry': {'type': 'LineString', 'coordinates': [(-103.8036896, 44.4258771)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'FID', 0.0)])}

But you don't see anything in a GIS because

from shapely.geometry import shape
shape(first['geometry'])
....
ValueError: LineStrings must have at least 2 coordinate tuples

If you consider a Point as a Vector with an origin (point (0,0)) and an end (coordinates x, y of the point)

enter image description here

You can use

line = LineString([(0, 0), (-103.8036896,44.4258771)])

(with the node of the original shapefile in blue)

enter image description here

But don't forget that if you work with EPSG 4269 coordinates, the unit is degree (you don't work in a Cartesian plane)

Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane.

  • okay this is good info, I am thinking about intersecting each linestring with a polygon that i know it is contained in to weed out that bad linestring – ziggy Nov 19 '16 at 17:15

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.