Sunday, 5 January 2020

python - geopandas.intersects() over each row


In trying to answer Writing selected features to shapefile using OGR and Python? I ran into problems with geopandas.GeoSeries intersects methods.


I have the following line and polygon:


import pandas as pd
import geopandas
import pickle


line = pickle.loads(
b'\x80\x03cgeopandas.geodataframe\nGeoDataFrame\nq\x00)\x81q\x01}q\x02(X\x05\x00\x00\x00_dataq\x03cpandas.core.internals.managers\nBlockManager\nq\x04)\x81q\x05(]q\x06(cpandas.core.indexes.base\n_new_Index\nq\x07cpandas.core.indexes.base\nIndex\nq\x08}q\t(X\x04\x00\x00\x00dataq\ncnumpy.core.multiarray\n_reconstruct\nq\x0bcnumpy\nndarray\nq\x0cK\x00\x85q\rC\x01bq\x0e\x87q\x0fRq\x10(K\x01K\x02\x85q\x11cnumpy\ndtype\nq\x12X\x02\x00\x00\x00O8q\x13K\x00K\x01\x87q\x14Rq\x15(K\x03X\x01\x00\x00\x00|q\x16NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK?tq\x17b\x89]q\x18(X\x03\x00\x00\x00FIDq\x19X\x08\x00\x00\x00geometryq\x1aetq\x1bbX\x04\x00\x00\x00nameq\x1cNu\x86q\x1dRq\x1eh\x07cpandas.core.indexes.range\nRangeIndex\nq\x1f}q (h\x1cNX\x05\x00\x00\x00startq!K\x00X\x04\x00\x00\x00stopq"K\x01X\x04\x00\x00\x00stepq#K\x01u\x86q$Rq%e]q&(h\x0bh\x0cK\x00\x85q\'h\x0e\x87q(Rq)(K\x01K\x01K\x01\x86q*h\x12X\x02\x00\x00\x00i8q+K\x00K\x01\x87q,Rq-(K\x03X\x01\x00\x00\x00Rq?(K\x01K\x01\x85q@h\x15\x89]qAh\x19atqBbh\x1cNu\x86qCRqDh\x07h\x08}qE(h\nh\x0bh\x0cK\x00\x85qFh\x0e\x87qGRqH(K\x01K\x01\x85qIh\x15\x89]qJh\x1aatqKbh\x1cNu\x86qLRqMe}qNX\x06\x00\x00\x000.14.1qO}qP(X\x04\x00\x00\x00axesqQh\x06X\x06\x00\x00\x00blocksqR]qS(}qT(X\x06\x00\x00\x00valuesqUh)X\x08\x00\x00\x00mgr_locsqVcbuiltins\nslice\nqWK\x00K\x01K\x01\x87qXRqYu}qZ(hUh4hVhWK\x01K\x02K\x01\x87q[Rq\\ueustq]bX\x04\x00\x00\x00_typq^X\t\x00\x00\x00dataframeq_X\t\x00\x00\x00_metadataq`]qa(X\x03\x00\x00\x00crsqbX\x15\x00\x00\x00_geometry_column_nameqcehb}qdX\x04\x00\x00\x00initqeX\t\x00\x00\x00epsg:4326qfshch\x1aub.'
)

