Is there a python option to "join attributes by location"?

Question

I am trying to perform the function join attributes by location as found on the QGIS menu Vector>Data Management Tools. I am seeking an open source python option for this. I know arcpy has a spatial join function but I'm trying to perform this outside of the ESRI environment.

I would suggest to look at the source code of the actual Join attributes by location command from the fTools plugin: doSpatialJoin.py, particularly the compute() method. Shouldn't be too hard to eliminate any UI code from that and strip it down to a plain python function. — Lukas Graf, Oct 25 '13 at 22:32
Hi i have a little bit different problem , i want to check if there is a joint between two layer!I am seeking an open source python option for this. I want to test if the s.join function has been used and I'm trying to perform this outside of the ESRI environment.. Can anyone help me please ! — Rania ben othmen, Nov 30 '17 at 11:27

score 9 · Accepted Answer · answered Oct 27 '13 at 19:11

9

You may want to take a look at Shapely and Fiona. Fiona is a wrapper for gdal to make spatial file import and export easy. Shapely provides geometry functionality. Here is a very simple example to give you the idea. It joins polygon attributes to all points within that polygon.

The example data I have used are these polygons and these points.

import fiona
from shapely.geometry import shape
from copy import deepcopy

with fiona.open("planning_neighborhoods.shp", "r") as n: 

    with fiona.open("Schools_Private_Pt.shp", "r") as s:

        # create a schema for the attributes
        outSchema =  deepcopy(s.schema)
        outSchema['properties'].update(n.schema['properties'])

        with fiona.open ("Schools_withNbhd.shp", "w", s.driver, outSchema, s.crs) as output:

            for school in s: 
                for neighborhood in n:
                    # check if point is in polygon and set attribute
                    if shape(school['geometry']).within(shape(neighborhood['geometry'])):  
                        school['properties']['neighborho'] = neighborhood['properties']['neighborho'] 
                    # write out
                        output.write({                                 
                            'properties': school['properties'], 
                            'geometry': school['geometry']
                        })

answered Oct 27 '13 at 19:11

cengel

3,274
2
18
30

Thanks @cengel. This looks like it will put me on the right track! I'm actually interested in the join with lines and polygons (specifically finding where rivers intersect model cells) and I think that will work following your example. – mishaF Oct 28 '13 at 15:26
@cengel Do qgis stations using plugins employing these methods require installation of modules and gdal? – user25976 May 05 '14 at 22:36
@user25976 sorry, not sure I quite understand your question. My code example is a standalone python script. Both fiona and shapely require gdal. – cengel May 07 '14 at 00:36
@cengel Excuse me, let me clarify (I'm new to programming). Regarding standalone python scripts: you mean that a plugin written with fiona and shapely imports can be used by QGIS users even if they don't have python or the modules installed on their computer? – user25976 May 07 '14 at 00:47
@user25976 They do need the modules installed on their computer. See for example here – cengel May 07 '14 at 06:16
Would it be possible to have the neighborhoods as the output? I'm trying to do something where I want to use points to select polygons http://gis.stackexchange.com/questions/125043/finding-an-ogr2ogr-equivalent-of-selectlayerbylocation not using ArcGIS. I tried switching the input and the output shapefiles in this code without luck. – GIS Danny Dec 11 '14 at 04:01
Think I got it - see here http://gis.stackexchange.com/questions/125043/finding-an-ogr2ogr-equivalent-of-selectlayerbylocation/125227#125227 – GIS Danny Dec 11 '14 at 05:09
Anyone looking at this (fine) example, be aware that for bigger analyses, you will want to use a spatial index like rtree. There is plenty of information about this on GISSA for your enjoyment. :) – bugmenot123 Jul 30 '18 at 09:24

score 3 · Answer 2 · answered Aug 28 '17 at 15:56

Although still a bit rough around the edges, especially when it comes to documentation and examples, but geopandas' future looks bright. It basically combines the power of pandas dataframes with geospatial capabilities of shapely.

the function you look for is called sjoin

Make sure your machine/instance has enough memory to perform the operation

import geopandas as gpd
import pandas as pd
import os

gdfLeft = gpd.read_file(os.path.join(PATH,INPUT_FILE_NAME_1))
gdfRight = gpd.read_file(os.path.join(PATH,INPUT_FILE_NAME_2))

gdfJoined = gpd.sjoin(gdfLeft, gdfRight, how="left", op='intersects')

This piece of code performs the spatial join but the attributes of the joined shapefile is empty. Any pointers? — April, Oct 16 '19 at 11:27

Is there a python option to "join attributes by location"?

2 Answers2

Linked