4

I have two geopandas data frames. For each row in the left frame, I would like to find which rows in the right frame overlap that row spatially, and by how much. Once I have this information, I'll be able to do a spatial join based on the degree of overlap.

Unfortunately, I'm doing this across a large number of polygons: all the US census tracts in a state (Texas has 5,265 of them) and a large number of polygons similar in size (but not coincident with) US census blocks (Texas has ~914,231 of these).

I'm looking for a way to do this faster. My code so far is below.

The datasets used can be acquired from the US Census: blocks data, tracts data.

#!/usr/bin/env python3

import geopandas as gpd
import geopandas_fast_sjoin as gpfsj
import time
import os
import pickle
import sys

os.environ["GDAL_DATA"] = "/usr/share/gdal"

TRACT_FILE    = "./data/tracts/tl_2010_{fips}_tract10.shp"
BLOCK_FILE    = "./data/blocks/tabblock2010_{fips}_pophu.shp"
PROJECTION    = '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

print("Reading data...")
start_time = time.time()
tracts = gpd.read_file(TRACT_FILE.format(fips=48))
blocks = gpd.read_file(BLOCK_FILE.format(fips=48))
print('Time: ', time.time()-start_time )

print("Converting coordinate reference systems...")
start_time = time.time()
tracts = tracts.to_crs(PROJECTION)
blocks = blocks.to_crs(PROJECTION)
print('Time: ', time.time()-start_time )

print("Performing spatial join...")
start_time = time.time()
joined = gpd.sjoin(tracts, blocks, how='left')
print('Time: ', time.time()-start_time )

print("Calculating areas of intersection...")
start_time = time.time()
joined['area_of_intersect'] = [row['geometry'].intersection(blocks.loc[row['index_right']]['geometry']).area for i,row in joined.iterrows()]
print('Time: ', time.time()-start_time )

1 Answer 1

7

There are several optimizations that would make this operation faster: doing all of the work in C++ without involving Python, using a spatial index to quickly identify candidates for intersections, using Prepared Geometries to quickly check the candidates, and parallelizing the entire operation over your available cores.

All of this could be accomplished in Python with the downside of some overhead, except that in my tests using the multiprocessing module on your multigigabyte dataset led Python to exhaust the available memory. On Windows this would probably be inevitable, on Linux copy-on-write should have prevented it. Maybe it can be done with careful programming, but the whole point of using Python is to not have to worry about these kinds of details. Therefore, I chose to move the computation to C++.

To accomplish this, I used pybind11 to build a new Python module that accepts lists of geometries from geopandas and produces an output of three lists: (1) Row indices of left-hand geometries; (2) Row indices of right-hand geometries; (3) Area of overlap between the two (only if >0).

For example, for an input with geometries left=[A,B,C,D] and right=[E,F,G,H], let:

  • E be entirely in A
  • F overlap with A and B
  • G and H overlap nothing

Then the return looks like:

List1 List2 List3
A     E     Area(E)
A     F     AreaIntersection(A,F)
B     F     AreaIntersection(B,F)

On my machine, the sjoin operation took 73s and calculating the intersections took 1,066s, for a total of 1139s (19 minutes)

On my 12 core machine, the code below takes 50s to do all of this.

So, for a spatial join in which the areas of intersection are needed, this saves only a little time. However, for a spatial join where the areas of intersection are needed, this saves a lot of time. Put another way, calculating all of those intersections takes a lot of work!

In further testing, calculating the areas of intersection without using the prepared geometries for acceleration took 287s on 12 cores. So parallelizing the intersections results in a 4x speed-up and parallelizing with prepared geometries results in a 23x speed-up.

Makefile

all:
    $(CXX) -O3 -g -shared -std=c++11 -I include `python3-config --cflags --ldflags --libs` quick_join.cpp -o geopandas_fast_sjoin.so -fPIC -Wall -lgeos -fopenmp 

quick_join.cpp

#define GEOS_USE_ONLY_R_API 1

#include <geos/geom/Geometry.h>
#include <geos/geom/prep/PreparedGeometry.h>
#include <geos/geom/prep/PreparedGeometryFactory.h>
#include <geos/index/strtree/STRtree.h>

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <memory>

#ifdef _OPENMP
#include <omp.h>
#else
#define omp_get_max_threads() 1
#define omp_get_thread_num() 0
#endif

