Files
shanhe-puzzle/scripts/prepare-maps.py
T

122 lines
7.2 KiB
Python

import json, math, pathlib, collections
ROOT=pathlib.Path(__file__).resolve().parent.parent
AREAS=json.load(open(ROOT/'research/areas.json'));LOOKUP={a['adcode']:a for a in AREAS}
def area(r):return abs(sum(r[i][0]*r[(i+1)%len(r)][1]-r[(i+1)%len(r)][0]*r[i][1] for i in range(len(r)))/2)
def dist(p,a,b):
dx,dy=b[0]-a[0],b[1]-a[1];n=dx*dx+dy*dy
t=max(0,min(1,((p[0]-a[0])*dx+(p[1]-a[1])*dy)/n)) if n else 0
return (p[0]-a[0]-t*dx)**2+(p[1]-a[1]-t*dy)**2
def dp(p,eps):
if len(p)<3:return p
v,i=max((dist(p[k],p[0],p[-1]),k) for k in range(1,len(p)-1))
if v<=eps*eps:return [p[0],p[-1]]
return dp(p[:i+1],eps)[:-1]+dp(p[i:],eps)
def project(p):return [p[0]*.82,-math.log(math.tan(math.pi/4+math.radians(p[1])/2))*180/math.pi*.82]
def make_level(features,code,name,province):
fs=[f for f in features if f.get('properties',{}).get('name') and f['properties'].get('adcode')]
if len(fs)<2:return None
reg=[];boundary=[]
for f in fs:
prop=f['properties'];g=f['geometry'];polys=g['coordinates'] if g['type']=='MultiPolygon' else [g['coordinates']]
rings=[project_ring(p[0]) for p in polys if len(p[0])>=4]
rings.sort(key=area,reverse=True)
if not rings:continue
main=rings[0];boundary.append({(round(x,3),round(y,3)) for r in rings for x,y in r})
# Keep island geometry, simplify each independently; no invented boundaries.
span=max(max(p[0] for p in main)-min(p[0] for p in main),max(p[1] for p in main)-min(p[1] for p in main))
simple=[]
for r in rings:
localspan=max(max(p[0] for p in r)-min(p[0] for p in r),max(p[1] for p in r)-min(p[1] for p in r))
eps=max(min(span*.005,localspan*.01),.00000001);q=dp(r,eps)
while len(q)>90:eps*=1.4;q=dp(r,eps)
if len(q)<4:q=r
simple.append(q[:-1] if q[0]==q[-1] else q)
cx=sum(x for x,y in main[:-1])/max(1,len(main)-1);cy=sum(y for x,y in main[:-1])/max(1,len(main)-1)
reg.append({'id':str(prop['adcode']),'name':prop['name'],'rings':simple,'center':[cx,cy]})
if len(reg)<2:return None
edges=[]
for i in range(len(reg)):
for j in range(i):
if len(boundary[i]&boundary[j])>=2:edges.append([j,i,0])
# Join disconnected island components with explicit game-only sea routes.
def components():
groups=[{i} for i in range(len(reg))]
for a,b,_ in edges:
ga=next(g for g in groups if a in g);gb=next(g for g in groups if b in g)
if ga is not gb:ga.update(gb);groups.remove(gb)
return groups
while len(components())>1:
gs=components();candidates=[]
for a in gs[0]:
for g in gs[1:]:
for b in g:
pa,pb=reg[a]['center'],reg[b]['center'];candidates.append(((pa[0]-pb[0])**2+(pa[1]-pb[1])**2,a,b))
_,a,b=min(candidates);edges.append([a,b,1])
allp=[p for r in reg for ring in r['rings'] for p in ring]
minx,miny=min(p[0] for p in allp),min(p[1] for p in allp);maxx,maxy=max(p[0] for p in allp),max(p[1] for p in allp)
span=max(maxx-minx,maxy-miny);scale=1000/span
for r in reg:
r['center']=[round((r['center'][0]-minx)*scale,2),round((r['center'][1]-miny)*scale,2)]
r['rings']=[[[round((x-minx)*scale,4),round((y-miny)*scale,4)] for x,y in ring] for ring in r['rings']]
return {'id':str(code),'name':name,'province':str(province),'width':round((maxx-minx)*scale,2),'height':round((maxy-miny)*scale,2),'regions':reg,'edges':edges,'projection':[minx,miny,scale]}
def project_ring(r):return [project(p) for p in r]
levels=[]
cn=json.load(open(ROOT/'research/china.json'));nation=make_level(cn['features'],100000,'中国 · 省级拼图',100000);levels.append(nation)
for p in sorted((ROOT/'data/raw').glob('*.json')):
code=int(p.stem);info=LOOKUP.get(code)
if not info:continue
fs=json.load(open(p)).get('features',[])
# Do not turn city-level outlines into county puzzles.
fs=[f for f in fs if f.get('properties',{}).get('level')=='district']
prov=code if info['level']=='province' else info['parent']
level=make_level(fs,code,info['name'],prov)
if level:levels.append(level)
provinces=[{'id':str(a['adcode']),'name':a['name']} for a in AREAS if a['level']=='province']
result={'source':'DataV.GeoAtlas','downloadDate':'2026-10-16','provinces':provinces,'levels':levels}
(ROOT/'data/county-materials').mkdir(parents=True,exist_ok=True)
(ROOT/'data/county-materials/maps.js').write_text('window.MAP_DATA='+json.dumps(result,ensure_ascii=False,separators=(',',':'))+';\n')
summary={'levels':len(levels),'countyLevels':len(levels)-1,'provincePieces':len(nation['regions']),'countyPieces':sum(len(l['regions']) for l in levels[1:]),'availableProvinces':len({l['province'] for l in levels[1:]}),'seaRoutes':sum(e[2] for l in levels for e in l['edges'])}
(ROOT/'data/county-materials/summary.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2));print(summary)
# Formal release: only province and prefecture-level units, no county pieces.
citylevels=[nation]
excluded=[]
for p in sorted((ROOT/'data/provinces-raw').glob('*.json')):
code=int(p.stem);info=LOOKUP.get(code)
if not info:continue
fs=json.load(open(p)).get('features',[])
cityfs=[f for f in fs if f.get('properties',{}).get('level')=='city' and int(f['properties'].get('adcode',0))%100==0]
excluded.extend(f['properties']['name'] for f in fs if f not in cityfs)
lv=make_level(cityfs,code,info['name'],code)
if lv:
# Missing county-only zones may disconnect prefectures: these are not land neighbors.
for edge in lv['edges']:
if edge[2]:edge[2]=2
citylevels.append(lv)
def inside(p,ring):
x,y=p;hit=False
for i,a in enumerate(ring):
b=ring[(i+1)%len(ring)]
if (a[1]>y)!=(b[1]>y) and x<(b[0]-a[0])*(y-a[1])/(b[1]-a[1])+a[0]:hit=not hit
return hit
for lv in citylevels:
for r in lv['regions']:
ring=r['rings'][0];minx,maxx=min(p[0] for p in ring),max(p[0] for p in ring);miny,maxy=min(p[1] for p in ring),max(p[1] for p in ring)
candidates=[[minx+(maxx-minx)*i/14,miny+(maxy-miny)*j/14] for i in range(1,14) for j in range(1,14)]
candidates=[p for p in candidates if inside(p,ring)]
if candidates:r['label']=max(candidates,key=lambda p:min(dist(p,ring[k],ring[(k+1)%len(ring)]) for k in range(len(ring))))
else:r['label']=r['center']
r['label']=[round(v,3) for v in r['label']]
# Preserve the source's unnamed maritime boundary features as a separate layer.
minx,miny,scale=nation['projection'];decorations=[]
for feature in cn['features']:
if feature['properties'].get('adcode')=='100000_JD':
geom=feature['geometry'];polys=geom['coordinates'] if geom['type']=='MultiPolygon' else [geom['coordinates']]
for poly in polys:decorations.append([[round((project(p)[0]-minx)*scale,4),round((project(p)[1]-miny)*scale,4)] for p in poly[0]])
nation['decorations']=decorations
for lv in citylevels:lv.pop('projection',None)
result={'source':'DataV.GeoAtlas','downloadDate':'2026-10-16','provinces':provinces,'levels':citylevels}
(ROOT/'src/maps.js').write_text('window.MAP_DATA='+json.dumps(result,ensure_ascii=False,separators=(',',':'))+';\n')
summary={'levels':len(citylevels),'prefectureLevels':len(citylevels)-1,'provincePieces':len(nation['regions']),'prefecturePieces':sum(len(l['regions']) for l in citylevels[1:]),'availableProvinces':len(citylevels)-1,'excludedCountyUnits':excluded,'countyMaterials':{'levels':len(levels)-1,'pieces':sum(len(l['regions']) for l in levels[1:]),'path':'data/county-materials/maps.js'}}
(ROOT/'data/map-summary.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2));print('Release:',{k:v for k,v in summary.items() if k!='excludedCountyUnits'})