"""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}))