2

I am new to working with geodata.

Setup: GeoPandas, Python3.6

I have a dataset of square polygons of 10x10km for countries. However, some of these squares are already completely in a neighboring country. To identify and later remove the ones outside the country border. Hence I have two GeoDataFrames:

  1. squares >> contains all squares plus a country identifier.
  2. countries >>> contains the country borders as polygon or multipolygon

What puzzles me is that the below code works fast for some countries, but takes forever for others. The number of square polygons cannot be the issue as the processing times for countries with approximately 1200 squares can be 45 seconds or 30 minutes. I suspect that countries with an intricate border (e.g. Sweden and its archipelago) take longer to process.

The code never stops, but it just takes forever. I am thus now running each country one by one rather than the complete list of squares. Each country is saved as a pickle file after processing to clean out memory as it goes up to 7-8GB per country. Below is my code for Sweden which has some 23,000 squares.

How can I improve speed? I am assigning a spatial index with df.sindex. I also looked at R-Tree as used by Geoff Boeing, but that joins points in polygons. I could not make it work for polygons in polygons. And maybe there is even a better way of doing this.

import pandas as pd
import geopandas as gp
from matplotlib import pyplot as plt
import numpy as np

### Data overview
countries.head(2)
     COUNTRY                                           geometry
0         SE  MULTIPOLYGON (((24.15513 65.81603, 24.12993 65...
1         FI  MULTIPOLYGON (((28.92968 69.05190, 28.82917 69...

# Sample data for Sweden
squares[squares['COUNTRY'] == 'SE'].head(2)
       COUNTRY                                           geometry
717995      SE  POLYGON ((9.64431 58.11107, 9.64343 58.20099, ...
717996      SE  POLYGON ((9.64343 58.20099, 9.64254 58.29091, ...

### Processing
# Subset country
country = 'SE' # SE for Sweden
mask = squares['COUNTRY'] == country 
squares_inter = squares.loc[mask, :]
# Spatial join
shape_inter = gp.sjoin(shapeisect, countries, how='inner', op='intersects')
# Amending columns
shape_inter.rename(columns={'COUNTRY_left':'COUNTRY_sq',
                            'COUNTRY_right':'COUNTRY_c'}, inplace=True)  
# Get match
shape_inter['match'] = np.where(shape_inter['COUNTRY_sq'] == shape_inter['COUNTRY_c'], 1, np.nan)
# Save to pickle
filename = 'country  + '.pkl'
shape_inter.to_pickle(filename)

### Plotting squares and country borders BEFORE processing
mask = countries['COUNTRY'] == 'SE'
base = countries.loc[mask, ['geometry']].plot(figsize=(10,10), facecolor="none", edgecolor='black')
mask = (squares['COUNTRY'] == 'SE')
squares.loc[mask, ['geometry']].plot(ax=base, facecolor='none', edgecolor='green')

Below is an image that shows how the squares go beyond country borders. enter image description here

0

It seems it was a data problem after all. Countries with archipelagoes did take longer. At least I could speed up the process a little bit by limiting .sjoin to just the country I was assessing. Here is the snippet I added:

# Spatial join
mask = countries['CNTR_CODE'] == country
countryisect = countries.loc[mask,['geometry', 'CNTR_CODE']]
shape_inter = gp.sjoin(shapeisect, countryisect , how='inner', op='intersects')
| improve this answer | |

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.