投影之间 · 算法技术文档
18 METHODS / FORWARD & INVERSE / FORMAL MATHEMATICS
回到文档顶部离线公式版 · 2026-10-19

中国与全球地图 18 种投影教学实验：坐标转换算法技术文档

版本：2026-10-19 · 全面核查修复版 v4。本文面向本项目的 HTML + PHP 教学实现，说明坐标模型、18 种投影正反算、数值求解与教学解释。核心实现以 src/Projection.php 为准。

阅读提示：本版的分式、求和、上下标和矩阵使用 LaTeX 数学公式；推荐打开同目录的《投影算法技术文档_公式版.html》查看完整排版。18种投影分别列出正算与反算关系，数值求逆明确给出迭代方程。

1. 坐标转换的范围

本项目将经纬度转换为地图平面坐标，也将平面坐标反算为经纬度。这里的“转换”是投影运算，不包含 WGS84、CGCS2000、北京 54、西安 80 之间的基准转换，不包含 GCJ-02 或 BD-09 的偏移处理，不转换高程、历元或垂直基准。

输入经纬度按十进制度解释，轴序固定为 [经度, 纬度]，输出轴序固定为 [东向 X, 北向 Y]，单位为米。我国测量教材和部分测绘成果使用“X 为北坐标、Y 为东坐标”的约定；导入这些数据前必须交换到本网站的顺序。

1.1 地球模型

方法 本项目计算模型 常数
高斯－克吕格 tmerc WGS84 旋转椭球 长半轴 a = 6378137 m，扁率 f = 1 / 298.257223563，第一偏心率平方 e² = f(2−f)
其余 17 种 统一教学球面 R = 6371007.180918 m

17 种球面算法直接使用输入纬度作为球面纬度，不先进行等积纬度或保角纬度转换。因此，选择了接近 WGS84 等积球的半径，并不意味着已实现对应的 WGS84 椭球投影。球面结果与椭球结果的差异属于模型差异，不能解释成代码计算误差。尤其不能把本项目 merc 输出当作 EPSG:3857：Web Mercator 使用半径 6378137 m，本项目使用上述统一球半径。

高斯－克吕格的椭球正形性质定义在 WGS84 椭球上；其余方法的正形、等积、等距性质定义在上述球面上。跨模型比较时应先说明这一点。

1.2 通用符号与运算约定

符号 含义
λ, φ 输入经度、纬度，进入公式前换算为弧度
λ₀, φ₀ 中央经线与纬度原点，或方位投影中心
φ₁, φ₂ 第一、第二标准纬线
L 归一化经差 wrap(λ−λ₀)，通常落在 [−π, π]
x₀, y₀ 假东距、假北距，米
k₀ 中央比例因子，无量纲，必须大于 0
atan2(y,x) 能区分象限的反正切
hypot(x,y) √(x²+y²)

所有三角函数均使用弧度。输出前对公式结果加假东距、假北距；反算前先扣除它们。下文除高斯－克吕格外，u=(X−x₀)/R、v=(Y−y₀)/R 表示无量纲平面坐标；球面立体另外乘 k₀，所以其反算归一化应除以 Rk₀。本项目只有高斯－克吕格与球面立体使用输入的 k0；墨卡托的尺度由标准纬线 lat1 决定。

下列球面公式先算无量纲坐标，再恢复米制坐标；经纬度必须先从度换成弧度：

λrad=π180λdeg,φrad=π180φdeg,L=wrap⁡(λ−λ0),u=X−x0R,v=Y−y0R,X=x0+Ru,Y=y0+Rv.\begin{aligned} \lambda_{\mathrm{rad}}&=\frac{\pi}{180}\lambda_{\mathrm{deg}}, &\varphi_{\mathrm{rad}}&=\frac{\pi}{180}\varphi_{\mathrm{deg}},\\ L&=\operatorname{wrap}(\lambda-\lambda_0), &u&=\frac{X-x_0}{R},\quad v=\frac{Y-y_0}{R},\\ X&=x_0+Ru, &Y&=y_0+Rv. \end{aligned}(1)

球面立体投影使用 Rk₀ 替代归一化与恢复公式中的 R；高斯－克吕格直接使用其椭球米制公式。

当前程序的统一输入约束比某些方法的数学定义域更窄：经度 [-180°,180°]，输入纬度 (-89.999999°,89.999999°)，参数纬度 (-89.9°,89.9°)，中央经线 [-180°,180°]，0<k₀≤10，假东距、假北距绝对值不大于 10⁹ m，正算输出及反算输入X、Y绝对值均不大于 10¹² m。即使某种投影数学上可表达极点，本网站当前版本仍统一拒绝精确极点；以下极点讨论解释几何极限，不表示界面接受极点输入。

逆三角函数只对浮点舍入引起的微小超限做夹取，例如把 1+10⁻¹⁵ 限制为 1；明显超出定义域的输入应报错，不能夹成一个貌似有效的结果。经度在极点、方位投影的对跖点等位置可能不唯一，这属于投影本身的几何性质。

2. 方法目录与固定顺序

序号 中文名 标识 本项目主要性质
1 墨卡托 merc 球面正形
2 高斯－克吕格 tmerc 椭球正形，分带使用
3 兰伯特正形圆锥 lcc 球面正形
4 阿尔伯斯等积圆锥 aea 球面等积
5 正射方位 ortho 球面透视，半球可见
6 方位等距 aeqd 保持中心至其他点的球面距离
7 等距圆柱 eqc 沿经线与标准纬线保距
8 圆柱等积 cea 球面等积
9 米勒圆柱 mill 折衷
10 等距圆锥 eqdc 沿经线与标准纬线保距
11 兰伯特方位等积 laea 球面等积
12 球面立体／体视 stere 球面正形
13 球心／日晷 gnom 球面大圆映为直线
14 摩尔威德 moll 球面等积
15 正弦 sinu 球面等积
16 等积地球 eqearth 球面等积
17 罗宾森 robin 折衷，系数插值
18 温克尔三重 wintri 折衷，组合映射

前六项保持用户最初提出的顺序。分类用于解释数学性质，不代表所有方法都适合作为中国全国地图的生产投影。USGS 投影手册给出了常见投影的正反算与应用背景。

3. 18 种投影的正算与反算

3.1 墨卡托 Mercator（merc）

本项目采用球面公式，以第一标准纬线确定尺度 k=cosφ₁：

正算与反算公式

k=cos⁡φ1,正算:u=kL,v=kln⁡tan⁡(π4+φ2),反算:λ=λ0+uk,φ=2arctan⁡ ⁣(ev/k)−π2.\begin{aligned} k&=\cos\varphi_1,\\ \text{正算:}\qquad u&=kL,& v&=k\ln\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right),\\ \text{反算:}\qquad \lambda&=\lambda_0+\frac{u}{k},& \varphi&=2\arctan\!\left(e^{v/k}\right)-\frac{\pi}{2}. \end{aligned}(2)

要求 |φ| < π/2；两极趋于无穷远。它保持无穷小图形的角度，局部比例因子为 k/cosφ，面积比例为其平方。正形不等于全国轮廓完全不变，也不代表地图任意距离正确。PROJ：Mercator

3.2 高斯－克吕格 Gauss–Krüger（tmerc）

高斯－克吕格是横轴墨卡托的一种常见使用形式，一般 k₀=1。本实现采用 WGS84 椭球及六阶第三扁率 Krüger 级数，避免用只在很窄经差内表现较好的低阶经差展开覆盖中国全境。

正算与反算公式

