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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
gonzalez.ivan90
  • 283
  • 3
  • 11

1 Answers1

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

gene
  • 54,868
  • 3
  • 110
  • 187
  • 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? – gonzalez.ivan90 Jul 13 '18 at 15:25