Merge counties represented by multiple Polygon features into a single MultiPolygon feature

join_county.py

#! /usr/bin/python

import simplejson as json

def join_by_fips(counties, fips):
    """Joins a county by the given FIPS code into a single multipolygon geometry,
       it is assumed that each feature has properties as seen in the the us-atlas 
       shapefiles see: https://github.com/mbostock/us-atlas
    """

    joined_state = {}
    joined_state["type"] = counties["type"]
    joined_state["features"] = []
    joined_properties = {
            'AREA': 0,
            'COUNTY': '',
            'COUNTYP010':0,
            'FIPS':fips,
            'PERIMETER':0,
            'SQUARE_MIL':0,
            'STATE':'',
            'STATE_FIPS':''
            }
    cords = []

    # loop through collect coordinates for duplicate counties
    for feature in counties['features']:
        if feature['properties']['FIPS'] == fips:

            cords.append(feature['geometry']['coordinates'])
            joined_properties['AREA'] += feature['properties']['AREA']
            joined_properties['PERIMETER'] += feature['properties']['PERIMETER']
            joined_properties['SQUARE_MIL'] += feature['properties']['SQUARE_MIL']
            joined_properties['COUNTY'] = feature['properties']['COUNTY']
            joined_properties['COUNTYP010'] = feature['properties']['COUNTYP010']
            joined_properties['STATE'] = feature['properties']['STATE']
            joined_properties['STATE_FIPS'] = feature['properties']['STATE_FIPS']

        else:
            # rebuild other features
            joined_state["features"].append(feature)

    # create the newly joined county feature of type multipolygon
    joined_county = {
            "type": "Feature",
            "geometry": {
                "type": "MultiPolygon",
                "coordinates": cords
                },
            "properties": joined_properties,
            }
    # add it to the feature list
    joined_state["features"].append(joined_county)
    return joined_state

if __name__ == '__main__':
    import sys

    county_file = sys.argv[1]
    fips = sys.argv[2]
    
    with open(county_file, 'r') as f:
        counties = json.load(f)

    joined_counties = join_by_fips(counties, fips)
    with open("joined_"+county_file, 'w+') as f:
        json.dump(joined_counties, f)