0

I'd tried replicating intersecting two shapefiles from Python or command line by anna-magdalena which requires intersecting two shapes. I used @gene answer for that:

import os
from shapely.geometry import mapping, shape
import fiona

os.chdir('/media/dir')
# Create schema of values that new shape will contain
schema =  {'geometry': 'Polygon','properties': {'area': 'float:13.3', 'id1': 'str', 'id3': 'str'}}

with fiona.open('13fiona.shp', 'w',driver='ESRI Shapefile', schema=schema) as output:
    for c1 in fiona.open('1.shp'):
        for c2 in fiona.open('3.shp'):
           if shape(c1['geometry']c1).intersects(shape(c2['geometry'])):     
              area = shape(c1['geometry']).intersection(shape(c2['geometry'])).area
              prop = {'area': area, 'id1' : c1['id'],'id3': c2['id']} 
              geom = mapping(shape(c1['geometry']).intersection(shape(c2['geometry'])))
              output.write({'geometry':geom,'properties': prop})

And I got this error:

ValueError: Record's geometry type does not match collection schema's geometry type: 'GeometryCollection' != 'Polygon'

After check the structure of objects c1, c2 and geom I got this:

Structure of c1 polygon:

{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.71213794799996, 1.2098971930000744), (-74.71 ...

and c2 polygon:

{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.69988807399994, 1.2222...

I found that intersection polygon assigned to geom objetc have this structure:

{'type': 'GeometryCollection', 'geometries': [{'type': 'Point', 'coordinates': (-74.71039521399996, 1.2098918560000698)}, {'type': 'Point', 'coordinates': (-74.70929715499994, 1.2102729000000636)}, ...

The c1 and c2 shapes cand be found here I'm using shapely 1.6b2 version

Any advice?

0

1 Answer 1

1

The problem is that the intersection of two polygons is not always a polygon (look at Intersecting two shape problem using geopandas)

For example, with one polygon in each shapefile

c = fiona.open('1.shp')
d = fiona.open('3.shp')
u = shape(c[0]['geometry']) # first polygon  of 1.shp
v  = shape(d[3728]['geometry']) # polygon 3728 of 3.shp
print(u.intersection(v)) 
GEOMETRYCOLLECTION (POINT (-74.71039521399996 1.20989185600007)....

The intersection is not a single polygon but a GEOMETRYCOLLECTION with 6 points, 21 lines and 10 polygons

for i in  u.intersection(v):
    print(i.wkt)
POINT (-74.71039521399996 1.20989185600007)
POINT (-74.70929715499994 1.210272900000064)
POINT (-74.70912409499994 1.213558722000073)
POINT (-74.70873115599994 1.210839667000073)
POINT (-74.70725587999993 1.212713845000053)
POINT (-74.70721854399994 1.212594696000053)
LINESTRING (-74.70901529499997 1.214883536000059, -74.70901765999997 1.214711119000071)
LINESTRING (-74.70901656099994 1.214098235000051, -74.70901886199994 1.214013153000053)
LINESTRING (-74.70923533299998 1.213327715000048, -74.70925547599995 1.213245727000071)
LINESTRING (-74.70925547599995 1.213245727000071, -74.70926798099998 1.21317024800004)
LINESTRING (-74.70927560399997 1.213063977000047, -74.70927258999996 1.213044866000075)
LINESTRING (-74.70899918999999 1.212640957000076, -74.70893900099998 1.212590409000029)
LINESTRING (-74.70893900099998 1.212590409000029, -74.70887306699996 1.212557479000054)
LINESTRING (-74.70873137199999 1.21253157700005, -74.70865583399996 1.212528399000064)
LINESTRING (-74.70843252699996 1.212550415000067, -74.70836507699994 1.212569773000041)
LINESTRING (-74.70836507699994 1.212569773000041, -74.70830443299997 1.212597207000044)
LINESTRING (-74.70830443299997 1.212597207000044, -74.70825156899997 1.212633270000026)
LINESTRING (-74.70809120499996 1.212741707000021, -74.70802753899994 1.212767495000037)
LINESTRING (-74.70802753899994 1.212767495000037, -74.70789617999998 1.212813200000028)
LINESTRING (-74.70777746899995 1.212865851000061, -74.70772889399996 1.212898795000058)
LINESTRING (-74.70759040599995 1.212969122000061, -74.70754020999993 1.212967823000042)
LINESTRING (-74.70885317899996 1.210700220000035, -74.70890060899995 1.210660536000034)
LINESTRING (-74.70898691799994 1.210572569000021, -74.70902578099998 1.210524269000075)
LINESTRING (-74.70902578099998 1.210524269000075, -74.70906696999998 1.210478309000052)
LINESTRING (-74.70906696999998 1.210478309000052, -74.70911341099998 1.210437633000026)
LINESTRING (-74.70939977699999 1.210201349000045, -74.70946099199995 1.210175537000055)
LINESTRING (-74.71006820699995 1.209912458000076, -74.71014630399998 1.209903631000032)
POLYGON ((-74.70901384319636 1.214919839732185, -74.70901529499997 1.214883536000059, -74.70901384399997 1.21491984000005, -74.70901384319636 1.214919839732185))
POLYGON ((-74.70902363705481 1.21392356305488, -74.70903007099997 1.213849638000056, -74.70904105599999 1.213773292000042, -74.70903007199996 1.213849638000056, -74.70902363705481 1.21392356305488))
POLYGON ((-74.70926682999993 1.212989782000022, -74.70924596999998 1.212939768000069, -74.70921041999998 1.212884735000046, -74.70924597099997 1.212939768000069, -74.70926682999993 1.212989782000022))
POLYGON ((-74.70851186799996 1.212536940481429, -74.70850468199995 1.212537629000053, -74.70843252699996 1.212550415000067, -74.70850468099997 1.212537629000053, -74.70851186799996 1.212536940481429))
POLYGON ((-74.70825156899997 1.212633270000026, -74.70820193899993 1.212673080000059, -74.70819125845584 1.212680792666712, -74.70820193699996 1.212673081000048, -74.70825156899997 1.212633270000026))
POLYGON ((-74.70772889399996 1.212898795000058, -74.70768445199997 1.212931320000052, -74.70763935899998 1.212956696000049, -74.70759040599995 1.212969122000061, -74.70763935799994 1.212956696000049, -74.70768445099998 1.212931320000052, -74.70772889399996 1.212898795000058))
POLYGON ((-74.70722619012928 1.212430557464625, -74.7072293166775 1.212416530935534, -74.70722931699999 1.212416531000031, -74.70722619012928 1.212430557464625))
POLYGON ((-74.70890060899995 1.210660536000034, -74.70894613799999 1.210618940000074, -74.70898691799994 1.210572569000021, -74.70894613799999 1.210618941000064, -74.70890060899995 1.210660536000034))
POLYGON ((-74.70992560399998 1.209943787000043, -74.70999416199999 1.209925364000071, -74.71006820699995 1.209912458000076, -74.70999416299998 1.209925364000071, -74.70992560399998 1.209943787000043))
POLYGON ((-74.71056735967858 1.209890249003069, -74.71056735999997 1.209890249000068, -74.71135265399997 1.209893721000071, -74.71056735967858 1.209890249003069))

Points and lines in blue, polygons in green

enter image description here

1
  • Thanks, Its possible to get polygons and create a shapefile only with these gemoetries? i.e. how subset the result and generate the desired polygon shapefile output? Commented Jul 13, 2018 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.