4

First of all, some info about my system characteristics:

  • Python 2.7
  • GeoPandas 0.3.0

I am working on a program, where I filter out the data of two shapefiles -link with data- (each shapefile is a vector layer with polygons) and organize it according some statements. The two shapefiles look like the picture below, the picture is just a small part of the whole map.

Link to the picture

After doing some filtering, I need to find in the remaining data those polygons inside the big polygons (in the picture the purple and the green polygon) and change some information if there is an intersection. I have tried using some of the GeoPandas functions for this task, but the results are not the ones that I expected.

Try 1:

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

#Small parts of the shape files
la_g = la.groupby('fclass').get_group('residential')
b_r = b[100:130]  # here I know there is an intersection with the big polygon

# Find the polygons inside big polygons 
inter = gpd.overlay(b_r, la_g, how='intersection')

Here the work is well done, but the function returns me another shapefile and it takes a lot computing time if a use the hole data in the shapefiles.

Try 2:

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

# Find the polygons inside big polygons 
inter = b.within(la.loc[0, 'geometry'])

In “try 2” I got a Boolean series with everything false, although I know that there are polygons inside the big ones. I tried to use instead the polygons their centroid and make the same calculations, but the result was the same, a Boolean series with all false

# Find the polygons inside big polygons 
b_c = b.centroid
inter = b_c.within(la.loc[0, 'geometry'])

I would prefer a boolean series, because I just need to change some information in the shapefile if there is an intersection. The idea is something like this (the code below is a matter illustration from my idea):

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

# Find intersections
inter = b.within(la.loc[0, 'geometry'])

# Comparation
for r in range(len(b))  
    if inter[r] == True:
        b[r]['type'] = 'residential'

How can I do that or have a better idea that achieves the result that I expect?

2

I think what you want is a "spatial join", where we add information of the second dataframe (the land use dataframe, information on the type of land use) to the first dataframe (the buildings) based on their spatial relationship:

gpd.sjoin(b, la[['fclass', 'geometry']], how='left', op='intersects')

This adds a new column to the b dataframe with the land use class in which each building is located.


Some additional notes about your attempts:

Try 1: a overlay can be suited here if you want to create a new geometry (e.g. in the case where a building is only partly inside a big land use polygon, then overlay will actually clip the original builing's geometry). But if you don't want to edit the geometries (take the intersection), the the spatial join (sjoin) shown above is more appropriate.
It is true that overlay is currently very show, as you noticed. There is work underway to improve that in the next version of geopandas: https://github.com/geopandas/geopandas/pull/429.

Try 2: this should also work, and actually for me it does not give a full series of False values. For the exact example you gave (inter = b.within(la.loc[0, 'geometry'])), I get 18 True values (from inter.sum()). However, this is actually equivalent to the spatial join, apart from that the spatial join does this automatically for all elements of the land use dataset.

| improve this answer | |
  • Thank you very much!!!! Your idea with “spatial join” is better suited to my goal extracting the information from the second layer (“la”) base on the position plus it creates another column that I can use later in my program. I have no idea why I get everything in false (I am a beginner in the python world), but in any case, i can know continue working on my program. – Yiyi Jun 19 '18 at 7:52

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.