n=f2−f,e=f(2−f),A=a1+n(1+n24+n464+n6256),χ=arctan⁡ ⁣{sinh⁡ ⁣[arsinh⁡(tan⁡φ)−eartanh⁡(esin⁡φ)]},ξ′=atan2⁡(tan⁡χ,cos⁡L),η′=arsinh⁡ ⁣(sin⁡Ltan⁡2χ+cos⁡2L),ξ=ξ′+∑j=16αjsin⁡(2jξ′)cosh⁡(2jη′),η=η′+∑j=16αjcos⁡(2jξ′)sinh⁡(2jη′),X=x0+k0Aη,Y=y0+k0A(ξ−ξ0).\begin{aligned} n&=\frac{f}{2-f},&e&=\sqrt{f(2-f)},\\ A&=\frac{a}{1+n}\left(1+\frac{n^2}{4}+\frac{n^4}{64}+\frac{n^6}{256}\right),\\ \chi&=\arctan\!\left\{\sinh\!\left[\operatorname{arsinh}(\tan\varphi) -e\operatorname{artanh}(e\sin\varphi)\right]\right\},\\ \xi'&=\operatorname{atan2}(\tan\chi,\cos L),\\ \eta'&=\operatorname{arsinh}\!\left(\frac{\sin L}{\sqrt{\tan^2\chi+\cos^2 L}}\right),\\ \xi&=\xi'+\sum_{j=1}^{6}\alpha_j\sin(2j\xi')\cosh(2j\eta'),\\ \eta&=\eta'+\sum_{j=1}^{6}\alpha_j\cos(2j\xi')\sinh(2j\eta'),\\ X&=x_0+k_0A\eta,\\ Y&=y_0+k_0A(\xi-\xi_0). \end{aligned}(3)

α1=n2−2n23+5n316+41n4180−127n5288+7891n637800,α2=13n248−3n35+557n41440+281n5630−1983433n61935360,α3=61n3240−103n4140+15061n526880+167603n6181440,α4=49561n4161280−179n5168+6601661n67257600,α5=34729n580640−3418889n61995840,α6=212378941n6319334400.\begin{aligned} \alpha_1={}&\frac n2-\frac{2n^2}{3}+\frac{5n^3}{16}+\frac{41n^4}{180} -\frac{127n^5}{288}+\frac{7891n^6}{37800},\\ \alpha_2={}&\frac{13n^2}{48}-\frac{3n^3}{5}+\frac{557n^4}{1440} +\frac{281n^5}{630}-\frac{1983433n^6}{1935360},\\ \alpha_3={}&\frac{61n^3}{240}-\frac{103n^4}{140} +\frac{15061n^5}{26880}+\frac{167603n^6}{181440},\\ \alpha_4={}&\frac{49561n^4}{161280}-\frac{179n^5}{168} +\frac{6601661n^6}{7257600},\\ \alpha_5={}&\frac{34729n^5}{80640}-\frac{3418889n^6}{1995840},\\ \alpha_6={}&\frac{212378941n^6}{319334400}. \end{aligned}(4)

反算（六阶逆级数，与当前代码一致）

设带号为 N，带号前缀开关 s 取0或1；扩展演示模式固定 s=0。先扣除有效假东距，再反解辅助坐标：

x0∗=x0+sN 106,η=X−x0∗k0A,ξ=Y−y0k0A+ξ0,ξ′=ξ−∑j=16βjsin⁡(2jξ)cosh⁡(2jη),η′=η−∑j=16βjcos⁡(2jξ)sinh⁡(2jη),L=atan2⁡(sinh⁡η′,cos⁡ξ′),χ=atan2⁡ ⁣(sin⁡ξ′,sinh⁡2η′+cos⁡2ξ′).\begin{aligned} x_0^*&=x_0+sN\,10^6,& \eta&=\frac{X-x_0^*}{k_0A},& \xi&=\frac{Y-y_0}{k_0A}+\xi_0,\\ \xi'&=\xi-\sum_{j=1}^{6}\beta_j\sin(2j\xi)\cosh(2j\eta),& \eta'&=\eta-\sum_{j=1}^{6}\beta_j\cos(2j\xi)\sinh(2j\eta),\\ L&=\operatorname{atan2}(\sinh\eta',\cos\xi'),& \chi&=\operatorname{atan2}\!\left(\sin\xi',\sqrt{\sinh^2\eta'+\cos^2\xi'}\right). \end{aligned}(5)

β1=n2−2n23+37n396−n4360−81n5512+96199n6604800,β2=n248+n315−437n41440+46n5105−1118711n63870720,β3=17n3480−37n4840−209n54480+5569n690720,β4=4397n4161280−11n5504−830251n67257600,β5=4583n5161280−108847n63991680,β6=20648693n6638668800.\begin{aligned} \beta_1={}&\frac n2-\frac{2n^2}{3}+\frac{37n^3}{96}-\frac{n^4}{360}-\frac{81n^5}{512}+\frac{96199n^6}{604800},\\ \beta_2={}&\frac{n^2}{48}+\frac{n^3}{15}-\frac{437n^4}{1440}+\frac{46n^5}{105}-\frac{1118711n^6}{3870720},\\ \beta_3={}&\frac{17n^3}{480}-\frac{37n^4}{840}-\frac{209n^5}{4480}+\frac{5569n^6}{90720},\\ \beta_4={}&\frac{4397n^4}{161280}-\frac{11n^5}{504}-\frac{830251n^6}{7257600},\\ \beta_5={}&\frac{4583n^5}{161280}-\frac{108847n^6}{3991680},\\ \beta_6={}&\frac{20648693n^6}{638668800}. \end{aligned}(6)

最后由保角纬度恢复大地纬度。令 q=arsinh(tanχ)，从 φ₀=χ 开始作一维牛顿迭代：

g(φ)=arsinh⁡(tan⁡φ)−eartanh⁡(esin⁡φ),g′(φ)=1−e2cos⁡φ(1−e2sin⁡2φ),φi+1=φi−g(φi)−qg′(φi),λ=wrap⁡(λ0+L).\begin{aligned} g(\varphi)&=\operatorname{arsinh}(\tan\varphi)-e\operatorname{artanh}(e\sin\varphi),\\ g'(\varphi)&=\frac{1-e^2}{\cos\varphi(1-e^2\sin^2\varphi)},\\ \varphi_{i+1}&=\varphi_i-\frac{g(\varphi_i)-q}{g'(\varphi_i)},& \lambda&=\operatorname{wrap}(\lambda_0+L). \end{aligned}(7)

最多迭代12次，纬度改变量小于 10⁻¹⁴ rad 时停止，随后执行统一的正算回代。这样不再需要在经差边界外取二维差分，也不再受旧版 |φ|<1.56 rad 候选纬度限制。极区经度恢复会放大浮点舍入；代码只允许与 cosχ 条件数相匹配的微小边界夹取，之后仍需通过米制回代。

ξ₀ 由同样算法在 φ=φ₀, L=0 求得；六个级数系数已完整列于上式。

上式采用代码中的保角纬度 χ 表达，与常见的辅助正切纬度表达等价。正反算采用六阶 Krüger 系数；保角纬度反求大地纬度采用上式一维牛顿法。与PROJ的系数阶数一致，不表示其内部全部数值代码相同。

界面提供三种模式：扩展宽带演示（中央经线±45°）、3°分带（±1.5°）、6°分带（±3°）。正算、反算、地图预览与导出均使用同一分带范围。宽带是观察变形的教学窗口，不是正式测绘90°分带，更不是完整全球投影。

3°带中央经线采用 3N°，6°带为 (6N−3)°。zone_width=0/3/6；zone_number=0 表示按中央经线推算，正整数则由带号设置中央经线。带号按东经0°至360°编号，3°带支持1～120（中央经线0°对应120带），6°带支持1～60；中央经线归一化为[-180°,180°)。这是本教学程序约定，不能自动解释各国不同带号体系。

zone_prefix=1 时正算东坐标再加 N×10⁶ m，反算扣除同一前缀；x0仍是用户设置的假东距，常用500000m。正算上式X应使用有效假东距 x₀*。例如郑州点113.6254°E、34.7466°N：3°第38带、中央经线114°、k₀=1、x₀=500000m、y₀=0时，东坐标465698.278m，北坐标3846545.058m；开启前缀后东坐标为38465698.278m。

内置中国分带图会重建带内闭合多边形；上传数据严格导出遇到域外顶点仍报错。本功能没有自动识别或转换CGCS2000、北京54、西安80等基准。全国单中央经线图不能作为全国统一的常规分带测量坐标。

参考：PROJ：Transverse Mercator、Karney 的原始论文与算法资料。文献精度结论属于相应算法与测试域，不直接作为本网站已经达到的精度承诺。

3.3 兰伯特正形圆锥 Lambert Conformal Conic（lcc）

设 t(φ)=tan(π/4+φ/2)。

正算与反算公式

t(φ)=tan⁡(π4+φ2),n={ln⁡(cos⁡φ1/cos⁡φ2)ln⁡[t(φ2)/t(φ1)],φ1≠φ2,sin⁡φ1,φ1=φ2,F=cos⁡φ1 t(φ1)nn,ρ=Ft(φ)n,ρ0=Ft(φ0)n,θ=nL,u=ρsin⁡θ,v=ρ0−ρcos⁡θ.\begin{aligned} t(\varphi)&=\tan\left(\frac\pi4+\frac\varphi2\right),\\ n&=\begin{cases} \displaystyle\frac{\ln(\cos\varphi_1/\cos\varphi_2)} {\ln[t(\varphi_2)/t(\varphi_1)]},&\varphi_1\ne\varphi_2,\\ \sin\varphi_1,&\varphi_1=\varphi_2, \end{cases}\\ F&=\frac{\cos\varphi_1\,t(\varphi_1)^n}{n},& \rho&=\frac F{t(\varphi)^n},&\rho_0&=\frac F{t(\varphi_0)^n},\\ \theta&=nL,&u&=\rho\sin\theta,&v&=\rho_0-\rho\cos\theta. \end{aligned}(8)

s=sgn⁡(n),ρ=su2+(ρ0−v)2,θ=atan2⁡ ⁣(su,s(ρ0−v)),λ=λ0+θn,φ=2arctan⁡ ⁣[(Fρ)1/n]−π2.\begin{aligned} s&=\operatorname{sgn}(n),&\rho&=s\sqrt{u^2+(\rho_0-v)^2},\\ \theta&=\operatorname{atan2}\!\bigl(su,s(\rho_0-v)\bigr),\\ \lambda&=\lambda_0+\frac\theta n,& \varphi&=2\arctan\!\left[\left(\frac F\rho\right)^{1/n}\right]-\frac\pi2. \end{aligned}(9)

要求标准纬线不在极点、圆锥常数 n≠0，两条标准纬线不能互为相反数。圆锥顶点需要按极点极限处理。它保角，标准纬线附近尺度较合适；不保面积。PROJ：LCC

3.4 阿尔伯斯等积圆锥 Albers Equal Area（aea）

正算与反算公式

n=sin⁡φ1+sin⁡φ22,C=cos⁡2φ1+2nsin⁡φ1,ρ=C−2nsin⁡φn,ρ0=C−2nsin⁡φ0n,θ=nL,u=ρsin⁡θ,v=ρ0−ρcos⁡θ.\begin{aligned} n&=\frac{\sin\varphi_1+\sin\varphi_2}{2},& C&=\cos^2\varphi_1+2n\sin\varphi_1,\\ \rho&=\frac{\sqrt{C-2n\sin\varphi}}{n},& \rho_0&=\frac{\sqrt{C-2n\sin\varphi_0}}{n},\\ \theta&=nL,&u&=\rho\sin\theta,&v&=\rho_0-\rho\cos\theta. \end{aligned}(10)

s=sgn⁡(n),ρ=su2+(ρ0−v)2,θ=atan2⁡ ⁣(su,s(ρ0−v)),λ=λ0+θn,φ=arcsin⁡ ⁣(C−n2ρ22n).\begin{aligned} s&=\operatorname{sgn}(n),&\rho&=s\sqrt{u^2+(\rho_0-v)^2},\\ \theta&=\operatorname{atan2}\!\bigl(su,s(\rho_0-v)\bigr),\\ \lambda&=\lambda_0+\frac\theta n,& \varphi&=\arcsin\!\left(\frac{C-n^2\rho^2}{2n}\right). \end{aligned}(11)

要求 n≠0，根式与反正弦在定义域内。相同标准纬线可以构成切圆锥；关于赤道对称的标准纬线会使本圆锥形式退化。保持面积意味着球面区域与其平面像的面积成统一比例，不意味着每个方向长度都不变。PROJ：Albers

3.5 正射方位 Orthographic（ortho）

方位投影共用下列量：

正算与反算公式

q=sin⁡φ0sin⁡φ+cos⁡φ0cos⁡φcos⁡L=cos⁡c,Ah=cos⁡φsin⁡L,Bh=cos⁡φ0sin⁡φ−sin⁡φ0cos⁡φcos⁡L.\begin{aligned} q&=\sin\varphi_0\sin\varphi+\cos\varphi_0\cos\varphi\cos L=\cos c,\\ A_h&=\cos\varphi\sin L,\\ B_h&=\cos\varphi_0\sin\varphi-\sin\varphi_0\cos\varphi\cos L. \end{aligned}(12)

u=Ah,v=Bh,q≥0.u=A_h,\qquad v=B_h,\qquad q\ge0.(13)

r=u2+v2,c=arcsin⁡r,0≤r≤1.r=\sqrt{u^2+v^2},\qquad c=\arcsin r,\qquad 0\le r\le1.(14)

r=u2+v2,φ=arcsin⁡ ⁣(cos⁡csin⁡φ0+vsin⁡ccos⁡φ0r),λ=λ0+atan2⁡ ⁣(usin⁡c, rcos⁡φ0cos⁡c−vsin⁡φ0sin⁡c).\begin{aligned} r&=\sqrt{u^2+v^2},\\ \varphi&=\arcsin\!\left(\cos c\sin\varphi_0+ \frac{v\sin c\cos\varphi_0}{r}\right),\\ \lambda&=\lambda_0+\operatorname{atan2}\!\left( u\sin c,\ r\cos\varphi_0\cos c-v\sin\varphi_0\sin c\right). \end{aligned}(15)

当 r=0 时，直接取 λ=λ₀、φ=φ₀。方向量 q、A_h、B_h 的定义见正射方位一节。

c 是投影中心到该点的球心角。正射正算为 u=A_h, v=B_h。只绘制 q≥0 的正面半球；背面点不能投到同一张地图冒充可见点。

反算令 r=hypot(u,v)，要求 r≤1，c=asin(r)，再使用第 3.6 节给出的方位通用逆式。r=0 时直接返回中心。地平线是 r=1，边缘附近的反算条件变差。PROJ：Orthographic

3.6 方位等距 Azimuthal Equidistant（aeqd）

用第 3.5 节的 q,A_h,B_h，令 c=acos(q)、K=c/sin(c)。中心处按极限令 K=1。

正算与反算公式

c=arccos⁡q,K={c/sin⁡c,c≠0,1,c=0,u=KAh,v=KBh,反算使用c=r=u2+v2,0≤r<π.\begin{aligned} c&=\arccos q,& K&=\begin{cases}c/\sin c,&c\ne0,\\1,&c=0,\end{cases}\\ u&=KA_h,&v&=KB_h,\\ \text{反算使用}\qquad c&=r=\sqrt{u^2+v^2},&0&\le r<\pi. \end{aligned}(16)

r=u2+v2,φ=arcsin⁡ ⁣(cos⁡csin⁡φ0+vsin⁡ccos⁡φ0r),λ=λ0+atan2⁡ ⁣(usin⁡c, rcos⁡φ0cos⁡c−vsin⁡φ0sin⁡c).\begin{aligned} r&=\sqrt{u^2+v^2},\\ \varphi&=\arcsin\!\left(\cos c\sin\varphi_0+ \frac{v\sin c\cos\varphi_0}{r}\right),\\ \lambda&=\lambda_0+\operatorname{atan2}\!\left( u\sin c,\ r\cos\varphi_0\cos c-v\sin\varphi_0\sin c\right). \end{aligned}(17)

当 r=0 时，直接取 λ=λ₀、φ=φ₀。方向量 q、A_h、B_h 的定义见正射方位一节。

中心处单独处理；精确对跖点 c=π 的方向不唯一，应拒绝。逆算要求 r<π。其性质是 中心到点的平面距离 = R c，不保证任意两个非中心点之间保距。上述逆式也用于后面的其他方位投影，只需更换 c(r)。PROJ：Azimuthal Equidistant

3.7 等距圆柱 Equidistant Cylindrical（eqc）

正算与反算公式

正算:u=Lcos⁡φ1,v=φ−φ0,反算:λ=λ0+ucos⁡φ1,φ=v+φ0.\begin{aligned} \text{正算:}\quad u&=L\cos\varphi_1,&v&=\varphi-\varphi_0,\\ \text{反算:}\quad\lambda&=\lambda_0+\frac u{\cos\varphi_1},& \varphi&=v+\varphi_0. \end{aligned}(18)

要求 cosφ₁≠0。当 φ₁=0 时就是常见的 Plate Carrée。沿经线方向的弧长和标准纬线上的弧长按球半径正确表达，斜向任意点对的直线距离通常错误。它一般不正形、不等积。PROJ：Equidistant Cylindrical

3.8 圆柱等积 Cylindrical Equal Area（cea）

正算与反算公式

k=cos⁡φ1,正算:u=kL,v=sin⁡φk,反算:λ=λ0+uk,φ=arcsin⁡(kv).\begin{aligned} k&=\cos\varphi_1,\\ \text{正算:}\quad u&=kL,&v&=\frac{\sin\varphi}{k},\\ \text{反算:}\quad\lambda&=\lambda_0+\frac uk,& \varphi&=\arcsin(kv). \end{aligned}(19)

要求 k>0，标准纬线不能是极点；反算要求 |kv|≤1。φ₁=0°、30°、45° 分别对应常见的兰伯特圆柱等积、Behrmann 与 Gall–Peters 实例。标准纬线改变形状分配，但保持等积性质。PROJ：CEA

3.9 米勒圆柱 Miller Cylindrical（mill）

正算与反算公式

正算:u=L,v=54ln⁡tan⁡(π4+2φ5),反算:λ=λ0+u,φ=52[arctan⁡ ⁣(e4v/5)−π4].\begin{aligned} \text{正算:}\quad u&=L,& v&=\frac54\ln\tan\left(\frac\pi4+\frac{2\varphi}{5}\right),\\ \text{反算:}\quad\lambda&=\lambda_0+u,& \varphi&=\frac52\left[\arctan\!\left(e^{4v/5}\right)-\frac\pi4\right]. \end{aligned}(20)

可表达两极的有限纵坐标，但不正形、不等积。它对墨卡托的纬度变换做了压缩，使高纬区的拉伸缓和。PROJ：Miller

3.10 等距圆锥 Equidistant Conic（eqdc）

正算与反算公式

n={cos⁡φ1−cos⁡φ2φ2−φ1,φ1≠φ2,sin⁡φ1,φ1=φ2,G=cos⁡φ1n+φ1,ρ=G−φ,ρ0=G−φ0,θ=nL,u=ρsin⁡θ,v=ρ0−ρcos⁡θ.\begin{aligned} n&=\begin{cases} \displaystyle\frac{\cos\varphi_1-\cos\varphi_2}{\varphi_2-\varphi_1}, &\varphi_1\ne\varphi_2,\\ \sin\varphi_1,&\varphi_1=\varphi_2, \end{cases}\\ G&=\frac{\cos\varphi_1}{n}+\varphi_1,& \rho&=G-\varphi,&\rho_0&=G-\varphi_0,\\ \theta&=nL,&u&=\rho\sin\theta,&v&=\rho_0-\rho\cos\theta. \end{aligned}(21)

s=sgn⁡(n),ρ=su2+(ρ0−v)2,θ=atan2⁡ ⁣(su,s(ρ0−v)),λ=λ0+θn,φ=G−ρ.\begin{aligned} s&=\operatorname{sgn}(n),&\rho&=s\sqrt{u^2+(\rho_0-v)^2},\\ \theta&=\operatorname{atan2}\!\bigl(su,s(\rho_0-v)\bigr),\\ \lambda&=\lambda_0+\frac\theta n,&\varphi&=G-\rho. \end{aligned}(22)

要求 n≠0。经线方向尺度为 1，标准纬线方向尺度为 1；通常既不保角也不等积，不保持所有点对的直线距离。PROJ：Equidistant Conic

3.11 兰伯特方位等积 Lambert Azimuthal Equal Area（laea）

使用第 3.5 节的共用量：

正算与反算公式

K=21+q,u=KAh,v=KBh,反算使用r=u2+v2,c=2arcsin⁡r2,0≤r<2.\begin{aligned} K&=\sqrt{\frac2{1+q}},&u&=KA_h,&v&=KB_h,\\ \text{反算使用}\quad r&=\sqrt{u^2+v^2},&c&=2\arcsin\frac r2,&0&\le r<2. \end{aligned}(23)

r=u2+v2,φ=arcsin⁡ ⁣(cos⁡csin⁡φ0+vsin⁡ccos⁡φ0r),λ=λ0+atan2⁡ ⁣(usin⁡c, rcos⁡φ0cos⁡c−vsin⁡φ0sin⁡c).\begin{aligned} r&=\sqrt{u^2+v^2},\\ \varphi&=\arcsin\!\left(\cos c\sin\varphi_0+ \frac{v\sin c\cos\varphi_0}{r}\right),\\ \lambda&=\lambda_0+\operatorname{atan2}\!\left( u\sin c,\ r\cos\varphi_0\cos c-v\sin\varphi_0\sin c\right). \end{aligned}(24)

当 r=0 时，直接取 λ=λ₀、φ=φ₀。方向量 q、A_h、B_h 的定义见正射方位一节。

中心按极限处理。精确对跖点使 1+q=0，不能唯一映射；逆算要求 r<2。它保持球面面积，在中心附近局部形状较好，离中心较远形状变形增大。PROJ：LAEA

3.12 球面立体 Stereographic（stere）

正算与反算公式

K=21+q,u=KAh,v=KBh,X=x0+Rk0u,Y=y0+Rk0v,反算先取u=X−x0Rk0,v=Y−y0Rk0,r=u2+v2,c=2arctan⁡r2.\begin{aligned} K&=\frac2{1+q},&u&=KA_h,&v&=KB_h,\\ X&=x_0+Rk_0u,&Y&=y_0+Rk_0v,\\ \text{反算先取}\quad u&=\frac{X-x_0}{Rk_0},&v&=\frac{Y-y_0}{Rk_0},\\ r&=\sqrt{u^2+v^2},&c&=2\arctan\frac r2. \end{aligned}(25)

r=u2+v2,φ=arcsin⁡ ⁣(cos⁡csin⁡φ0+vsin⁡ccos⁡φ0r),λ=λ0+atan2⁡ ⁣(usin⁡c, rcos⁡φ0cos⁡c−vsin⁡φ0sin⁡c).\begin{aligned} r&=\sqrt{u^2+v^2},\\ \varphi&=\arcsin\!\left(\cos c\sin\varphi_0+ \frac{v\sin c\cos\varphi_0}{r}\right),\\ \lambda&=\lambda_0+\operatorname{atan2}\!\left( u\sin c,\ r\cos\varphi_0\cos c-v\sin\varphi_0\sin c\right). \end{aligned}(26)

当 r=0 时，直接取 λ=λ₀、φ=φ₀。方向量 q、A_h、B_h 的定义见正射方位一节。

精确对跖点发散；本项目使用球面正形形式，而非某个椭球斜轴立体 EPSG CRS。随着距中心增大，线性尺度增大；正形不保面积。若采用方法支持的中心比例因子，需要先在正算乘、反算除该因子。PROJ：Stereographic

3.13 球心 Gnomonic（gnom）

正算与反算公式

K=1q,u=KAh,v=KBh,q>0,反算使用r=u2+v2,c=arctan⁡r.\begin{aligned} K&=\frac1q,&u&=KA_h,&v&=KB_h,&q&>0,\\ \text{反算使用}\quad r&=\sqrt{u^2+v^2},&c&=\arctan r. \end{aligned}(27)

r=u2+v2,φ=arcsin⁡ ⁣(cos⁡csin⁡φ0+vsin⁡ccos⁡φ0r),λ=λ0+atan2⁡ ⁣(usin⁡c, rcos⁡φ0cos⁡c−vsin⁡φ0sin⁡c).\begin{aligned} r&=\sqrt{u^2+v^2},\\ \varphi&=\arcsin\!\left(\cos c\sin\varphi_0+ \frac{v\sin c\cos\varphi_0}{r}\right),\\ \lambda&=\lambda_0+\operatorname{atan2}\!\left( u\sin c,\ r\cos\varphi_0\cos c-v\sin\varphi_0\sin c\right). \end{aligned}(28)

当 r=0 时，直接取 λ=λ₀、φ=φ₀。方向量 q、A_h、B_h 的定义见正射方位一节。

只允许 q>0，即与中心相距小于 90°；接近 90° 时坐标趋于无穷远。球面上大圆弧的可见部分映成直线，这是路径教学的关键性质；不等于椭球测地线在一般球心投影中严格成直线。PROJ：Gnomonic

3.14 摩尔威德 Mollweide（moll）

先求辅助角 θ：

正算与反算公式

2θ+sin⁡2θ=πsin⁡φ,−π2≤θ≤π2,正算:u=22πLcos⁡θ,v=2sin⁡θ,反算:θ=arcsin⁡v2,φ=arcsin⁡ ⁣(2θ+sin⁡2θπ),λ=λ0+πu22cos⁡θ.\begin{aligned} 2\theta+\sin2\theta&=\pi\sin\varphi, &-\frac\pi2&\le\theta\le\frac\pi2,\\ \text{正算:}\quad u&=\frac{2\sqrt2}{\pi}L\cos\theta, &v&=\sqrt2\sin\theta,\\ \text{反算:}\quad\theta&=\arcsin\frac v{\sqrt2}, &\varphi&=\arcsin\!\left(\frac{2\theta+\sin2\theta}{\pi}\right),\\ \lambda&=\lambda_0+\frac{\pi u}{2\sqrt2\cos\theta}. \end{aligned}(29)

θj+1=θj−2θj+sin⁡2θj−πsin⁡φ2+2cos⁡2θj,θ0=φ.\theta_{j+1}=\theta_j- \frac{2\theta_j+\sin2\theta_j-\pi\sin\varphi}{2+2\cos2\theta_j}, \qquad\theta_0=\varphi.(30)

辅助角从 θ=φ 开始牛顿迭代，最多 30 次，步长绝对值小于 10⁻¹³ rad 时停止。精确极点由通用输入检查拒绝；从公式极限看，整条极点纬圈坍缩为一点，经度不唯一。逆算除检查 Y 外，还通过恢复经差的合法性及正算回代检查椭圆外框。PROJ：Mollweide

3.15 正弦 Sinusoidal（sinu）

正算与反算公式

正算:u=Lcos⁡φ,v=φ,反算:φ=v,λ=λ0+ucos⁡φ.\begin{aligned} \text{正算:}\quad u&=L\cos\varphi,&v&=\varphi,\\ \text{反算:}\quad\varphi&=v,&\lambda&=\lambda_0+\frac u{\cos\varphi}. \end{aligned}(31)

两极 cosφ=0，经度不唯一，需单独处理。它保持面积，中央经线与每条纬线的尺度正确；远离中央经线时角度变形明显。“纬线长度正确”不表示跨纬度点对距离都正确。PROJ：Sinusoidal

3.16 等积地球 Equal Earth（eqearth）

采用常数 A₁=1.340264、A₂=−0.081106、A₃=0.000893、A₄=0.003796，M=√3/2。

正算与反算公式

M=32,θ=arcsin⁡(Msin⁡φ),D(θ)=A1+3A2θ2+7A3θ6+9A4θ8,P(θ)=θ(A1+A2θ2+A3θ6+A4θ8),正算:u=Lcos⁡θMD(θ),v=P(θ),反算:P(θ)=v,φ=arcsin⁡sin⁡θM,λ=λ0+MuD(θ)cos⁡θ.\begin{aligned} M&=\frac{\sqrt3}{2},&\theta&=\arcsin(M\sin\varphi),\\ D(\theta)&=A_1+3A_2\theta^2+7A_3\theta^6+9A_4\theta^8,\\ P(\theta)&=\theta(A_1+A_2\theta^2+A_3\theta^6+A_4\theta^8),\\ \text{正算:}\quad u&=\frac{L\cos\theta}{MD(\theta)},&v&=P(\theta),\\ \text{反算:}\quad P(\theta)&=v,&\varphi&=\arcsin\frac{\sin\theta}{M},\\ \lambda&=\lambda_0+\frac{MuD(\theta)}{\cos\theta}. \end{aligned}(32)

θj+1=θj−P(θj)−vD(θj),θ0=clip⁡ ⁣(vA1,−π3,π3).\theta_{j+1}=\theta_j-\frac{P(\theta_j)-v}{D(\theta_j)}, \qquad\theta_0=\operatorname{clip}\!\left(\frac v{A_1},-\frac\pi3,\frac\pi3\right).(33)

反算的有效 θ 位于 [−π/3,π/3]；初值取 v/A₁ 并夹入这个区间。导数就是 D(θ)，最多牛顿迭代 30 次，步长绝对值小于 10⁻¹³ rad 时停止，再检查反三角函数定义域和正算回代。等积地球保面积，不保角度；本项目使用球面版本。PROJ：Equal Earth、PROJ 算法源码

3.17 罗宾森 Robinson（robin）

罗宾森没有一个简单的初等闭式公式。本实现采用按绝对纬度每 5° 一段的三次多项式系数表。

正算与反算公式

d=180π∣φ∣,i=min⁡ ⁣(18,⌊d5+10−12⌋),t=d−5i,Px,i(t)=ax,i+bx,it+cx,it2+dx,it3,Py,i(t)=ay,i+by,it+cy,it2+dy,it3,u=0.8487LPx,i(t),v=1.3523sgn⁡(φ)Py,i(t).\begin{aligned} d&=\frac{180}{\pi}|\varphi|,& i&=\min\!\left(18,\left\lfloor\frac d5+10^{-12}\right\rfloor\right),\\ t&=d-5i,\\ P_{x,i}(t)&=a_{x,i}+b_{x,i}t+c_{x,i}t^2+d_{x,i}t^3,\\ P_{y,i}(t)&=a_{y,i}+b_{y,i}t+c_{y,i}t^2+d_{y,i}t^3,\\ u&=0.8487L P_{x,i}(t),\\ v&=1.3523\operatorname{sgn}(\varphi)P_{y,i}(t). \end{aligned}(34)

H(ψ)=1.3523Py,i(ψ) ⁣(t(ψ))−∣v∣=0,ψ∈[0,π2],mj=ℓj+hj2,(ℓj+1,hj+1)={(mj,hj),H(mj)<0,(ℓj,mj),H(mj)≥0,φ=sgn⁡(v)ℓ50+h502,λ=λ0+u0.8487Px,i(t).\begin{aligned} H(\psi)&=1.3523P_{y,i(\psi)}\!\bigl(t(\psi)\bigr)-|v|=0, &\psi&\in\left[0,\frac\pi2\right],\\ m_j&=\frac{\ell_j+h_j}{2},\\ (\ell_{j+1},h_{j+1})&= \begin{cases}(m_j,h_j),&H(m_j)<0,\\(\ell_j,m_j),&H(m_j)\ge0,\end{cases}\\ \varphi&=\operatorname{sgn}(v)\frac{\ell_{50}+h_{50}}2,& \lambda&=\lambda_0+\frac{u}{0.8487P_{x,i}(t)}. \end{aligned}(35)

反算使用二分区间 ℓ₀=0、h₀=π/2；每次以候选纬度重新确定分段与多项式，参数 t 的单位始终为度。

完整系数表由 Projection::robin() 从 src/robin-table.php 读取；下表随代码生成，修改算法系数后需同步重新生成文档。系数来自 PROJ 的罗宾森球面实现，并保留其单精度系数转换后的数值。它不是把19个原始节点简单连直线；插值形式不同会得到不同结果。

反算先检查 |v|≤1.3523（容差 10⁻⁹），对单调的 Py 在 [0,π/2] 上用 50 次纬度二分搜索，恢复纬度符号后再由 L=u/[0.8487Px] 求经度。纬度分段边界采用一致的归段方式。它既不等积也不正形，适合作为折衷投影对照。PROJ：Robinson、PROJ 系数表与源码

罗宾森分段系数表（两组，各19行）

下表按5°纬度分段，四列依次为三次多项式的常数、一次、二次、三次系数。数值用9位有效数字表示源表的IEEE-754单精度常数；代码中的完整双精度展开值仍保留在 src/robin-table.php。最后一行是90°边界节点。

横向系数 Px

纬度起点 (°) a b c d
0 1 2.21989998e-17 -7.1551498e-05 3.11029999e-06
5 0.998600006 -0.000482242991 -2.48969991e-05 -1.33089998e-06
10 0.995400012 -0.000831030018 -4.48604987e-05 -9.86701025e-07
15 0.99000001 -0.00135363999 -5.96610007e-05 3.67770008e-06
20 0.982200027 -0.00167441997 -4.49546997e-06 -5.72411e-06
25 0.97299999 -0.00214867992 -9.03570981e-05 1.87359994e-08
30 0.959999979 -0.00305085001 -9.00761006e-05 1.64917003e-06
35 0.942700028 -0.00382791995 -6.53386014e-05 -2.61540004e-06
40 0.921599984 -0.00467746006 -0.000104569997 4.81242978e-06
45 0.896200001 -0.00536222989 -3.23831009e-05 -5.43431997e-06
50 0.867900014 -0.0060936301 -0.000113898001 3.32484001e-06
55 0.834999979 -0.00698324991 -6.40252983e-05 9.34959019e-07
60 0.798600018 -0.00755337998 -5.0000901e-05 9.3532401e-07
65 0.7597 -0.0079832403 -3.59709993e-05 -2.27625992e-06
70 0.718599975 -0.00851367041 -7.01148965e-05 -8.63029982e-06
75 0.673200011 -0.00986209046 -0.000199568996 1.91973995e-05
80 0.621299982 -0.0104179997 8.83922985e-05 6.24051017e-06
85 0.5722 -0.00906600989 0.000182000003 6.24051017e-06
90 0.532199979 -0.00677797012 0.000275608007 6.24051017e-06

纵向系数 Py

纬度起点 (°) a b c d
0 -5.20417001e-18 0.0124000004 1.21431003e-18 -8.45283982e-11
5 0.061999999 0.0124000004 -1.26792998e-09 4.22642005e-10
10 0.123999998 0.0124000004 5.07171016e-09 -1.60603997e-09
15 0.186000004 0.0123998998 -1.90189002e-08 6.00152017e-09
20 0.247999996 0.0124001997 7.10039032e-08 -2.24000001e-08
25 0.310000002 0.0123992004 -2.64996999e-07 8.35986e-08
30 0.372000009 0.0124028996 9.88982947e-07 -3.11994e-07
35 0.433999985 0.0123893004 -3.69093004e-06 -4.3562099e-07
40 0.495799989 0.0123198004 -1.02251997e-05 -3.45523006e-07
45 0.557099998 0.0121916002 -1.54081008e-05 -5.8228801e-07
50 0.617600024 0.0119938003 -2.41424004e-05 -5.25326982e-07
55 0.676900029 0.011713 -3.20223007e-05 -5.16405009e-07
60 0.734600008 0.0113541 -3.97684016e-05 -6.09051995e-07
65 0.790300012 0.0109107001 -4.89041995e-05 -1.04739001e-06
70 0.843500018 0.0103430999 -6.46150002e-05 -1.40374001e-09
75 0.893599987 0.00969685987 -6.46359986e-05 -8.54700011e-06
80 0.939400017 0.00840947032 -0.000192841006 -4.21059985e-06
85 0.976100028 0.00616527023 -0.000255999999 -4.21059985e-06
90 1 0.00328946998 -0.000319159008 -4.21059985e-06

3.18 温克尔三重 Winkel Tripel（wintri）

它是艾托夫投影和等距圆柱投影的坐标平均。标准形式通常取 cosφ₁=2/π，本项目默认以近似值 φ₁=50.467° 实现；与参考实现比较时应传入相同值，不要把近似默认值和更长有效位的标准纬线混用。

正算与反算公式

α=arccos⁡ ⁣(cos⁡φcos⁡L2),sinc⁡α={sin⁡α/α,α≠0,1,α=0,xA=2cos⁡φsin⁡(L/2)sinc⁡α,yA=sin⁡φsinc⁡α,u=Lcos⁡φ1+xA2,v=φ+yA2.\begin{aligned} \alpha&=\arccos\!\left(\cos\varphi\cos\frac L2\right),\\ \operatorname{sinc}\alpha&=\begin{cases}\sin\alpha/\alpha,&\alpha\ne0,\\1,&\alpha=0,\end{cases}\\ x_A&=\frac{2\cos\varphi\sin(L/2)}{\operatorname{sinc}\alpha},& y_A&=\frac{\sin\varphi}{\operatorname{sinc}\alpha},\\ u&=\frac{L\cos\varphi_1+x_A}{2},&v&=\frac{\varphi+y_A}{2}. \end{aligned}(36)

(u(L,φ)v(L,φ))=((X−x0)/R(Y−y0)/R),λ=λ0+L.\begin{pmatrix}u(L,\varphi)\\v(L,\varphi)\end{pmatrix} =\begin{pmatrix}(X-x_0)/R\\(Y-y_0)/R\end{pmatrix}, \qquad \lambda=\lambda_0+L.(37)

反算没有采用闭式表达，用第 4 节二维阻尼牛顿法求 λ,φ。应限制合法经差和纬度，并回代正算检查残差。名称中的“三重”不能理解为面积、角度、距离都严格保持；这是折衷投影。PROJ：Winkel Tripel、PROJ 艾托夫与温克尔实现

4. 数值算法和异常处理

4.1 二维阻尼牛顿反算

温克尔三重使用正算函数 F(λ,φ) 反求经纬度。高斯－克吕格现使用3.2节的逆级数，不再经过本节二维迭代。给定目标 (X*,Y*)：

迭代求解必须使用实际正算函数及其雅可比，不能把“反算”留作一句概念说明：

z=(Lφ),r(z)=F(z)−(X∗Y∗),JF=(∂X/∂L∂X/∂φ∂Y/∂L∂Y/∂φ)=(abcd),det⁡JF=ad−bc,ΔL=d rX−b rYad−bc,Δφ=−c rX+a rYad−bc,zj+1=zj−sj(ΔLΔφ).\begin{aligned} \mathbf z&=\begin{pmatrix}L\\\varphi\end{pmatrix},& \mathbf r(\mathbf z)&=F(\mathbf z)-\begin{pmatrix}X^*\\Y^*\end{pmatrix},\\ J_F&=\begin{pmatrix} \partial X/\partial L&\partial X/\partial\varphi\\ \partial Y/\partial L&\partial Y/\partial\varphi \end{pmatrix}=\begin{pmatrix}a&b\\c&d\end{pmatrix},\\ \det J_F&=ad-bc,\\ \Delta L&=\frac{d\,r_X-b\,r_Y}{ad-bc},& \Delta\varphi&=\frac{-c\,r_X+a\,r_Y}{ad-bc},\\ \mathbf z_{j+1}&=\mathbf z_j-s_j\begin{pmatrix}\Delta L\\\Delta\varphi\end{pmatrix}. \end{aligned}(38)

∂F∂zk(z)≈F(z+hek)−F(z−hek)2h,h=10−6 rad.\frac{\partial F}{\partial z_k}(\mathbf z) \approx\frac{F(\mathbf z+h\mathbf e_k)-F(\mathbf z-h\mathbf e_k)}{2h}, \qquad h=10^{-6}\ \mathrm{rad}.(39)

  1. 从投影中心和近似经纬度开始，计算残差 r=F(λ,φ)−(X*,Y*)。代码直接调用相对经度形式的 raw()，再施加比例与假坐标，不在迭代内部反复包装绝对经度。
  2. 用小角度差分近似两个偏导，组成2×2雅可比矩阵 J。内部使用中心差分，靠近经差和纬度边界时使用域内截断差分，并按实际采样间距计算导数。
  3. 求解线性方程 J Δ=r。若行列式过小，说明局部求逆不稳定，应换初值或报告失败。
  4. 试取 (λ,φ)−sΔ，从 s=1 开始逐次减半，只有残差下降才接受。
  5. 保持纬度与经差在允许域内；达到迭代上限或持续不下降时返回错误。
  6. 用求出的经纬度再次正算，检查米制残差，确认不是“迭代停止了但坐标不对”。

当前实现最多60次外层迭代，雅可比名义差分步长 10⁻⁶ rad，二维正算残差小于 10⁻⁵ m 时接受迭代结果；阻尼步长从1逐次减半，尝试至大于 1/1024。候选纬度限制与统一输入纬度一致，即 |φ|<89.999999°；经差 |L|≤π。越界候选投回允许区间后，仍必须让残差下降才接受。

全部反算另经过统一正算回代。普通坐标接受门限为0.1m；接近奇点、平面坐标非常大时，增加相对浮点容差：

∥F(λ,φ)−(X,Y)∥2≤max⁡ ⁣(0.1 m,10−11(X−x0∗)2+(Y−y0)2).\|F(\lambda,\varphi)-(X,Y)\|_2\le \max\!\left(0.1\ \mathrm m,10^{-11}\sqrt{(X-x_0^*)^2+(Y-y_0)^2}\right).(40)

这不是定位精度承诺。回代测试必须配合独立正算参考及反算经纬度检查，不能只看程序有没有抛出异常。接缝经差在 10⁻¹² rad 舍入邻域统一到−π分支，避免同一边界在正反算中跳到地图另一侧；这只统一数值分支，不改变180°经线的地理位置。

这类反算在奇点、折叠或域边缘附近可能失败。失败应在结果中显式呈现，不得把上一行坐标或零值当作有效结果。

4.2 一维辅助角求解

摩尔威德求辅助角，等积地球反求 θ，罗宾森反求纬度，都是一维问题。优先利用单调性和有界区间。二分法收敛稳定，牛顿法收敛较快，但靠近零导数或定义域边缘时必须保护步长。

当前代码中，等积地球采用牛顿法，罗宾森采用二分法。摩尔威德仍解同一个辅助角方程，但接近两极时改解余角 z=π/2−|θ|，避免 πsinφ−π 的数值消减：

δ=π2−∣φ∣,2z−sin⁡2z=2πsin⁡2δ2,2z−sin⁡2z=43z3−415z5+8315z7−42835z9+8155925z11+O(z13).\begin{aligned} \delta&=\frac\pi2-|\varphi|,&2z-\sin2z&=2\pi\sin^2\frac\delta2,\\ 2z-\sin2z&=\frac43z^3-\frac4{15}z^5+\frac8{315}z^7-\frac4{2835}z^9+\frac8{155925}z^{11}+O(z^{13}). \end{aligned}(41)

代码在 z<0.1 rad 时用上式多项式计算左端，导数为 4sin²z，以三次项给出初值，最多30次迭代。它仍是原摩尔威德投影，不是替换了投影定义。

方位类投影接近中心或对跖点时也要避免只使用 acos(q) 或直接相减 1+q。令 q=sinφ₀sinφ+cosφ₀cosφcosL、E=cosφsinL、B=cosφ₀sinφ−sinφ₀cosφcosL、t=hypot(E,B)，代码使用：

c=atan2⁡(t,q),kaeqd=ct,klaea=2(1−q)t(q<0),kstere=2(1−q)t2(q<0).\begin{aligned} c&=\operatorname{atan2}(t,q),&k_{\mathrm{aeqd}}&=\frac ct,\\ k_{\mathrm{laea}}&=\frac{\sqrt{2(1-q)}}t\quad(q<0),& k_{\mathrm{stere}}&=\frac{2(1-q)}{t^2}\quad(q<0). \end{aligned}(42)

这些是常见公式的代数等价稳定表达；中心取连续极限。当前仍拒绝 q≤−1+10⁻¹² 的对跖点邻域，以及超过 10¹² m 的平面分量。全球地图的显示窗口比坐标算法支持域更窄，二者应分别解释。

4.3 误差应分开说明

类型 例子 正确评价方法
模型差异 球面墨卡托与 WGS84 椭球墨卡托不同 对照时使用相同 R 或椭球
浮点与级数截断误差 六阶 Krüger 与参考实现的差值 使用独立 PROJ 参考值，在给定测试域内统计
反算数值误差 牛顿迭代未充分收敛 查看正算回代残差与经纬度差
数据误差 简化国界、坐标偏移、源精度不足 核查源数据，不用换投影掩盖
显示误差 曲线顶点过稀、画布缩放 加密顶点并检查绘制流程

往返测试 inverse(forward(P))≈P 能发现不一致，但如果正反算共享同一错误，往返仍可能通过。因此应同时用独立参考实现验证正算值。测试报告中的最大误差只覆盖所列样本与参数，不代表所有输入、所有投影边缘都获得同样精度。

5. 从坐标转换到地图转换

地图转换不是对屏幕截图做拉伸。其数据流程应是：

经纬度几何 → 坐标有效性检查 → 逐点投影
→ 处理不可见点、奇点和不连续处 → 生成平面几何
→ 平移缩放到画布 → 绘制边界、经纬网与辅助图形

点、线、多边形都可以逐顶点投影；但要让一条原本直的经纬度线段在投影中正确弯曲，需要在源空间增加采样点。遇到投影裁切边界时，多边形还需要重建闭合环。教学实现的逐顶点转换不能自动等同于完整的 GIS 拓扑裁切引擎。跨 180° 经线、接近地平线、穿越球心投影无穷远边界的几何，应单独检查。

当前路径为 Geometry::project() 递归处理几何层级，position() 调用 Projection::forward()；地图预览将不可投影顶点记为 null 并统计数量，使绘图端可以断开线段。严格转换模式遇到不可投影顶点即报错。输入几何默认顶点上限为 60000；第三维数值若存在，仅原样保留，不实施高程变换。上传几何不会自动加密边界，也不会自动修复拓扑。

5.1 GeoJSON 的坐标语义

标准 RFC 7946 第 4 节规定 GeoJSON 使用 WGS84 经纬度。把 [经度,纬度] 改为米制 [X,Y] 后，即便几何结构仍像 GeoJSON，也已不再是符合 RFC 7946 的经纬度 GeoJSON。

因此，项目导出的米制 JSON 必须保留投影标识、全部参数、单位和模型，并明确标注为投影坐标 JSON／GeoJSON 结构的教学交换文件。接收方必须主动设置正确 CRS。若要得到标准 GeoJSON，应该把结果反算回 WGS84 经纬度后再导出。输入文件中的旧式 crs 字段也不能代替真实的 CRS 识别与变换。

5.2 数据来源与完整性

内置中国轮廓用于课堂观察投影几何，不应把它当作经审定的标准行政区划成果。实际数据来源、许可、简化级别与缺少的细节以项目数据说明为准。行政边界的精度、岛屿完整性与投影算法是不同问题；任一算法通过数值测试，不会提高源边界数据的权威性或精度。

5.3 中国与全球底图切换

网站保留中国示例，并增加 region=world 的全球陆地与海岸线。全球源数据是 Natural Earth 1:110m land（公共领域，含南极洲）；不使用它表达国家政治边界。中国与全球分别保留自己的投影参数；全球初始中央经线与纬度原点为0°，其余参数按同一投影目录确定。

坐标正反算仍由同一个 Projection 类执行。新增 WorldMap 只处理全球地图的显示范围、接缝与绘制几何，不改变任何投影定义或正反算算法。

经度接缝处理：预先把陆地分成不超过5°的经纬度小片，并按最大0.75°步长加密。对非方位类方法，在相对中央经线的经度空间平移、裁切，然后逐点调用PHP正算，避免把左右图边直接连接。

L~=λ−λ0+2kπ,−π+ε≤L~≤π−ε,k∈{−1,0,1}.\widetilde L=\lambda-\lambda_0+2k\pi,\qquad -\pi+\varepsilon\le\widetilde L\le\pi-\varepsilon, \qquad k\in\{-1,0,1\}.(43)

方位类显示范围：以球心角限制可显示区域，利用边界求交裁切小片，再调用相同PHP投影。

q=sin⁡φ0sin⁡φ+cos⁡φ0cos⁡φcos⁡(λ−λ0),q≥cos⁡cmax⁡.q=\sin\varphi_0\sin\varphi+ \cos\varphi_0\cos\varphi\cos(\lambda-\lambda_0),\qquad q\ge\cos c_{\max}.(44)

方法 全球教学显示范围 原因
墨卡托 纬度±85° 两极纵坐标发散
高斯－克吕格 经差按模式取±45°、±1.5°或±3°；纬度±85° 宽带仅供变形演示；3°/6°模式显示一个分带
兰伯特正形圆锥 纬度±80° 避免对侧极区发散压缩其余画面
正射方位 距中心89.9°以内 接近可见半球，背面不可见
球心 距中心80°以内 90°边界发散
球面立体 距中心150°以内 避免对跖点附近无限放大
方位等距、兰伯特方位等积 距中心178°以内 对跖点邻域方向奇点
其他方法 经度全范围、纬度±89.5° 避开精确极点与经度不唯一问题

这是全球数据在各投影下的有限、明确注明的课堂显示窗口，不表示这些投影都能无变形、无裁切地表示全球。全球墨卡托、正射、球心的差异本身也是重要的教学内容。

高斯单带图为何只显示部分地球：单带扩展模式限制为中央经线左右各45°，合计90°经度窗口；3°／6°模式只绘制当前一个分带。这不是“高斯只能投影半个地球”的数学结论。当前全球模式新增默认“全球分带展开”，具体算法见5.4节；也可切回单带进行坐标实验。横轴墨卡托的其他实现可扩展更大的定义域，但不能据此把本程序的六阶算法任意放宽而不复核边界精度。

小片裁切在经纬度空间中完成；弧形显示边界用细分短线近似。海岸线独立投影，遇到接缝或不可显示点就断开；陆地填充不绘制小片内部网格。此流程用于教学显示，不是测绘级拓扑裁切，也不应把裁切小片当作新的行政或自然边界。

全球模式的JSON导出是当前显示范围内裁切后投影的陆地与海岸线，包含投影、参数、模型和显示范围说明；不能称为未裁切的完整全球成果。中国与用户上传数据的严格导出规则保持原样。

全球经纬网覆盖各经度与南北纬；辅助圆统一使用300km球面半径，中国示例继续使用100km。比较不同投影时保持同一实验区域、同一圆半径和同一尺度设置。区域切换不改变单独坐标转换的数学定义域。

来源：Natural Earth Land、Natural Earth Terms of Use。预处理脚本 tools/prepare-world.py 使用Shapely生成静态数据；运行网站不依赖Python或Shapely。

5.4 高斯全球分带展开算法

在全球模式下，高斯默认展示6°、60带展开，也可切换3°、120带；单图、对比和18图总览均使用同一数据包。其余投影与中国高斯仍保持原有行为。展开覆盖全部经度，纬度裁切为±89.5°以避开精确极点经度不唯一；这不是未经裁切的完整全球测绘成果。

分带与独立坐标。按既有约定，6°带号为1～60，3°为1～120；经度归一化到[-180°,180°)，按中央经线从西到东排列。3°第60带以-180°为中心，跨日期变更线，其两侧数据归入同一带。带边界归属采用左闭右开，避免同一样点重复归属。

N=360/w,w∈{3∘,6∘},λ0,n={6n−3,w=6∘,3n,w=3∘,Ln=wrap⁡[−180,180)(λ−λ0,n),−w/2≤Ln<w/2.\begin{aligned} N&=360/w,\qquad w\in\{3^\circ,6^\circ\},\\ \lambda_{0,n}&=\begin{cases}6n-3,&w=6^\circ,\\3n,&w=3^\circ,\end{cases}\\ L_n&=\operatorname{wrap}_{[-180,180)}(\lambda-\lambda_{0,n}),\qquad -w/2\le L_n<w/2. \end{aligned}(45)

每个带调用现有WGS84椭球六阶Krüger正反算；固定纬度原点0°、中央比例因子1、假东距500000m、假北距0、不加百万位带号前缀。带号与中央经线是坐标系元数据；不为这一组自定义参数虚构一个统一EPSG编号。侧栏参数用于当前单带坐标实验，展开自身采用上述固定参数。点击侧栏“应用参数”或由助手应用高斯参数时，切回单带视图，避免把自定义参数误认为已用于全部展开带。

布局只作平移。令h为某带赤道东侧边界的东坐标减500000m，统一取P=2h×1.08为带间中心距；将带按中央经线排序，序号i=0,…,N−1。保留每带相同的米制尺度，不压缩高纬带宽，也不拉伸纬向。

h=Xn(λ0,n+w/2,0)−500000,P=2h(1+0.08),Oi=(i−N−12)P,Xlayout=Xn−500000+Oi,Ylayout=Yn.\begin{aligned} h&=X_n(\lambda_{0,n}+w/2,0)-500000,\\ P&=2h(1+0.08),\qquad O_i=\left(i-\frac{N-1}{2}\right)P,\\ X_{\mathrm{layout}}&=X_n-500000+O_i,\qquad Y_{\mathrm{layout}}=Y_n. \end{aligned}(46)

全球底图在每带经纬度范围内先裁切、再投影；经纬网、城市及辅助图层以同样平移进入画布。带框保留高纬收窄的形状。图顶显示带号与中央经线；总览时自动减少标签，放大后逐步显示更多标签，完整目录仍可从“查看分带”选取。点击“进入该带坐标实验”可切到该带单图，参数与坐标算法一致。

取点与导出。先把屏幕点击还原到布局米坐标，用布局槽位确定候选带，再扣除平移并调用该带反算；若位于间隔、带框外或显示纬度外则拒绝。取点成功时保留缩放和平移，将侧栏带号、中央经线、经纬度同步到工作台。

Xn=Xlayout−Oi+500000,Yn=Ylayout,(λ,φ)=Fn−1(Xn,Yn).X_n=X_{\mathrm{layout}}-O_i+500000,\qquad Y_n=Y_{\mathrm{layout}},\qquad (\lambda,\varphi)=F_n^{-1}(X_n,Y_n).(47)

SVG／PNG保留展开排版；投影JSON导出为bands[]，每个带分别包含原始本带几何、完整参数、轴序X东Y北、单位米和仅供绘图的display_offset_x。不会把布局坐标冒充一个统一平面CRS，也不把不同带东坐标直接拼接为测绘成果。比例尺只用于带内的投影长度，带间空隙不代表地面距离，不能跨带量距。辅助圆在带边界裁切，中心路径分别以各带原点生成。

实现位置：src/GaussAtlas.php、API操作atlas／atlas_export和前端高斯全球视图。独立检查见tests/verify_gauss_atlas.py：覆盖全部经度唯一归属、60／120个带号、日期变更线、北南高纬、逐带与pyproj正算对照，以及导出本带坐标与绘图平移的分离。

6. 变形、距离与显示尺度

6.1 局部变形

在球面上，以一点的东、北单位地面方向为基，投影局部雅可比可写成：

Jground=(∂X/∂λRcos⁡φ∂X/∂φR∂Y/∂λRcos⁡φ∂Y/∂φR).J_{\mathrm{ground}}= \begin{pmatrix} \displaystyle\frac{\partial X/\partial\lambda}{R\cos\varphi}& \displaystyle\frac{\partial X/\partial\varphi}{R}\\[6pt] \displaystyle\frac{\partial Y/\partial\lambda}{R\cos\varphi}& \displaystyle\frac{\partial Y/\partial\varphi}{R} \end{pmatrix}.(48)

J 的两个奇异值是最大、最小局部线性尺度。正形投影的两者相等；等积投影的行列式绝对值为 1（统一比例因子为1时）。极点的 cosφ=0 需要特殊处理。对于高斯－克吕格，应把球面分母换成椭球的卯酉圈半径 N cosφ 和子午圈曲率半径 M，否则会把模型差异误当作投影角度变形。

代码 Projection::factors() 在经纬度方向以 10⁻⁵° 为名义差分步长，直接对相对经度 raw() 求导，不跨越绝对经度±180°或地图接缝。内部使用中心差分；边界处用域内单侧差分，接近不可见域时减小步长。高斯－克吕格采用上述椭球分母，其余采用球面分母。输出包括面积尺度、最大／最小尺度、东西／南北尺度与最大角度变形；其中最大角度变形为 2asin[(s₁−s₂)/(s₁+s₂)]，界面以度显示。奇点邻域的差分指标仍可能不稳定，应与所列显示域和模型一起解释。

可进一步写出面积倍率、主尺度、角度变形及椭球曲率半径：

μA=∣det⁡Jground∣,T=tr⁡(JgroundTJground),s1=T+T2−4μA22,s2=T−T2−4μA22,ωmax⁡=2arcsin⁡ ⁣(s1−s2s1+s2),N=a1−e2sin⁡2φ,M=a(1−e2)(1−e2sin⁡2φ)3/2.\begin{aligned} \mu_A&=|\det J_{\mathrm{ground}}|,& T&=\operatorname{tr}(J_{\mathrm{ground}}^\mathsf T J_{\mathrm{ground}}),\\ s_1&=\sqrt{\frac{T+\sqrt{T^2-4\mu_A^2}}{2}},& s_2&=\sqrt{\frac{T-\sqrt{T^2-4\mu_A^2}}{2}},\\ \omega_{\max}&=2\arcsin\!\left(\frac{s_1-s_2}{s_1+s_2}\right),\\ N&=\frac a{\sqrt{1-e^2\sin^2\varphi}},& M&=\frac{a(1-e^2)}{(1-e^2\sin^2\varphi)^{3/2}}. \end{aligned}(49)

其中 s₁≥s₂；椭球雅可比把球面式的 Rcosφ、R 分别替换为 Ncosφ、M。

画出有限半径地理小圆，再投影到平面，是易理解的教学近似。只有半径趋近于零时才严格接近蒂索指示椭圆；较大的地理圆可能投成非椭圆曲线。

当前 Geometry::mapLayers() 生成共同球面上的小圆：中国模式半径100km、全球模式半径300km，每圆48段，布置在多个固定经纬度位置。该图层应称为“变形辅助圆／图形”；高斯－克吕格画面中的辅助圆仍从共同球面构造，不能把它当作精确的椭球蒂索椭圆。严格局部性质以 factors() 的相应模型结果解释。

6.2 距离

球面参考距离可用大圆公式：

h=sin⁡2φ2−φ12+cos⁡φ1cos⁡φ2sin⁡2λ2−λ12,dsphere=2Ratan2⁡(h,1−h),dmap=(X2−X1)2+(Y2−Y1)2,εd=dmapdsphere−1.\begin{aligned} h&=\sin^2\frac{\varphi_2-\varphi_1}{2} +\cos\varphi_1\cos\varphi_2\sin^2\frac{\lambda_2-\lambda_1}{2},\\ d_{\mathrm{sphere}}&=2R\operatorname{atan2}(\sqrt h,\sqrt{1-h}),\\ d_{\mathrm{map}}&=\sqrt{(X_2-X_1)^2+(Y_2-Y_1)^2},\\ \varepsilon_d&=\frac{d_{\mathrm{map}}}{d_{\mathrm{sphere}}}-1. \end{aligned}(50)

它和 WGS84 椭球测地距离不是同一个量。展示高斯－克吕格的测量级误差时应使用椭球测地距离作为参考。球心图中的大圆直线性质、方位等距图中的中心保距性质，都应以匹配的球面模型解释。

6.3 自适应与统一比例尺

各图自适应画框会独立缩放，便于看清轮廓，却会隐藏总体尺度差。统一比例尺应使所有图使用相同的“像素／米”值，并保持 X、Y 相同缩放，不能把每幅图分别拉满宽高后宣称是统一比例尺。假东距、假北距只移动坐标，不改变形状或变形性质。

v4绘图先把投影坐标转换到局部画布坐标，避免带号前缀或较大假坐标超过浏览器SVG长度支持范围。几何、城市、示意圆、标签和样点使用同一转换；SVG/PNG导出使用画布坐标，数值JSON仍输出原投影米制坐标。设平面图心为 (Xc,Yc)，每画布单位代表 m 米，则：

xcanvas=X−Xcm,ycanvas=−Y−Ycm,X=Xc+mxcanvas,Y=Yc−mycanvas.x_{\mathrm{canvas}}=\frac{X-X_c}{m},\qquad y_{\mathrm{canvas}}=-\frac{Y-Y_c}{m},\qquad X=X_c+mx_{\mathrm{canvas}},\quad Y=Y_c-my_{\mathrm{canvas}}.(51)

点击地图时先执行右侧恢复关系，再调用对应投影的反算。统一尺度模式保持各图的 m 相同，平移和画布归一化不改变投影性质。

6.3.1 动态比例尺与缩放更新算法

地图左下角的“比例尺·投影平面”表示当前投影坐标中的长度。它支持中国／全球、单图／对比／总览，并在放大、缩小、平移、重置或窗口尺寸变化后更新。普通地图还显示“1 px≈…”，指一个浏览器CSS像素对应多少投影米；小图省略这行以节省空间。没有使用固定的1000km线段，也没有把屏幕像素直接等同于印刷厘米。

1. 当前分辨率。 渲染时保存 M₀=data-mpp，表示一单位局部SVG画布坐标对应的投影米数。SVG的 getScreenCTM() 把当前画布坐标转换到屏幕CSS像素，其线性部分为矩阵A。当前地图只做统一缩放与平移，不做非等比例拉伸，因此：

A=(acbd),s=a2+b2[px/画布单位],μ=M0s[m/px].A=\begin{pmatrix}a&c\\b&d\end{pmatrix},\qquad s=\sqrt{a^2+b^2}\quad[\mathrm{px}/\text{画布单位}],\qquad \mu=\frac{M_0}{s}\quad[\mathrm m/\mathrm{px}].(52)

这里的屏幕像素是CSS像素，不是显示器硬件像素。直接读取SVG实际变换矩阵，能够把viewBox缩放、居中留白和画框实际尺寸纳入计算，而不是假设画框宽度永远等于初始viewBox宽度。

2. 选择整洁刻度。 普通图目标尺长至多120px，总览等小图至多90px，同时根据画框宽度自适应。设画框屏幕宽度为W，目标尺长为B；代码以36px作为目标的下限：

B=min⁡ ⁣(Bmax⁡,max⁡(36,0.26W)),Bmax⁡={120,普通图,90,总览或小图,T=μB,n=⌊log⁡10T⌋,q=T10n,c={5,q≥5,2,2≤q<5,1,1≤q<2,D=c 10n,Lpx=Dμ.\begin{aligned} B&=\min\!\left(B_{\max},\max(36,0.26W)\right), &B_{\max}&=\begin{cases}120,&\text{普通图},\\90,&\text{总览或小图},\end{cases}\\ T&=\mu B,&n&=\lfloor\log_{10}T\rfloor,\quad q=\frac{T}{10^n},\\ c&=\begin{cases}5,&q\ge5,\\2,&2\le q<5,\\1,&1\le q<2,\end{cases} &D&=c\,10^n,\qquad L_{\mathrm{px}}=\frac D\mu. \end{aligned}(53)

因此完整刻度D取 1、2、5×10ⁿ 米系列，实际尺长不会超过目标B。黑白尺分成四等份，在起点标0、中央标D/2，终点旁标完整长度与单位。地图宽度小于280px、高度小于220px或处于总览时使用紧凑布局。

3. 单位换算。 刻度值先保持米制计算，最终只对显示数值和单位作换算；分辨率说明独立使用同一规则：

(v,u)={(D/1000,km),D≥1000 m,(D,m),1≤D<1000 m,(100D,cm),0.01≤D<1 m,(1000D,mm),0<D<0.01 m.(v,u)=\begin{cases} (D/1000,\mathrm{km}),&D\ge1000\ \mathrm m,\\ (D,\mathrm m),&1\le D<1000\ \mathrm m,\\ (100D,\mathrm{cm}),&0.01\le D<1\ \mathrm m,\\ (1000D,\mathrm{mm}),&0<D<0.01\ \mathrm m. \end{cases}(54)

显示数字通常保留至多3位有效数字，逗号用于千位分组；过大或过小数值使用科学计数法。单位变化不改变实际的D或尺长。无效、非有限或非正的分辨率不生成标尺。

4. 随缩放更改。 代码以系数f同时乘viewBox的宽和高。点击放大时 f=1/1.3，点击缩小时 f=1.3；滚轮分别使用 1/1.12 或 1.12。在画框尺寸不变时：

Vw′=fVw,Vh′=fVh,s′=sf,μ′=fμ,D′=NiceFloor⁡(μ′B),Lpx′=D′μ′.V_w'=fV_w,\quad V_h'=fV_h,\quad s'=\frac sf,\quad \mu'=f\mu,\qquad D'=\operatorname{NiceFloor}(\mu'B),\quad L_{\mathrm{px}}'=\frac{D'}{\mu'}.(55)

放大后每像素代表的米数减少，比例尺显示长度相应减小；由于采用1、2、5系列，标注值会分级跳变，而不必每次缩放都改变数字。例如 μ=1000m/px、B=120px 时选择D=100km、尺长100px；放大到2倍后 μ=500m/px，选择D=50km、尺长仍为100px。如果刻度暂未跳变，尺条本身的像素长度会变化，仍满足 D=μLpx。

5. 固定在画面左下角。 取画框左边和底边各14px的屏幕锚点P，用当前SVG完整变换矩阵C的逆矩阵恢复画布坐标。标尺组再施加反向缩放，确保文字和尺条厚度不会随地图放大：

p=C−1 ⁣(xleft+14ybottom−141),G=Translate⁡(px,py)Scale⁡(1/s).\mathbf p=C^{-1}\!\begin{pmatrix}x_{\mathrm{left}}+14\\y_{\mathrm{bottom}}-14\\1\end{pmatrix}, \qquad G=\operatorname{Translate}(p_x,p_y)\operatorname{Scale}(1/s).(56)

标尺内部按CSS像素布局。平移时C的平移分量改变，但s和μ不变；重新计算p即可让标尺留在画框左下角，刻度值不变。重置时恢复初始viewBox，重新计算C、μ、D和位置。响应式尺寸变化与区域、投影、统一尺度切换后，地图重绘并重复上述流程；对比图的每幅地图各自生成标尺。

6. 导出与地面距离解释。 SVG／PNG保留比例尺图形和长度标注，它们与地图一同缩放，D与图上几何的对应关系不变。导出时移除“1 px≈…”这行屏幕分辨率说明，因为文件可能被再次缩放、以1.5倍像素栅格化或打印；它不能继续当作导出文件每一像素的米数。标尺通过SVG图形绘制，不依赖仅在网页可见的CSS叠层。

当前标尺严格表示投影平面长度。投影局部线性尺度k随位置、方向和投影改变；在足够小的范围内，相关方向上可近似 D地面≈D平面/k，不能把同一平面标尺解释为整幅全球地图处处相同的真实地面距离。网站没有额外进行这种局部地面修正，也不根据未知显示器物理尺寸标称精确的1∶N印刷比例尺。

实现对应 public/assets/scale-bar.js 的刻度、单位选择与 app.js 的 updateScaleBar()。108项数学检查覆盖固定算例、缩放关系、刻度系列、长度单位、不超过目标尺长及无效输入；浏览器检查另验证实际矩阵、位置、缩放重置和导出。

6.3.2 视图保留与固定屏幕标记

取点仅重新正算并更新样点图层，不重建地图SVG，因此保留当前viewBox、缩放、平移和比例尺。连续取点使用请求序号避免较早的异步响应覆盖新标记。图层重绘及工作台返回地图时，在投影、参数、区域、数据和尺度模式未变的情况下保留缩放倍率与投影坐标中心；切换投影、区域或应用参数则重新适配。

城市点、取点样点及上传GeoJSON点是屏幕符号。设目标半径为CSS像素r_px，由6.3.1节屏幕变换取得s，则SVG半径取r_svg=r_px/s。城市点半径2.5px，取点半径5px，上传点半径3px；按钮和滚轮缩放后重新计算。变形辅助圆表示地理半径，仍随地图放大，以保留其教学含义。

6.4 完整算例：郑州点的墨卡托正算与反算

使用与网站一致的球面与默认墨卡托参数：R=6371007.180918 m、λ₀=105°、φ₁=0°、x₀=y₀=0。输入郑州示例点 λ=113.6254°、φ=34.7466°。

L=(113.6254−105)π180,φ=34.7466π180,k=cos⁡0=1,X=RL≈959101.801307 m,Y=Rln⁡tan⁡(π4+φ2)≈4124882.023450 m,λdeg=(XR+105π180)180π≈113.6254000000∘,φdeg=[2arctan⁡ ⁣(eY/R)−π2]180π≈34.7466000000∘.\begin{aligned} L&=(113.6254-105)\frac\pi{180},\\ \varphi&=34.7466\frac\pi{180},\qquad k=\cos0=1,\\ X&=R L\approx959101.801307\ \mathrm m,\\ Y&=R\ln\tan\left(\frac\pi4+\frac\varphi2\right) \approx4124882.023450\ \mathrm m,\\ \lambda_{\mathrm{deg}}&=\left(\frac XR+105\frac\pi{180}\right)\frac{180}{\pi} \approx113.6254000000^\circ,\\ \varphi_{\mathrm{deg}}&=\left[2\arctan\!\left(e^{Y/R}\right)-\frac\pi2\right]\frac{180}{\pi} \approx34.7466000000^\circ. \end{aligned}(57)

这些数值用于检查单位、轴序和正反算流程；小数位数不表示原始位置或边界数据具有同等现实精度。

7. AI 与算法的职责

AI 将用户语言转换成受限的操作对象：选择已有投影、填入已知参数与坐标、切换显示选项、解释推荐理由。参数必须经后端白名单、类型和范围校验，再进入投影计算。坐标结果由 PHP 算法返回，不由语言模型猜测。

当前操作字段包括 projection、params、mode、compare、layers、tab，以及 coordinates、direction、same_scale、region。region只允许china或world，助手能够切换区域并基于当前区域推荐投影。coordinates 最多填写 1000 个有限数值对；正向时为 [经度,纬度]，反向时为 [X,Y]，方向只允许 forward 或 inverse，统一尺度只接受布尔值。AI 将这些内容填入坐标表单后，由用户点击“开始转换”才运行 PHP 计算；填写成功本身不是转换已完成。

用户说“比较省份面积”，可以优先推荐阿尔伯斯；说“从郑州看各城市距离”，可以推荐以郑州为中心的方位等距；说“得到真实高斯坐标”，则需要明确中央经线、分带、假东距、基准和轴序。AI 缺少这些信息时不应冒充已知。未配置外部模型时，本地规则只能覆盖预设语言模式，不能称为已连接的通用大模型。

7.1 kimi-k3真实模型接入

2026-10-19本机已接入https://sub2.guuiy.cn/v1/chat/completions、模型kimi-k3。接口实测要求temperature=1，使用非流式Chat Completions与JSON对象输出，超时90秒。此前本地规则验收记录属于当时版本，当前接入结果以本节及验收记录最新节为准。

请求由PHP服务端发送，仅发送用户本次需求、最近8条对话、当前实验状态和18种投影目录；不发送上传地图文件、CSV文件或密钥给浏览器。密钥从环境变量PROJECTION_AI_KEY或app/config/ai.local.php读取，私有文件位于public之外，不进入Git及交付包；无密钥自动使用原有本地规则。

模型输出为reply和action。服务端校验投影ID、适用参数、模式、图层、区域、坐标、方向、统一尺度，以及高斯gauss_atlas和atlas_width。通过后前端更新实验选项并提供撤销；模型仅填写坐标，正反算仍由PHP完成。自然语言的未知中心缺少经纬度时应追问，不直接猜测。模型服务错误不会伪装成本地规则推荐成功。

真实测试包含全球森林面积→等积地球、郑州中心距离→方位等距、墨卡托与阿尔伯斯对比、高斯3°第38带坐标填表、全球120带展开及未知中心追问。最后一项曾出现一次服务错误，重试后通过；不能据有限测试保证每次自然语言推荐都完全一致。原始响应和浏览器实际自动填参结果分别见tests/ai-first-live-result.json、tests/ai-live-results.json、tests/ai-browser-results.json，不包含密钥。

8. 验证与复核方法

建议至少覆盖以下有效域内样本：投影中心、经纬网样本、中国东西南北代表点、南半球合法参数、正反算往返、定义域外坐标和退化参数。

独立 PROJ 比对时，17 种球面方法须显式使用 +R=6371007.180918，高斯－克吕格须用 +ellps=WGS84；中央经线、原点、标准纬线、比例因子和假坐标必须一致。不要用 PROJ 默认的 GRS80 椭球结果检验本项目球面实现。

球面示例：+proj=aea +R=6371007.180918 +lon_0=105 +lat_0=0 +lat_1=25 +lat_2=47
椭球示例：+proj=tmerc +ellps=WGS84 +lon_0=105 +lat_0=0 +k_0=1 +x_0=500000

8.1 已执行的验证结果

2026-10-19 已完成以下本地验证，原始结果保存在项目测试目录：

验证 证据文件／程序 实际结果
独立正算与往返参考比对 tests/reference-results.json、tests/verify_reference.py 18 种投影，3 组参数，共 2159 个有效样点，失败 0
PHP HTTP API 功能检查 tests/api-results.json、tests/verify_api.py 80 项通过，失败 0
AI 操作校验与负中央经线回归 tests/verify_validation.php 10 项通过，失败 0

独立参考版本为 pyproj 3.7.2 / PROJ 9.5.1。样本包含每组参数的 5 个固定点及 35 个固定随机种子的区域样点；随机点覆盖经度 [75°,134°]、纬度 [5°,54°]，另含纬度 0°、4°等固定点。三组中央经线为 105°、112°、90°，并按方法匹配原点、标准纬线、比例因子、假东距和地球模型；高斯－克吕格排除超出所设测试经差的一个点。完整参数和选点规则见测试脚本。

v4重新运行这2159个基础样本，最大正算平面差为 4.8617×10⁻⁹ m（米勒圆柱），最大正反算经纬度差为 3.1293×10⁻¹¹°（温克尔三重，按经度差与纬度差绝对值的较大者统计）。测试通过门限为正算差 0.05 m、往返经纬度差 10⁻⁷°。这些只描述基础区域样本；全球极区结果另见8.2节。

这些极小数值描述的是本项目公式实现与同模型参考程序在这批样本上的数值一致性，不代表真实地图、源数据或实际定位达到纳米精度，也不证明全部全球区域、极点邻域、定义域边缘、全部参数组合或所有上传几何均已验证。17 种球面方法仍存在相对于真实椭球的模型近似；内置边界精度也不由这些数值决定。

正算行结果在坐标有效但局部差分或辅助反算受边界影响时，仍保留正算输出，并以 factors_notice / roundtrip_notice 分别说明，避免把辅助指标失败当作正算失败。

API 检查覆盖 18 种默认投影的中国地图、坐标正反算，部分行错误、退化圆锥参数、未知投影、上传几何及严格导出、不可投影顶点的预览断开、本地助手常见意图和私有源码访问隔离。PHP 校验检查包含拒绝非法操作，以及负中央经线的高斯－克吕格往返回归。它们不等同于已验证外部大模型服务；真实 AI 接口仍须在用户填写配置后联调。

8.2 全面核查修复版 v4

2026-10-19进一步扩展到全球、南半球、非零原点、偏移、比例因子、地图接缝与高纬度，并保留修复前结果。测试门限和选点均写在脚本中：

验证 样本／检查 修复后结果 证据
18种投影扩展参考与往返 9906个有效域样点 失败0 tests/extremes-results.json
高斯专项 1676个样点，含±45°和极区 失败0 tests/tmerc-audit-results.json
对跖点邻域专项 72个案例 支持域内通过；按预定保护条件拒绝12例 tests/antipode-results.json
修复、分带、导出、AI填表、变形性质 332项检查 通过332项 tests/repairs-results.json

扩展参考使用相同球面或椭球模型，与PROJ统一经差接缝分支。正算接受差为0.05m或参考坐标长度的 10⁻¹²，取较大者；往返经纬度通常要求 10⁻⁶°，极接近极点（|纬度|>89.99°）放宽至 2×10⁻⁵°，因为极区经度的条件数很大。摩尔威德在89.9999°附近与本次PROJ版本约有厘米级辅助角数值差，不能把区域样点的纳米级一致性外推到极区。

对跖点专项使用 pyproj.Geod(a=R,f=0) 独立计算距离与方位，再按方位投影的径向关系构造参考值。原因是同一PROJ版本的部分旧投影公式在极近对跖点时也会出现数值消减，单独与其投影函数一致不能证明几何结果正确。实际检查包括反算位置，而不只有回代残差。12例拒绝包括9例体视坐标超过±10¹²m、3例位于明确的对跖点保护邻域，不能算成成功转换。

修复前的扩展报告 tests/extremes-before.json 有215个失败样本；这些是重复触发的样本数，不是215种不同错误。改动原因、复现点与范围见《18种投影核查与修复报告_20261019.md》。这些检查仍是有限样本验证，不是对全部参数、全部奇点或测绘生产精度的证明。

9. 主要资料

各方法小节已附 PROJ 官方文档链接；数学表达按本项目算法路径组织。基础理论和常见球面公式可进一步查阅 Snyder 的 USGS 手册。更新的系数和数值方法参照 PROJ 与 Karney 的原始资料。

本技术文档区分“投影定义成立”“代码实现该算法”“实际样本通过验证”三种证据层次。算法说明用于理解和复核；生产测绘任务还应采用经验证的完整 CRS 引擎与合格源数据。

公式已预排版，字体全部内嵌；无需联网或JavaScript即可阅读。可使用浏览器打印功能输出。数学排版：KaTeX 0.19.0。
数学排版组件许可
The MIT License (MIT)

Copyright (c) 2013-2020 Khan Academy and other contributors

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.