City-Lifeline-Cockpit/tools/city-massing/build-candidate-imagery-plates.py

"""Local TMS imagery audit plates; preserve unmarked image next to surveyed outlines."""
import json, math, hashlib
from pathlib import Path
from PIL import Image, ImageDraw, ImageFont

ROOT = Path('.scratch/cockpit-shell-panel-layout/verify/artifacts/t35')
AUDIT = ROOT / 'large-building-audit/audit.json'
OUT = ROOT / 'large-building-audit/imagery-plates'
SOURCE = Path('D:/GIS/GIS/yuanlingPic3')
ZOOM = 19
N = 2 ** ZOOM

def pixel(lon, lat):
    return ((lon + 180) / 360 * N * 256,
            (1 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2 * N * 256)

def rings(geometry):
    if geometry['type'] == 'Polygon': return geometry['coordinates']
    if geometry['type'] == 'MultiPolygon': return [r for p in geometry['coordinates'] for r in p]
    raise ValueError(geometry['type'])

OUT.mkdir(parents=True, exist_ok=True)
candidates = sorted((b for b in json.loads(AUDIT.read_text())['candidates'] if b['tier']=='urban'), key=lambda b: (-b['areaSqm'], b['id']))[:10]
font = ImageFont.truetype('C:/Windows/Fonts/arial.ttf', 18)
records = []
for b in candidates:
    paths = [list(map(lambda p: pixel(*p[:2]), r)) for r in rings(b['geometry'])]
    points = [p for r in paths for p in r]
    x0 = math.floor(min(p[0] for p in points)) - 70
    y0 = math.floor(min(p[1] for p in points)) - 70
    x1 = math.ceil(max(p[0] for p in points)) + 70
    y1 = math.ceil(max(p[1] for p in points)) + 70
    minx, miny = x0//256, y0//256
    maxx, maxy = (x1-1)//256, (y1-1)//256
    mosaic = Image.new('RGB', ((maxx-minx+1)*256, (maxy-miny+1)*256))
    sources = []
    for x in range(minx, maxx+1):
        for y in range(miny, maxy+1):
            path = SOURCE / str(ZOOM) / str(x) / f'{N-1-y}.png'
            data = path.read_bytes()
            with Image.open(path) as tile:
                assert tile.size == (256,256)
                mosaic.paste(tile.convert('RGB'), ((x-minx)*256,(y-miny)*256))
            sources.append({'path': str(path), 'sha256': hashlib.sha256(data).hexdigest(), 'xyz': [ZOOM,x,y], 'tmsY':N-1-y})
    image = mosaic.crop((x0-minx*256,y0-miny*256,x1-minx*256,y1-miny*256))
    marked = image.copy()
    draw = ImageDraw.Draw(marked)
    for r in paths: draw.line([(x-x0,y-y0) for x,y in r], fill='#00ffff',width=3)
    w,h = image.size
    plate = Image.new('RGB',(w*2+12,h+65),'#101927')
    plate.paste(image,(0,65)); plate.paste(marked,(w+12,65))
    draw = ImageDraw.Draw(plate)
    draw.text((8,5),f"{b['id']} | {b['areaSqm']} m2 | survey relief {b.get('terrainRelief')} m",fill='white',font=font)
    draw.text((8,34),'Satellite source (north up)',fill='white',font=font)
    draw.text((w+20,34),'Same source + cyan JMD survey outline',fill='white',font=font)
    path = OUT / f"{b['id']}.png"
    plate.save(path)
    records.append({'id':b['id'],'areaSqm':b['areaSqm'],'terrainRelief':b.get('terrainRelief'),'plate':str(path),'pixelBBox':[x0,y0,x1,y1],'sourceTiles':sources,'geometrySha256':hashlib.sha256(json.dumps(b['geometry'],sort_keys=True).encode()).hexdigest(),'interpretation':None})
(OUT/'manifest.json').write_text(json.dumps({'sourceAuditSha256':hashlib.sha256(AUDIT.read_bytes()).hexdigest(),'zoom':ZOOM,'candidateCount':len(records),'plates':records,'modelChanged':False,'satelliteCaptureDateKnown':False},indent=2),encoding='utf8')
print(json.dumps({'passed':True,'plates':len(records),'output':str(OUT),'modelChanged':False}))