City-Lifeline-Cockpit/tools/city-massing/export-cad-layers.py

# -*- coding: utf-8 -*-
"""Report extents and area statistics for chosen CAD layers, and export them to shapefile.

Run with ArcGIS Desktop Python 2.7:
    python.exe tools/city-massing/export-cad-layers.py <cad.gdb> <outDir> [layer,layer,...]

Only the named layers are exported, so the delivered shapefile is a deliberate subset rather than
a dump of the whole survey drawing.
"""
import sys
import os
import io
import json
import arcpy


def polygon_area_stats(layer):
    areas = []
    with arcpy.da.SearchCursor(layer, ['SHAPE@AREA']) as cursor:
        for row in cursor:
            value = row[0]
            # Drawing-space entities (title blocks, symbols) have no meaningful area here.
            if value is None:
                continue
            areas.append(float(value))
    if not areas:
        return {}
    areas.sort()
    pick = lambda p: round(areas[int(p * (len(areas) - 1))], 1)
    return {
        'count': len(areas),
        'minSqm': round(areas[0], 1),
        'p25Sqm': pick(0.25),
        'medianSqm': pick(0.5),
        'p75Sqm': pick(0.75),
        'p95Sqm': pick(0.95),
        'maxSqm': round(areas[-1], 1),
        'under15Sqm': sum(1 for a in areas if a < 15),
        'over3000Sqm': sum(1 for a in areas if a > 3000),
    }


def main():
    if len(sys.argv) < 3:
        print('usage: export-cad-layers.py <cad.gdb> <outDir> [layer,layer,...] [window=XMIN,YMIN,XMAX,YMAX]')
        return 2
    gdb = os.path.abspath(sys.argv[1])
    out_dir = os.path.abspath(sys.argv[2])
    wanted = None
    window = None
    for arg in sys.argv[3:]:
        if arg.startswith('window='):
            window = [float(v) for v in arg[len('window='):].split(',')]
        elif wanted is None:
            wanted = arg.split(',')
    if not os.path.isdir(out_dir):
        os.makedirs(out_dir)
    arcpy.env.workspace = gdb

    targets = []
    for dataset in arcpy.ListDatasets() or []:
        for name in arcpy.ListFeatureClasses(feature_dataset=dataset) or []:
            targets.append(name)

    # The drawing carries no coordinate system, but its survey coordinates are CGCS2000
    # 3-degree Gauss-Kruger zone CM 111E (EPSG:4546), the same system as the sibling shapefiles.
    SURVEY_SR = arcpy.SpatialReference(4546)
    report = {'gdb': gdb, 'assumedSpatialReference': 'CGCS2000_3_Degree_GK_CM_111E (EPSG:4546)', 'window': window, 'layers': []}
    print('assumed spatial reference: %s' % report['assumedSpatialReference'])

    def in_window(shape):
        if window is None:
            return True
        try:
            extent = shape.extent
            cx = (extent.XMin + extent.XMax) / 2.0
            cy = (extent.YMin + extent.YMax) / 2.0
        except Exception:
            # MultiPatch and other 3D/unsupported shapes fall back to their envelope.
            return True
        return (window[0] <= cx <= window[2]) and (window[1] <= cy <= window[3])

    for name in targets:
        describe = arcpy.Describe(name)
        entry = {'featureClass': name, 'featureType': describe.shapeType, 'featureCount': 0, 'inWindow': 0}
        areas = []
        with arcpy.da.SearchCursor(name, ['SHAPE@', 'Layer'] if 'Layer' in [f.name for f in arcpy.ListFields(name)] else ['SHAPE@']) as cursor:
            for row in cursor:
                entry['featureCount'] += 1
                shape = row[0]
                if shape is None or not in_window(shape):
                    continue
                entry['inWindow'] += 1
                if describe.shapeType == 'Polygon':
                    try:
                        if shape.area is not None:
                            areas.append(float(shape.area))
                    except Exception:
                        pass
        if areas:
            areas.sort()
            pick = lambda p: round(areas[int(p * (len(areas) - 1))], 1)
            entry['areaInWindow'] = {
                'count': len(areas), 'minSqm': round(areas[0], 1), 'p25Sqm': pick(0.25),
                'medianSqm': pick(0.5), 'p75Sqm': pick(0.75), 'p95Sqm': pick(0.95),
                'maxSqm': round(areas[-1], 1),
                'under15Sqm': sum(1 for a in areas if a < 15),
                'over3000Sqm': sum(1 for a in areas if a > 3000),
            }
        report['layers'].append(entry)
        print('%-12s %-9s total=%6d inWindow=%6d' % (name, entry['featureType'], entry['featureCount'], entry['inWindow']))
        if 'areaInWindow' in entry:
            print('             %s' % json.dumps(entry['areaInWindow']))

        if not wanted:
            continue
        fields = [f.name for f in arcpy.ListFields(name)]
        if 'Layer' not in fields:
            continue
        # MultiPatch repeats the polygon geometry in 3D and Annotation is text placement;
        # neither is a footprint source, so only true 2D geometry is exported.
        if describe.shapeType not in ('Polygon', 'Polyline', 'Point'):
            print('  skipping %s (%s)' % (name, describe.shapeType))
            continue
        for layer_name in wanted:
            safe = layer_name.encode('ascii', 'replace').replace('?', '_')
            out_shp = os.path.join(out_dir, '%s-%s.shp' % (name, safe))
            where = "Layer = '%s'" % layer_name.replace("'", "''")
            selected = os.path.join(out_dir, '_sel.shp')
            for stale in (selected, out_shp):
                if arcpy.Exists(stale):
                    arcpy.Delete_management(stale)
            arcpy.Select_analysis(name, selected, where)
            n = int(arcpy.GetCount_management(selected).getOutput(0))
            if n == 0:
                arcpy.Delete_management(selected)
                print('  no features for layer %s in %s' % (layer_name, name))
                continue
            # The drawing carries no CRS: declare the survey system, then reproject to WGS84.
            # Spatial clipping happens downstream in Node, where Turf already owns those rules.
            arcpy.DefineProjection_management(selected, SURVEY_SR)
            arcpy.Project_management(selected, out_shp, arcpy.SpatialReference(4326))
            arcpy.Delete_management(selected)
            final = int(arcpy.GetCount_management(out_shp).getOutput(0))
            print('  exported %-10s %-9s %6d features -> %s' % (layer_name, name, final, os.path.basename(out_shp)))
            entry.setdefault('exported', []).append({'layer': layer_name, 'features': final, 'shapefile': os.path.basename(out_shp)})

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


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