add an "area" field and a "centroid" field to feature.geometry for (multi)polygons in a geojson file

gj-area-centroid.py

#!/usr/bin/python
# coding=utf-8
# -*- coding: utf-8 -*-

# This script adds an "area" field (in square meters)
# and a "centroid" field to feature.geometry
# for Polygon's and MultiPolygon's in a geojson file.
# Usage example:
# ./gj-area-centroid.py university.geojson

# The computation is not exact. It is a planar approximation.
# Yet the error is typically less than 1/100000 for areas
# no larger than a few square kilometers in lower latitudes.
# (My 243 test cases are universities in Taiwan,
# between 22-25 degrees North.)

# If you have this true area computation package
# https://pypi.python.org/pypi/area
# installed in your system, you can add the -e flag, e.g.
# ./gj-area-centroid.py -e 100000 university.geojson
# to see the error rate (multiplied by the given number)
# as compared to the result of the (true) spherical formula.

import argparse, re, warnings, atexit, json, sys, math

def ring_area_centroid(ringlist):
    # Shift all coords to be relative to the approximate
    # middle of the list of rings so that the computation
    # will be more numerically stable.
    all_x = [v[0] for ring in ringlist for v in ring]
    all_y = [v[1] for ring in ringlist for v in ring]
    mid = [(max(all_x)+min(all_x))/2, (max(all_y)+min(all_y))/2]
    area = cx = cy = 0
    # no need to distinguish between the outer rings and
    # the inner rings because the inner rings' areas
    # will be negative anyway
    for c in ringlist:
        n = len(c)
        # last point is the same as the first,
        # so no need to wrap around.
        for i in range(0, n-1):
            x0 = c[i][0] - mid[0]
            y0 = c[i][1] - mid[1]
            x1 = c[i+1][0] - mid[0]
            y1 = c[i+1][1] - mid[1]
            cross = x0*y1 - x1*y0
            area += cross
            # https://en.wikipedia.org/wiki/Centroid#Centroid_of_a_polygon
            cx += (x0 + x1) * cross
            cy += (y0 + y1) * cross
    # https://en.wikipedia.org/wiki/World_Geodetic_System
    r = 6378137*math.pi/180
    return (
        area*r*r*math.cos(mid[1]*math.pi/180)/2,
        [mid[0]+cx/3/area, mid[1]+cy/3/area]
    )

def geom_area_centroid(geom):
    if geom['type'] == 'MultiPolygon':
        ringlist = [ring for poly in geom["coordinates"] for ring in poly]
    else:
        ringlist = geom["coordinates"]
    return ring_area_centroid(ringlist)

def r_area_centroid(data, areaFN='area', centroidFN='centroid'):
    # recursively compute area and centroid
    if type(data) is list:
        return [r_area_centroid(x, areaFN=areaFN, centroidFN=centroidFN) for x in data]
    elif type(data) is dict:
        if not 'type' in data:
            return data
        if 'Feature' in data['type']:   # "Feature" or "FeatureCollection"
            return { k: r_area_centroid(data[k], areaFN=areaFN, centroidFN=centroidFN) for k in data.keys() }
        elif 'Polygon' in data['type']: # "Polygon" or "MultiPolygon"
            (data[areaFN], data[centroidFN]) = geom_area_centroid(data)
            if args.error > 0:
                data['area_true'] = area.area(data)
                data['area_error'] = (data[areaFN] - data['area_true'])/data['area_true'] * args.error
            return data
        else:
            return data
    else:
        return data

parser = argparse.ArgumentParser(
    description='compute area of polygons in a json file',
    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('-a', '--area', type=str,
    default='area', help='field name of area')
parser.add_argument('-c', '--centroid', type=str,
    default='centroid', help='field name of centroid')
parser.add_argument('-e', '--error', type=float,
    default=-1, help='comput error rate and multiply it by this positive number. (e.g. "-e 100" for "percentage") needs the area package.')
parser.add_argument('FILE', help='the geojson file')
args = parser.parse_args()

if args.error > 0:
    import area
with open(args.FILE, 'r') as gjfile:
    data = json.load(gjfile)

# https://stackoverflow.com/a/5648769
# Python print unicode strings in arrays as characters, not code points
# https://stackoverflow.com/questions/2596714/why-does-python-print-unicode-characters-when-the-default-encoding-is-ascii
# print repr(r_area_centroid(data)).decode('unicode-escape')

# https://stackoverflow.com/questions/18337407/saving-utf-8-texts-in-json-dumps-as-utf8-not-as-u-escape-sequence
print json.dumps(r_area_centroid(data, areaFN=args.area, centroidFN=args.centroid), ensure_ascii=False, indent=2).encode('utf8')