174 lines
8.8 KiB
Python
174 lines
8.8 KiB
Python
"""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()
|