///Fast Spatial Joins
///
///@params gp_left  List of GEOS geometry pointers from the left data frame.
///                 The code works best if gp_left is comprised of relatively
///                 fewer and relatively larger geometries.
///@params gp_right List of GEOS geometry pointers from the right data frame
///                 The code works best if gp_right is comprised of relatively
///                 more numerous and relatively smaller geometries.
///
///The list of GEOS geometry pointers can be acquired with
///    geos_pointers = [x._geom for x in df['geometry']]
///
///A common task in GeoPandas is taking two dataframes and combining their
///contents based on how much the contents' geometries overlap. However, this
///operation is slow in GeoPandas because most of it is performed in Python.
///Here, we offload the entire computation to C++ and use a number of techniques
///to achieve good performance.
///

///Namely, we create a spatial index from the left-hand geometries. For each
///geometry from the right-hand side, this allows us to very quickly find which
///geometries on the left-hand side it might overlap with. For each geometry on
///the left-hand side, we create a "prepared geometry", this accelerates simple
///spatial queries, such as checking for containment or disjointedness, by an
///order of magnitude. Finally, we parallelize the entire computation across all
///of the computer's threads.
///
///@return Three lists: (1) Row indices of left-hand geometries; (2) Row indices
///        of right-hand geometries; (3) Area of overlap between the two (only 
///        if >0).
///
///For example, for an input with geometries left=[A,B,C,D] and right=[E,F,G,H],
///let: 
/// * E be entirely in A
/// * F overlap with A and B
/// * G and H overlap nothing
///
///Then the return looks like:
///List1 List2 List3
///A     E     Area(E)
///A     F     AreaIntersection(A,F)
///B     F     AreaIntersection(B,F)
pybind11::tuple fast_sjoin(pybind11::list gp_left, pybind11::list gp_right) {
  typedef geos::geom::Geometry Geometry;
  typedef geos::geom::prep::PreparedGeometry PGeometry;

  //These are our return values
  std::vector<size_t> lefts;  //Indices of geometries from the left
  std::vector<size_t> rights; //Geometries of the overlapping geometries from the right
  std::vector<double> areas;  //Area of the overlap (if it is >0)

  //If either list is empty, there is nothing to do
  if(gp_left.size()==0 || gp_right.size()==0)
    return pybind11::make_tuple(lefts, rights, areas);

  //Used to cache pointers to geometries so we don't have to constantly be doing
  //conversions
  std::vector<const Geometry* > lgeoms;
  std::vector<const Geometry* > rgeoms;

  //Prepared geometries can massively accelerate computation by provide quick
  //predicate checks on whether one geometry contains another. Here, we stash
  //the prepared geometries. Unfortunately, GEOS 3.6.2's prepared geometries are
  //not reentrant. So we make a private set of prepared geometries for each
  //thread
  std::vector<std::vector<const PGeometry*>> lpgeoms(omp_get_max_threads());

  //Used for creating prepared geometries
  geos::geom::prep::PreparedGeometryFactory preparer;

  //For each geometry on the left, convert the input (a Python integer) into a
  //geometry pointer, then create a prepared geometry
  for(size_t i=0;i<gp_left.size();i++){
    const size_t    ptr_add = gp_left[i].cast<size_t>(); //Item is an integer
    const Geometry* geom    = reinterpret_cast<Geometry*>(ptr_add);
    for(int i=0;i<omp_get_max_threads();i++)
      lpgeoms.at(i).push_back(preparer.prepare(geom));
    lgeoms.push_back(geom);
  }

  //For each geometry on the right, convert the input (a Python integer) into a
  //geometry pointer
  for(size_t i=0;i<gp_right.size();i++){
    const size_t    ptr_add = gp_right[i].cast<size_t>(); //Item is an integer
    const Geometry* geom    = reinterpret_cast<Geometry*>(ptr_add);
    rgeoms.push_back(geom);
  }

  //The STRtree spatial index stores rectangles and allows us to quickly find
  //all the rectangles that overlap with a query rectangle. We leverage this
  //here by inserting the bounding boxes of all of the left geometries into a
  //spatial index. The spatial index also allows us to store a pointer to a data
  //structure; this pointer is returned if a query finds a hit or hits. We abuse
  //this capability by using the pointer to store the array index containing the
  //left geometry. This allows us to quickly find both the left geometry and its
  //associated prepared geometry.
    geos::index::strtree::STRtree index;
  for(size_t i=0;i<lgeoms.size();i++)
    index.insert(lgeoms[i]->getEnvelopeInternal(), reinterpret_cast<void*>(i));

  //Once all of the geometries are inserted into the spatial index, the index
  //must be built. This must be done in serial since GEOS does not have
  //protection against multiple threads trying it. (This is logical since it
  //eliminates a lock that would otherwise slow queries, but a better design
  //would probably have been to throw an exception.) The GEOS spatial tree also
  //lacks an explicit command for building the tree (wtf), so here we perform a
  //meaningless, single-threaded query to ensure the tree gets built.
  {
    std::vector<void *> results;
    index.query( lgeoms[0]->getEnvelopeInternal(), results );
  }

  //Each query to the spatial index populates a predefined vector with the
  //results of the query. We define this vector here, outside of the loop, to
  //avoid the memory cost of reallocating on each iteration of the loop.
  std::vector<void *> results;

  //These are custom OpenMP reduction operators for combining vectors together
  //following the parallel section. We leverage them to have each thread build
  //its own private result vectors which are afterwards combined into a single
  //result.
  #pragma omp declare reduction(merge : std::vector<uint64_t> : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
  #pragma omp declare reduction(merge : std::vector<double>   : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))

  //Now we loop through all of the geometries on the right-hand side. We do this
  //in parallel because we assume there are many of them.
  #pragma omp parallel for default(none) shared(rgeoms,index,lgeoms,lpgeoms) private(results) reduction(merge:lefts) reduction(merge:rights) reduction(merge:areas)
  for(unsigned int r=0;r<rgeoms.size();r++){
    const Geometry *const rgeom = rgeoms.at(r);

    index.query( rgeom->getEnvelopeInternal(), results );
    for(const auto q: results){
      //results is a list of "pointers". But we abused the pointers by using
      //them to stash array indices. Let's unpack the "pointers" into indices
      //now.
      const size_t lnum = reinterpret_cast<size_t>(q);
      const Geometry  *const lgeom  = lgeoms.at(lnum);
      const PGeometry *const lpgeom = lpgeoms.at(omp_get_thread_num()).at(lnum);

      if(lpgeom->contains(rgeom)){
        //The right-hand geometry is entirely inside the left-hand geometry
        lefts.push_back(lnum);
        rights.push_back(r);
        areas.push_back(rgeom->getArea());
      } else if(lpgeom->disjoint(rgeom)){
        //The right-hand geometry is entirely outside the left-hand geometry
      } else {
        //The right-hand geometry is partially inside and partially outside the
        //left-hand geometry, so we get the area of intersection of the two.
        std::unique_ptr<Geometry> igeom(rgeom->intersection(lgeom)); 
        const auto    iarea = igeom->getArea();               
        lefts.push_back(lnum);
        rights.push_back(r);
        areas.push_back(iarea);
      }
    }

    //Clear the results vector so we're ready for the next iteration. Note that
    //clearing it does not release its memory, so after the first few iterations
    //we should no longer be performing allocations.
    results.clear();
  }

  return pybind11::make_tuple(lefts, rights, areas);
}



