2

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[0]):
    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
6

loc primarily "label of the index" based. So, for the first iteration gdf1['NAME_2'][0] returns muni1. Then you get error using gdf1.loc['muni1'], because Pandas interprets 'muni1'('Abad') as index label, but it's value of NAME_2

You should change munis = gdf1.loc[gdf1['NAME_2'][i]]
into gdf1.loc[gdf1['NAME_2'] == gdf1['NAME_2'][i]] in for loop.

| improve this answer | |
  • It worked! Thank you! I have so much to learn still in using for-loops. Thanks for the explanation too. Really heplful, @KadirSahbaz. – BallpenMan Jun 13 at 0:35

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.