I have two shapefiles, a province with 39 municipalities and a land cover of the province.
Goal: get the sqm of each land cover classification for each municipality. For example, how much open forest, closed forest, eroded areas, etc. are there in Muni 1.
I managed to write a code that does this but only for one of the 39 municipalities.
gdf1 = gpd.read_file('munis.shp') gdf2 = gpd.read_file('landcover.shp') muni1 = gdf1.loc[gdf1['NAME_2'] == 'muni1'] # NAME_2 is field containing all munis intersection = gpd.overlay(muni1, gdf2, how='intersection') area = intersection['geometry'].map(lambda p: p.area) print(area)
The result is the sqm of each land cover classification in
muni1. How to insert this in a for-loop so I can do the process for each 39 municipalities? I tried this but not working:
for i in range(gdf1.shape): munis = gdf1.loc[gdf1['NAME_2'][i]] intersection = gpd.overlay(munis, gdf2, how='intersection') area = intersection['geometry'].map(lambda p: p.area) print(area) Result KeyError: 'Abad' # first municipality if sorted alphabetically