4

I'm trying to snap two Linestrings shapefiles using Shapely/Geopandas. I tried :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1.to_crs({'init': 'epsg:2227'})
lines2 =gpd.GeoDataFrame.from_file('lines2.shp')
lines2.to_crs({'init': 'epsg:2227'})
res = snap(lines1, lines2, 10)
res.to_file('result.shp', driver="ESRI Shapefile")

But I got this error :

Traceback (most recent call last):
File "C:/Users/Documents/script.py", line 14, in <module>
    res = snap(lines1, lines2, 10)
File "C:\Program Files (x86)\Python27\lib\site-packages\shapely\ops.py", line 303, in snap
    return(geom_factory(lgeos.methods['snap'](g1._geom, g2._geom, tolerance)))
File "C:\Program Files (x86)\Python27\lib\site-packages\pandas\core\generic.py", line 2672, in __getattr__
    return object.__getattribute__(self, name)
AttributeError: 'GeoDataFrame' object has no attribute '_geom'

Process finished with exit code 1.

Any idea on what the problem is?

EDIT :

Updated code based on @caiohamamura's answer :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1 = lines1.to_crs({'init': 'epsg:2227'})
lines2 = gpd.GeoDataFrame.from_file('lines2.shp')
lines2 = lines2.to_crs({'init': 'epsg:2227'})
res = lines1
lines2_union = lines2.geometry.unary_union
res.geometry = res.geometry.apply(lambda x: snap(x, lines2_union, 10))
res.to_file('result.shp', driver="ESRI Shapefile")

Now I have this error :

Process finished with exit code -1073741819 (0xC0000005)

EDIT 2 :

I no longer get the above error. The script runs well, although the snapping doesn't seem to work that well. Here's an example of the result : Snapping two lines

The blue line is the result of the snapping

The black one is lines2

The red one is lines1

I tried changing the tolerance but it seems that if I give a value that's > to 15 feets I might get a wrong output. The lines do not get snapped totally to each other. In some places they do, in others they don't as demonstrated on the screenshot above.

  • This should work, I've executed that for a sample data I have here and it works fine. At what specific command do you get the error? What coordinate system are you using? If it is lat/long 10 will actually mean(1110 km). – caiohamamura Jul 21 '16 at 0:54
  • @caiohamamura I don't know at which specific command the error pops up now. This is the only shown message I get : Process finished with exit code -1073741819 (0xC0000005) As for the coordinates system I'm not using WGS84, that's why I'm doing a transformation at the beginning to 2227 (projected coordinate system of SF) – GeoSal Jul 21 '16 at 1:05
  • I think this error code is due to the size of your shapefile. It stands for the classical stackoverflow error. So you need to increase stack size. Read this. Also you can run the script in interactive terminal line by line so you can tell which command will throw the error. But I would go for a different tool as geopandas is not that mature. For this I would use PostGIS or SpatialLite. – caiohamamura Jul 21 '16 at 1:20
  • @caiohamamura I'm looking into the link u sent. I think that the problem is a bug in the script though. What unit does shapely use for the tolerance of the snapping? Since I used a projected coordinates system, I assumed it's "feet". Is it? – GeoSal Jul 21 '16 at 1:28
  • epsg:2227 seems to be in feets. But for the exit code you gave I'm pretty sure the error is due to exceeding limits for processing, I would go for another tool. As I said, geopandas is a new library. You could also try to execute the snapping for just the first row just to test. snap(lines1.geometry[0], lines2_union, 10). – caiohamamura Jul 21 '16 at 1:39
2

For this to work you will have to first unary_union lines2.geometry so that lines1 will snap to any line in lines2. And as snap isn't a vectorized function you will have to execute for each geometry using apply:

res = lines1
lines2_union = lines2.geometry.unary_union
res.geometry = lines1.geometry.apply(lambda x: snap(x, lines2_union, 10))
  • Thank u for ur quick answers! I tried what u suggested but i got nearly the same error as above except for the last line : AttributeError: 'GeoSeries' object has no attribute '_geom' Process finished with exit code 1 – GeoSal Jul 20 '16 at 18:20
  • Just saw ur edit! Will try what u suggested and get back to u – GeoSal Jul 20 '16 at 18:23
  • Didn't work :/ I'm still having that error AttributeError: 'GeoSeries' object has no attribute '_geom' – GeoSal Jul 20 '16 at 20:51
  • The error is in this line : res['geometry'] = snap(lines1['geometry'], lines2['geometry'], 10) – GeoSal Jul 20 '16 at 21:02
  • Yeah, that's because snap isn't vectorized, so you have to perform snap for each element of lines1 using apply as I edited. – caiohamamura Jul 21 '16 at 0:20

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.