City-Lifeline-Cockpit/tools/city-massing/cad-shp-to-geojson.py

# -*- coding: utf-8 -*-
"""Convert the exported CAD building layer to GeoJSON and report what actually lands in the scene.

Run with ArcGIS Desktop Python 2.7:
    python.exe tools/city-massing/cad-shp-to-geojson.py <inDir> <outDir> <west,south,east,north>

Only features whose envelope centre falls in the given window are written, so the GeoJSON stays a
scene-sized extract of the full survey rather than a county-wide dump.
"""
import sys
import os
import io
import json
import arcpy


def main():
    if len(sys.argv) < 4:
        print('usage: cad-shp-to-geojson.py <inDir> <outDir> <west,south,east,north>')
        return 2
    in_dir = os.path.abspath(sys.argv[1])
    out_dir = os.path.abspath(sys.argv[2])
    west, south, east, north = [float(v) for v in sys.argv[3].split(',')]
    if not os.path.isdir(out_dir):
        os.makedirs(out_dir)

    sources = [
        ('Polygon-JMD.shp', 'cad-buildings.geojson', 'polygon'),
        ('Annotation-JMD.shp', 'cad-building-labels.geojson', 'annotation'),
        ('Polyline-JMD.shp', 'cad-building-outlines.geojson', 'polyline'),
    ]
    report = {'window': [west, south, east, north], 'layers': []}
    for filename, out_name, kind in sources:
        path = os.path.join(in_dir, filename)
        if not os.path.exists(path):
            print('missing %s' % path)
            continue
        fields = [f.name for f in arcpy.ListFields(path)]
        text_field = 'TextString' if 'TextString' in fields else None
        shape_type = arcpy.Describe(path).shapeType
        features = []
        in_window = 0
        total = 0
        areas = []
        with arcpy.da.SearchCursor(path, ['SHAPE@'] + ([text_field] if text_field else [])) as cursor:
            for row in cursor:
                total += 1
                shape = row[0]
                if shape is None:
                    continue
                extent = shape.extent
                cx = (extent.XMin + extent.XMax) / 2.0
                cy = (extent.YMin + extent.YMax) / 2.0
                if not (west <= cx <= east and south <= cy <= north):
                    continue
                in_window += 1
                geometry = None
                if shape_type == 'Polygon' and kind == 'annotation':
                    # Annotation geometry is the label box; the point is the useful location.
                    anchor = shape.labelPoint if hasattr(shape, 'labelPoint') else None
                    geometry = {'type': 'Point', 'coordinates': [round(anchor.X, 7), round(anchor.Y, 7)]} if anchor else None
                elif shape_type == 'Polygon':
                    rings = []
                    for i in range(shape.partCount):
                        part = shape.getPart(i)
                        ring = [[round(p.X, 7), round(p.Y, 7)] for p in part if p is not None]
                        if len(ring) >= 4:
                            rings.append(ring)
                    if rings:
                        geometry = {'type': 'Polygon', 'coordinates': [rings[0]] + [[r] for r in rings[1:]]}
                        try:
                            areas.append(float(shape.area))
                        except Exception:
                            pass
                elif shape_type == 'Polyline':
                    lines = []
                    for i in range(shape.partCount):
                        part = shape.getPart(i)
                        line = [[round(p.X, 7), round(p.Y, 7)] for p in part if p is not None]
                        if len(line) >= 2:
                            lines.append(line)
                    if lines:
                        geometry = {'type': 'MultiLineString', 'coordinates': lines}
                if geometry is None:
                    continue
                props = {}
                if text_field:
                    props['text'] = row[1]
                features.append({'type': 'Feature', 'properties': props, 'geometry': geometry})

        out_path = os.path.join(out_dir, out_name)
        with io.open(out_path, 'w', encoding='utf-8') as handle:
            handle.write(unicode(json.dumps({'type': 'FeatureCollection', 'features': features}, ensure_ascii=False)))
        entry = {'source': filename, 'out': out_name, 'total': total, 'inWindow': in_window, 'written': len(features)}
        if areas:
            areas.sort()
            pick = lambda q: round(areas[int(q * (len(areas) - 1))], 1)
            entry['areaSqm'] = {'min': round(areas[0], 1), 'p25': pick(.25), 'median': pick(.5),
                                'p75': pick(.75), 'p95': pick(.95), 'max': round(areas[-1], 1),
                                'under15': sum(1 for a in areas if a < 15),
                                'over2000': sum(1 for a in areas if a > 2000)}
        report['layers'].append(entry)
        print('%-22s total=%6d inWindow=%6d written=%6d' % (filename, total, in_window, len(features)))
        if 'areaSqm' in entry:
            print('   area %s' % json.dumps(entry['areaSqm']))

    with io.open(os.path.join(out_dir, 'cad-extract-report.json'), 'w', encoding='utf-8') as handle:
        handle.write(unicode(json.dumps(report, ensure_ascii=False, indent=2, sort_keys=True)))
    return 0


if __name__ == '__main__':
    sys.exit(main())