Linking two GeoJSON data files using Python?

Question

I am new to using GIS data in Python. I have two datasets in GeoJSON format. One of them has states of a country with multi-polygons representing each states. The other has multiple points of latitude and longitude, each one describing a particular branch information (also in the same country).

I am looking to integrate these two GeoJSON files so that I can get something like all the branches or the number of branches located in a particular state.

How can this be done in Python?

I've attached the code so far and the results that I've obtained below.

import geopandas as gpd
import geojson, json

fname_states="/home/anaconda/gis_test/USAStates.json"
fname_branches="/home/anaconda/gis_test/branches.json"
with open(fname_states) as fs:
    data_states = json.load(fs)
with open(fname_branches) as fb:
    data_branches = json.load(fb)
for feature in data_states['features']:
    print (feature['geometry']['type'])
    print (feature['geometry']['coordinates'])

MultiPolygon [[[[-155.666192384, 18.921786345999976], ............, [-160.073803717, 22.00417732599999]]]]

MultiPolygon [[[[-122.40201532499998, 48.22521676600002],.............., [-122.96797834099999, 48.443794749999995]]]]

I've just attached two polygons output statement for reference. These coordinates indicate polygons which refer to a region.

for feature in data_branches['features']:
    print (feature['geometry']['type'])
    print (feature['geometry']['coordinates'])

Point [-84.11397141, 35.97929825] Point [-84.17377774, 35.89437917] Point [-83.84168579, 36.22807494] Point [-83.92360125, 36.08801205] Point [-84.02843334, 35.93056813] Point [-84.00351839, 36.04683377] Point [-83.97332803, 35.97456329]

This is the output of the branches file which also contain numerous other fields like branch id and other attributes.

So, I would like to add an extra field to the branches data as 'state'. I want to obtain this field by checking whether it lies in any of the regions that are covered by polygons in the states GeoJSON file. This is where I need help.

Please suggest an optimal method to achieve this.

Answer 1

It sounds like you want to do a geometry intersection to 'join' two datasets by location.

This website has an example of what you want to do, and it looks something like this:

import json
from shapely import shape, Point

# load GeoJSON file containing sectors
with ('sectors.json', 'r') as f:
    js = json.load(f)

# construct point based on lat/long returned by geocoder
point = Point(45.4519896, -122.7924463)

# check each polygon to see if it contains the point
for feature in js['features']:
    polygon = shape(feature['geometry'])
    if polygon.contains(point):
        print 'Found containing polygon:', feature

Check out the Shapely docs for more info.