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'])

print projected_points.head(5)

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?


2 Answers 2


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
# Now combine with interpolated point on line
np = line.interpolate(line.project(point))

If I use your original script with my data

enter image description here

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

enter image description here

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)

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
  • Thanks @gene for your answer. I tested the script on a small chunck and it didn't work either. I don't think there is a problem with the shapefile. I've made some changes in the code, please refer to my answer below. I basically had to loop over the points and lines' shapefiles. Not that good but it still works for me :/
    – GeoSal
    Aug 26, 2016 at 21:26

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]}


    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'])





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

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.