poly = pickle.loads(
b'\x80\x03cgeopandas.geodataframe\nGeoDataFrame\nq\x00)\x81q\x01}q\x02(X\x05\x00\x00\x00_dataq\x03cpandas.core.internals.managers\nBlockManager\nq\x04)\x81q\x05(]q\x06(cpandas.core.indexes.base\n_new_Index\nq\x07cpandas.core.indexes.base\nIndex\nq\x08}q\t(X\x04\x00\x00\x00dataq\ncnumpy.core.multiarray\n_reconstruct\nq\x0bcnumpy\nndarray\nq\x0cK\x00\x85q\rC\x01bq\x0e\x87q\x0fRq\x10(K\x01K\x02\x85q\x11cnumpy\ndtype\nq\x12X\x02\x00\x00\x00O8q\x13K\x00K\x01\x87q\x14Rq\x15(K\x03X\x01\x00\x00\x00|q\x16NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK?tq\x17b\x89]q\x18(X\x03\x00\x00\x00FIDq\x19X\x08\x00\x00\x00geometryq\x1aetq\x1bbX\x04\x00\x00\x00nameq\x1cNu\x86q\x1dRq\x1eh\x07cpandas.core.indexes.range\nRangeIndex\nq\x1f}q (h\x1cNX\x05\x00\x00\x00startq!K\x00X\x04\x00\x00\x00stopq"K\x06X\x04\x00\x00\x00stepq#K\x01u\x86q$Rq%e]q&(h\x0bh\x0cK\x00\x85q\'h\x0e\x87q(Rq)(K\x01K\x01K\x06\x86q*h\x12X\x02\x00\x00\x00i8q+K\x00K\x01\x87q,Rq-(K\x03X\x01\x00\x00\x00CM\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\xfa\xb9b\x96#?\xe7\xbf`\xf3T\n\xe9k\xce\xbf\xd4\xb8\x81Qg\x1a\xd7\xbf\xde\x822\xac\xd6[\xd0?\xd4\xb8\x81Qg\x1a\xd7\xbf@\xfa\x9a\xa7RH\xcf\xbf\xfa\xb9b\x96#?\xe7\xbf`\xf3T\n\xe9k\xce\xbfq?bh7)Rq@CM\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\x88\xb6\xbf\xc7\xee\xd0\xd6\xbf\xf8\xf7\xd8\x1d\xda\xfe\xce\xbf\x88\xb6\xbf\xc7\xee\xd0\xd6\xbf*\x85\xf45O\xa5\xd0?\xa0\xe5\xc8\xcf\x15\xb3\xac?4\x8d\x1b\x18u\xa6\xd1?\x88\xb6\xbf\xc7\xee\xd0\xd6\xbf\xf8\xf7\xd8\x1d\xda\xfe\xce\xbfqAbh7)RqBCM\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00<\xb4\xfd=v\x87\xd6\xbf\xa8\xf5\x16\x94a\xb5\xce\xbf\x00\xca\xb0ZoA\xa9?\x04\x84\x13\xf1\x92\x80\xd0?@\xa5\x90\xbe\xe6\xa9\xa4?@\xfa\x9a\xa7RH\xcf\xbf<\xb4\xfd=v\x87\xd6\xbf\xa8\xf5\x16\x94a\xb5\xce\xbfqCbetqDbe]qE(h\x07h\x08}qF(h\nh\x0bh\x0cK\x00\x85qGh\x0e\x87qHRqI(K\x01K\x01\x85qJh\x15\x89]qKh\x19atqLbh\x1cNu\x86qMRqNh\x07h\x08}qO(h\nh\x0bh\x0cK\x00\x85qPh\x0e\x87qQRqR(K\x01K\x01\x85qSh\x15\x89]qTh\x1aatqUbh\x1cNu\x86qVRqWe}qXX\x06\x00\x00\x000.14.1qY}qZ(X\x04\x00\x00\x00axesq[h\x06X\x06\x00\x00\x00blocksq\\]q](}q^(X\x06\x00\x00\x00valuesq_h)X\x08\x00\x00\x00mgr_locsq`cbuiltins\nslice\nqaK\x00K\x01K\x01\x87qbRqcu}qd(h_h4h`haK\x01K\x02K\x01\x87qeRqfueustqgbX\x04\x00\x00\x00_typqhX\t\x00\x00\x00dataframeqiX\t\x00\x00\x00_metadataqj]qk(X\x03\x00\x00\x00crsqlX\x15\x00\x00\x00_geometry_column_nameqmehl}qnX\x04\x00\x00\x00initqoX\t\x00\x00\x00epsg:4326qpshmh\x1aub.'
)

Which when plotted together look like this:


poly.plot(ax=line.plot(), color='red')


Line and polygons intersecting


I want to calculate and return which polygons in the triangular strip intersects with the line. I tried both of the following:


poly.intersects(line)
line.intersects(poly)

Both return True for the first item but False for the remaining 5 polygons:


0     True
1 False
2 False

3 False
4 False
5 False
dtype: bool

I had expected that it would return True for polygons 0, 1, and 2 and False for the rest but clearly this is not the case. From the GeoSeries intersects documentation I thought it would return a result based on each item in the Series.


I tried to instead iterate over each row:


for row in poly.iterrows():
row = geopandas.GeoSeries(row)
line.intersects(row)


But this returns the following error:


AttributeError: 'int' object has no attribute 'is_empty'

How should I calculate which polygons the line intersects with, so that my returned result is something like the following?


0     True
1 True
2 True
3 False
4 False

5 False

Answer



Not sure but I think you are running out of lines after the first comparision of polygon1-line1, second comparison will be polygon2-nothing.


If you for each polygon compare all lines (in your case one) you will get the results you want:


poly.geometry.map(lambda x: x.intersects(line.geometry.any()))

0 True
1 True
2 True
3 False

4 False
5 False

It is probably better to use Spatial Join.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...