PYBIND11_MODULE(geopandas_fast_sjoin,m){
  m.doc() = "Fast spatial joins";

  m.def("fast_sjoin", &fast_sjoin, "Performs a fast spatial join");
}

test.py

#!/usr/bin/env python3

import geopandas as gpd
import geopandas_fast_sjoin as gpfsj
import time
import os
import pickle
import sys

os.environ["GDAL_DATA"] = "/usr/share/gdal"

DATA_DIR      = "./data/"
TRACT_FILE    = "./data/tracts/tl_2010_{fips}_tract10.shp"
BLOCK_FILE    = "./data/blocks/tabblock2010_{fips}_pophu.shp"
PRECINCT_FILE = "./data/precincts/precincts2008/USA_precincts.shp"
STATES_FILE   = "./data/states/tl_2010_us_state10.shp"
PROJECTION    = '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

print("Reading data...")
start_time = time.time()
tracts = gpd.read_file(TRACT_FILE.format(fips=48))
blocks = gpd.read_file(BLOCK_FILE.format(fips=48))
print('Time: ', time.time()-start_time )

print("Converting coordinate reference systems...")
start_time = time.time()
tracts = tracts.to_crs(PROJECTION)
blocks = blocks.to_crs(PROJECTION)
print('Time: ', time.time()-start_time )

print("Performing spatial join...")
start_time = time.time()
joined = gpd.sjoin(tracts, blocks, how='left')
joined['area_of_intersect'] = [row['geometry'].intersection(blocks.loc[row['index_right']]['geometry']).area for i,row in joined.iterrows()]
print('Time: ', time.time()-start_time )

# pickle.dump( (blocks,tracts), open( "save.p", "wb" ) )
# sys.exit(-1)
# blocks, tracts = pickle.load( open( "save.p", "rb" ) )

print("Getting geometries...")
start_time = time.time()
tgeoms = [x._geom for x in tracts['geometry']]
bgeoms = [x._geom for x in blocks['geometry']]
print('Time: ', time.time()-start_time )

print("Running example")
start_time = time.time()
lefts, rights, areas = gpfsj.fast_sjoin(tgeoms,bgeoms)
print('Time: ', time.time()-start_time )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.