"""Prepare buildings exclusively from the user-supplied VWorld SHP archive.

python scripts/prepare-vworld-buildings.py PATH_TO_ZIP
Requires the existing .tools/terrain environment plus building-requirements.txt.
"""
import argparse
import collections
import copy
import hashlib
import json
from pathlib import Path
import struct
import sys
import zipfile

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT / '.tools/terrain'))
import shapefile
from pyproj import CRS, Transformer
from shapely.geometry import shape, box
from shapely.ops import transform
from shapely.strtree import STRtree

CONFIG = json.loads((ROOT / 'scripts/constants/gwacheon-buildings.json').read_text(encoding='utf-8'))
TO_LOCAL = Transformer.from_crs(4326, CONFIG['sourceCrs'], always_xy=True).transform
TO_GEO = Transformer.from_crs(CONFIG['sourceCrs'], 4326, always_xy=True).transform


def positive(value, maximum):
    try:
        number = float(value)
        return number if 0 < number <= maximum else None
    except (ValueError, TypeError):
        return None


def extract_archive(path):
    """Read SHPs sequentially and only decompress DBF up to selected records.

    pyshp.record() rewinds its DBF on every call; doing that on a compressed
    1.8 GB member repeatedly is prohibitive. Read the selected fixed-width rows
    in ascending order instead. Source files remain in the supplied ZIP.
    """
    bounds = transform(TO_LOCAL, box(*CONFIG['bbox'])).buffer(20).bounds
    features = []
    with zipfile.ZipFile(path) as archive:
        for name in archive.namelist():
            if not name.lower().endswith('.shp'):
                continue
            stem = name[:-4]
            source_crs = CRS.from_wkt(archive.read(stem + '.prj').decode())
            expected_crs = CRS(CONFIG['sourceCrs'])
            # ESRI WKT uses a different datum spelling; verify the actual TM parameters.
            if (source_crs.coordinate_operation.method_name != expected_crs.coordinate_operation.method_name
                or [p.value for p in source_crs.coordinate_operation.params] != [p.value for p in expected_crs.coordinate_operation.params]
                or source_crs.ellipsoid.semi_major_metre != expected_crs.ellipsoid.semi_major_metre
                or source_crs.ellipsoid.inverse_flattening != expected_crs.ellipsoid.inverse_flattening
                or any(axis.unit_conversion_factor != 1 for axis in source_crs.axis_info)):
                raise ValueError('Unexpected source CRS: ' + stem)
            with archive.open(name) as shp, archive.open(stem + '.shx') as shx, archive.open(stem + '.dbf') as dbf:
                reader = shapefile.Reader(shp=shp, shx=shx)
                selected = list(reader.iterShapes(bbox=bounds))
                header = dbf.read(32)
                header_size, row_size = struct.unpack_from('<HH', header, 8)
                fields = []
                for _ in range((header_size - 33) // 32):
                    field = dbf.read(32)
                    fields.append((field[:11].split(b'\0')[0].decode(), field[16]))
                for item in selected:
                    dbf.seek(header_size + item.oid * row_size)
                    row = dbf.read(row_size)
                    if row[:1] == b'*':
                        continue
                    properties, offset = {}, 1
                    for key, length in fields:
                        properties[key] = row[offset:offset + length].decode('cp949').strip(' \x00')
                        offset += length
                    geometry = transform(TO_GEO, shape(item.__geo_interface__))
                    features.append({'type': 'Feature', 'id': stem + ':' + str(item.oid),
                        'properties': properties, 'geometry': geometry.__geo_interface__})
                print(f'{stem}: {len(selected)} regional shapes', flush=True)
    return features


def official_features(raw):
    output = []
    area = box(*CONFIG['bbox'])
    for item in raw:
        p = item['properties']
        geometry = shape(item['geometry'])
        if not geometry.is_valid or geometry.is_empty:
            raise ValueError('Invalid source polygon: ' + item['id'])
        polygons = list(geometry.geoms) if geometry.geom_type == 'MultiPolygon' else [geometry]
        for part, polygon in enumerate(polygons):
            if not area.covers(polygon.centroid):
                continue
            tags = {'building': CONFIG['useTypes'].get(p.get('A9'), 'yes'),
                'source': 'VWorld', 'source:geometry': 'VWorld',
                'source:id': p.get('A1') or item['id'], 'source:record': item['id'],
                'source:date': p.get('A22', ''), 'source:use': p.get('A9', '')}
            if p.get('A24'):
                tags['name'] = p['A24']
            for field, tag, maximum in [('A16', 'height', 400), ('A26', 'building:levels', 100)]:
                value = positive(p.get(field), maximum)
                if value is not None and (field != 'A26' or value.is_integer()):
                    tags[tag] = f'{value:g}'
                    tags['source:' + tag] = 'VWorld'
            output.append({'type': 'Feature', 'id': 'vworld/' + (p.get('A1') or item['id']) + f'/{part}',
                'properties': tags, 'geometry': polygon.__geo_interface__})
    return output


def select_buildings(official):
    output = []
    official_shapes = [transform(TO_LOCAL, shape(f['geometry'])) for f in official]
    additions = STRtree(official_shapes)
    accepted = set()
    report = {'officialAdded': 0, 'officialOverlapsSkipped': []}
    # Prefer records with actual height/levels over empty duplicate footprints.
    order = sorted(range(len(official)), key=lambda i: (
        'height' not in official[i]['properties'], 'building:levels' not in official[i]['properties'], official[i]['id']))
    for index in order:
        feature, polygon = official[index], official_shapes[index]
        if any(int(i) in accepted and polygon.intersection(official_shapes[i]).area > CONFIG['overlapAreaMeters2']
               for i in additions.query(polygon)):
            report['officialOverlapsSkipped'].append(feature['id'])
            continue
        output.append(copy.deepcopy(feature))
        accepted.add(index)
        report['officialAdded'] += 1
    return output, report


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('archive', type=Path)
    args = parser.parse_args()
    with args.archive.open('rb') as source_file:
        digest = hashlib.file_digest(source_file, 'sha256').hexdigest()
    cache = ROOT / 'data/raw/gwacheon/vworld-extract.geojson'
    cache_hash = cache.with_suffix('.sha256')
    if cache.exists() and cache_hash.exists() and cache_hash.read_text() == digest:
        raw = json.loads(cache.read_text(encoding='utf-8'))['features']
    else:
        raw = extract_archive(args.archive)
        cache.write_text(json.dumps({'type': 'FeatureCollection', 'features': raw}, ensure_ascii=False), encoding='utf-8')
        cache_hash.write_text(digest)
    data, report = select_buildings(official_features(raw))
    counts = collections.Counter('height' if positive(f['properties'].get('height'), 400) else
        'levels' if positive(f['properties'].get('building:levels'), 100) else 'fallback' for f in data)
    manifest = {'bbox': CONFIG['bbox'], 'count': len(data),
        'types': dict(collections.Counter(f['properties']['building'] for f in data)),
        **CONFIG['licensing'],
        'heightTags': sum('height' in f['properties'] for f in data),
        'levelTags': sum('building:levels' in f['properties'] for f in data),
        'heightSources': dict(counts), 'selection': report,
        'official': {'archive': args.archive.name, 'sha256': digest, 'crs': CONFIG['sourceCrs'],
            'url': CONFIG['sourceUrl'], 'license': CONFIG['licensing']['licenses']['VWorld'],
            'rawRegionalFeatures': len(raw), 'fieldDefinitions': 'data/references/gwacheon/vworld-columns.xlsx'},
        'limitations': ['VWorld-only footprints and attributes. Snapshot can contain older footprints.',
            'Overlaps >0.5m² between VWorld records are excluded, preferring measured height or levels.',
            'Missing heights remain explicitly estimated; terrain is still Copernicus DSM, not surveyed ground.']}
    out = ROOT / 'public/buildings'
    (out / 'gwacheon.geojson').write_text(json.dumps({'type': 'FeatureCollection',
        'attribution': '국토교통부 GIS건물통합정보 (VWorld)',
        **CONFIG['licensing'], 'features': data}, ensure_ascii=False, separators=(',', ':')), encoding='utf-8')
    (out / 'metadata.json').write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding='utf-8')
    print(json.dumps({'count': len(data), 'heightSources': dict(counts), 'selection': {k: len(v) if isinstance(v, list) else v for k,v in report.items()}}, indent=2))


if __name__ == '__main__':
    sys.stdout.reconfigure(encoding='utf-8')
    main()