Weird output when projecting/interpolating points on lines

Question

I am trying to project a points' shapefile on a multilinestrings' shapefile. (Project each point to the closest line)

Here's my code :

import geopandas as gpd
import time

start_time = time.time()

points4326 = gpd.GeoDataFrame.from_file('./Outputs/points.shp')
points = points4326.to_crs({'init': 'epsg:2269'})
lines = gpd.GeoDataFrame.from_file('./Shapefiles/lines.shp')
projected_points = gpd.GeoDataFrame(columns=['ID','Arrow','Comments','ZoneName','geometry'], crs={'init' :'epsg:2269'})

for index, row in points.iterrows():

    pt = gpd.GeoSeries(row['geometry'])
    # Project
    projOnLines = lines.project(pt)
    # print projOnLines

    #  Interpolate
    ptInterpolate = lines.interpolate(projOnLines)
    # print ptInterpolate
    projected_points.set_value(index, 'geometry', ptInterpolate.iloc[0])
    # print projected_points

    # Set the value of each column of projected_points with the value of the corresponding row in the input points' shapefile
    projected_points.set_value(index, 'ID', row['OBJECTID_l'])
    projected_points.set_value(index,'Arrow',row['Arrow'])
    projected_points.set_value(index,'Comments',row['Comments'])
    projected_points.set_value(index,'ZoneName',row['ZoneName'])

print projected_points.head(5)
projected_points.to_file('./Outputs/pointsOnLines.shp')

print "time elapsed: {:.2f}s".format(time.time() - start_time)

If I just put :

projOnLines = lines.project(pt)
ptInterpolate = lines.interpolate(projOnLines)

I get this error :

ctypes.ArgumentError: argument 3: <type 'exceptions.TypeError'>: wrong type

The problem comes from the fact that projOnLines's result is a Geoseries, and interpolate expects one float. Here's what print projOnLines shows :

0        199.894
1          False
2          False
3          False
4          False
      ...   
69955      False
69956      False
69957      False
dtype: object

Now If I take only that first value and pass it to interpolate like this :

ptInterpolate = lines.interpolate(projOnLines[0])

The script takes a long time to execute and when I finally get to display the result on QGIS I see all points are overlapped on one same location! (the attribute table has all features but they're somehow gathered in one place although I'm looping through all the points in the shapefile)

Any ideas?

Answer 1

I think you have a problem with your points shapefile and/or the projections (lines and points).

These commands are those of Shapely used by GeoPandas

From Coordinate of the closest point on a line

# Length along line that is closest to the point
print(line.project(point))
# Now combine with interpolated point on line
np = line.interpolate(line.project(point))

If I use your original script with my data

The results are the lengths along lines that are closest to the points (without False and with lines.interpolate(projOnLines[0]) or lines.interpolate(projOnLines))

for index, row in points.iterrows():
    pt =  row['geometry']
    projOnLines = lines.project(pt)
    print  index, projOnLines[0]

 0 7569.65191878
 1 1749.26685563
 2 18261.9544162
 ....
 9 18415.2725214

You can use directly Shapely and Fiona (also used by GeoPandas) to compare the results

from shapely.geometry import shape
import fiona
lines = shape(fiona.open('multiline.shp').next()['geometry'])
for pt in  fiona.open('points.shp'):
  pt= shape(pt['geometry'])
  print lines.project(pt)
7569.65191878
1749.26685563
18261.9544162
...
18415.2725214

In your case, GeoPandas recognizes and uses only a single point (result of ctypes.ArgumentError: argument 3: <type 'exceptions.TypeError'>: wrong type) and not

points are overlapped on one same location!

0        199.894
1          False 
....
69957      False
dtype: object

Answer 2

So this is what I ended up doing and it works, although the script takes forever to run and I have very large shapefiles, therefore the loop over both the points and lines shapefiles is kinda ugly.. Would love to know if there is any idea to improve/optimize this code without having to loop over both shapefiles.

import geopandas as gpd
import time
start_time = time.time()

points4326 = gpd.GeoDataFrame.from_file('./points_test_zone.shp')
points = points4326.to_crs({'init': 'epsg:2269'})

lines4326 = gpd.GeoDataFrame.from_file('./lines_test_zone.shp')
lines =lines4326.to_crs({'init': 'epsg:2269'})

projected_points = gpd.GeoDataFrame(columns=['ID','Arrow','Comments','ZoneName','geometry'], crs={'init' :'epsg:2269'})

for indexp, rowp in points.iterrows():

    pt = gpd.GeoSeries(rowp['geometry'])

    listDist = []

    for indexl, rowl in lines.iterrows():

        ln = gpd.GeoSeries(rowl['geometry'])

        smallDic = {"lineID": indexl, "dis": ln.distance(pt)[0]}

        listDist.append(smallDic)

    lineID = sorted(listDist, key=lambda k: k['dis'])[0].values()[0]

for indexl, rowl in lines.iterrows(): 

    if indexl == lineID:

        selecLine = gpd.GeoSeries(rowl['geometry'])

        prjPt = selecLine.interpolate(selecLine.project(pt)[0]).iloc[0]

        projected_points.set_value(indexp, 'geometry', prjPt)

        projected_points.set_value(indexp, 'ID', rowp['OBJECTID_l'])

        projected_points.set_value(indexp,'Arrow',rowp['Arrow'])

        projected_points.set_value(indexp,'Comments',rowp['Comments'])

        projected_points.set_value(indexp,'ZoneName',rowp['ZoneName'])

projected_points.to_file('./pointsOnLines.shp')

print "time elapsed: {:.2f}s".format(time.time() - start_time)