"""Attach OSM attributes to actual GSI terrain, retaining source identifiers."""
from pathlib import Path
import numpy as np,json,gzip,struct,shapely
from shapely.geometry import Polygon,LineString,box
from shapely.ops import transform
from pyproj import Transformer
from scipy.ndimage import distance_transform_edt,map_coordinates
R=Path(__file__).resolve().parents[1];out=R/'site/public/models';raw=json.loads((R/'data/learning/osm.json').read_text())
for name,stride in [('tokyo',3),('tokyo-detail',1)]:
 a=np.load(R/'data'/name/'dem.npz');xs=a['x_m'][::stride];ys=a['y_m'][::stride];h=a['elevation_m'][::stride,::stride];mask=a['mask'][::stride,::stride];step=float(xs[1]-xs[0]);ny,nx=h.shape
 tx=Transformer.from_crs('EPSG:4326',str(a['crs']),always_xy=True);ll=Transformer.from_crs(str(a['crs']),'EPSG:4326',always_xy=True)
 cx=(xs[0]+xs[-1])/2;cy=(ys[0]+ys[-1])/2;scale=100/max(xs[-1]-xs[0],ys[-1]-ys[0]);clip=box(xs[0],ys[0],xs[-1],ys[-1]);grid=np.zeros(h.shape,np.uint32);features=[{'name':'この格子に対象属性なし','category':'unknown'}]
 nearest=distance_transform_edt(~mask,return_distances=False,return_indices=True);filled=h[tuple(nearest)]
 def point(x,y):
  elev=float(map_coordinates(filled,[[(y-ys[0])/step],[(x-xs[0])/step]],order=1,mode='nearest')[0])
  return [round((x-cx)*scale,5),round(elev*scale,6),round(-(y-cy)*scale,5)]
 elements=raw['elements']
 # Broad areas first; smaller parks/woods override broader residential polygons.
 shapes=[]
 for el in elements:
  tags=el.get('tags',{});coords=[(v['lon'],v['lat']) for v in el.get('geometry',[])];water=tags.get('waterway') in ['river','stream','canal']
  if len(coords)<(2 if water else 4) or (not water and coords[0]!=coords[-1]):continue
  geo=transform(tx.transform,LineString(coords) if water else Polygon(coords))
  if not geo.is_valid:geo=shapely.make_valid(geo)
  geo=geo.intersection(clip)
  if geo.is_empty:continue
  category='river' if water else 'forest' if tags.get('landuse')=='forest' or tags.get('natural')=='wood' else 'park' if tags.get('leisure')=='park' else 'farm' if tags.get('landuse') in ['farmland','orchard'] else 'urban'
  shapes.append((geo,el,category))
 shapes.sort(key=lambda t:-t[0].area)
 for geo,el,category in shapes:
  tags=el['tags'];fid=len(features);label={'river':'名称未登録の水路','forest':'名称未登録の森林','park':'名称未登録の公園','farm':'名称未登録の農地','urban':'名称未登録の市街地'}[category]
  info={'id':el['id'],'name':tags.get('name:ja',tags.get('name',label)),'category':category,'updated':el.get('timestamp',''),'tags':{k:v for k,v in tags.items() if k in ['name','name:ja','waterway','landuse','natural','leisure','leaf_type','leaf_cycle','intermittent','tunnel']}}
  center=geo.representative_point();info['center']=point(center.x,center.y)
  if category=='river':
   parts=list(geo.geoms) if hasattr(geo,'geoms') else [geo];lines=[]
   for part in parts:
    if part.geom_type!='LineString' or part.length<1:continue
    # Surface interpolation is for display only, never a measured water level.
    p=shapely.segmentize(part,min(step,60));lines.append([point(x,y) for x,y in p.coords])
   info['lines']=lines
  else:
   xmin,ymin,xmax,ymax=geo.bounds;ix0=max(0,int(np.floor((xmin-xs[0])/step)));ix1=min(nx,int(np.ceil((xmax-xs[0])/step))+1);iy0=max(0,int(np.floor((ymin-ys[0])/step)));iy1=min(ny,int(np.ceil((ymax-ys[0])/step))+1)
   xx,yy=np.meshgrid(xs[ix0:ix1],ys[iy0:iy1]);inside=shapely.contains_xy(geo,xx,yy);grid[iy0:iy1,ix0:ix1][inside]=fid
  features.append(info)
 valid=mask[:-1,:-1]&mask[:-1,1:]&mask[1:,:-1]&mask[1:,1:]
 for j in range(stride+1):
  for i in range(stride+1):valid &= a['mask'][j:j+(ny-1)*stride:stride,i:i+(nx-1)*stride:stride]
 rr,cc=np.where(valid);q=np.stack([rr*nx+cc,rr*nx+cc+1,(rr+1)*nx+cc+1,(rr+1)*nx+cc],1);tri=np.concatenate([q[:,[0,1,2]],q[:,[0,2,3]]]);used=np.unique(tri);ids=np.full(h.size,-1,dtype=np.int32);ids[used]=np.arange(len(used));ix=used%nx;iy=used//nx
 pos=np.stack([(xs[ix]-cx)*scale,h.ravel()[used]*scale,-(ys[iy]-cy)*scale],1).astype('<f4');indices=ids[tri].astype('<u4').ravel();attrs=grid.ravel()[used].astype('<u4');lon,lat=ll.transform(xs[ix],ys[iy]);geo=np.stack([lon,lat],1).astype('<f4')
 meta={'name':name,'vertices':len(pos),'indices':len(indices),'metersPerUnit':1/scale,'step_m':step,'exaggeration':15,'source':'GSI DEM10B','attributes':True,'coordinateOrder':'longitude,latitude','heightNote':'DEM10B interpolated grid; not survey-grade point measurement'}
 header=json.dumps(meta).encode();header+=b' '*((-len(header))%4);payload=struct.pack('<I',len(header))+header+pos.tobytes()+indices.tobytes()+attrs.tobytes()+geo.tobytes()
 (out/f'{name}.bin.gz').write_bytes(gzip.compress(payload,compresslevel=9,mtime=0));(out/f'{name}-features.json').write_text(json.dumps({'features':features,'snapshot':raw['osm3s']['timestamp_osm_base'],'license':'ODbL-1.0','attribution':'© OpenStreetMap contributors','scope':'Selected tagged ways only; relations omitted; no tag is not absence on the ground','bounds':[139.3,35.52,139.95,35.85]},ensure_ascii=False,separators=(',',':')))
 print(name,'vertices',len(pos),'features',len(features)-1,'categorized vertices',int((attrs>0).sum()),'MB',round((out/f'{name}.bin.gz').stat().st_size/1e6,2),flush=True)
# Provide the source OSM snapshot under ODbL with query and transformation scripts.
import zipfile
with zipfile.ZipFile(R/'site/public/osm-source.zip','w',zipfile.ZIP_DEFLATED) as z:
 z.write(R/'data/learning/osm.json','osm.json');z.write(R/'data/learning/query.overpassql','query.overpassql');z.write(__file__,'export_learning.py');z.writestr('LICENSE.txt','© OpenStreetMap contributors. Database: Open Database License 1.0. https://opendatacommons.org/licenses/odbl/1-0/\nFetched via Overpass. Only selected ways, not a complete map. Attribute rasterization and coordinate transformations are reproducible using export_learning.py and the GSI DEM scripts. GSI elevations are a separate source, not OSM observations.')
