地球表面是一個(gè)不規(guī)則的曲面，而海圖卻是平面的，從不規(guī)則曲面到平面的過程被稱為投影，這種投影必然在不同的地點(diǎn)會(huì)有不同的變形。傳統(tǒng)航海確定方向的主要工具是指南針，如果一個(gè)人照著指南針給定的固定角度航行，肯定被認(rèn)為是走的直線，這種”直線“的性質(zhì)也應(yīng)該體現(xiàn)在海圖上；同時(shí)，現(xiàn)實(shí)世界中前后左右的物體被投影到海圖其相對(duì)方位也應(yīng)該保持不變，即“等角”性質(zhì)。為了滿足航海人員對(duì)海圖的“直線”和“等角”的性質(zhì)，荷蘭地圖學(xué)家墨卡托在1569年擬定了一種投影規(guī)則 -- “等角正圓柱投影”。因此，利用等角正圓柱投影的海圖被稱作墨卡托海圖，世界上95%以上的海圖都是墨卡托海圖。
由上圖可知墨卡托海圖上：等經(jīng)差在海圖上等間距；等緯差在海圖上不等間距，且隨緯度增加，間距逐漸加大。
恒向線
比如，設(shè)定角度60度，一直沿著這個(gè)角度航行，真實(shí)航行軌跡是下圖的曲線，而非直線。該曲線被稱之為恒向線，是一條與所有經(jīng)線相交成恒定角度的、具有雙重曲率的球面螺旋線，它逐漸趨向地極，但永遠(yuǎn)達(dá)不到地極。
以地球?yàn)閳A球體簡(jiǎn)化計(jì)算可知恒向線的方程如下：
對(duì)恒向線上的微分三角形DM1M2進(jìn)行積分，即可得恒向線公式：
式中：起點(diǎn)(φ1, λ1)，終點(diǎn)(φ2, λ2)，航向角C
恒向線應(yīng)用場(chǎng)景：
- 已知起點(diǎn)、終點(diǎn)，求解兩點(diǎn)間恒向線航向，及線上所有點(diǎn)的坐標(biāo)；
- 已知起點(diǎn)、航向，求解線上其他點(diǎn)的坐標(biāo)
赤道里
赤道上經(jīng)度1'所對(duì)應(yīng)的弧長(zhǎng)，各國(guó)采用的地球橢圓體參數(shù)不同而略有不同，約為1855.36米。簡(jiǎn)化計(jì)算中，赤道里長(zhǎng)度 = 海里長(zhǎng)度 = 1852米。地球赤道周長(zhǎng) = 360 * 60 * 赤道里
等角性質(zhì) -- 緯度漸長(zhǎng)率
等角意味著每一投影點(diǎn)在任一方向上的變形是相同的。根據(jù)理論，若一點(diǎn)在經(jīng)線方向與緯線方向變形一致，那該點(diǎn)在任何方向的變形都一致，即滿足等角。因此，需要了解一個(gè)新的概念：緯度漸長(zhǎng)率。緯度漸長(zhǎng)率是指在墨卡托海圖上任一緯度線到赤道的子午線的長(zhǎng)度與圖上1赤道里的比值。
對(duì)于地球的微分區(qū)域ABCD，B點(diǎn)在經(jīng)線方向與緯線方向投影變形一致，即可推出緯度漸長(zhǎng)率公式。
式中：a為地球橢圓體長(zhǎng)半軸長(zhǎng)度。緯度漸長(zhǎng)率大小只與緯度有關(guān)，與墨卡托海圖比例尺無關(guān)。可以證明，恒向線在墨卡托海圖上是直線。
海圖繪制中，需先確定1赤道里的圖長(zhǎng)，即海圖繪制單位。然后計(jì)算任意兩緯度間的緯度漸長(zhǎng)率之差，乘以海圖繪制單位，即可得兩緯線間的圖長(zhǎng)。繪制出各經(jīng)緯線圖網(wǎng)后，將地理物標(biāo)填充進(jìn)圖網(wǎng)即可。
電子海圖投影步驟
- 一英尺=25.4毫米，一般電腦屏幕每英尺像素?cái)?shù)Dpi=96，每毫米像素?cái)?shù)=96/25.4；
- 赤道里的長(zhǎng)度為1852米，每赤道里像素?cái)?shù)
MinPP=1852*100*10*96/25.4； - 以經(jīng)度0度、緯度0度為原點(diǎn)，東北方向?yàn)檎?，可將地球按照墨卡托投影成平面?/li>
- 地理坐標(biāo)(λ, φ)與平面坐標(biāo)(x, y)間可由下列公式互轉(zhuǎn)：
x = λ * 60 * MinPPy = MP(φ) * MinPPφ = AMP(y / MinPP)
MP為緯度漸長(zhǎng)率函數(shù)；AMP為緯度漸長(zhǎng)率的反函數(shù)。
- 平面坐標(biāo)(x, y)右上為正，而屏幕坐標(biāo)(x', y')原點(diǎn)在左上角，左下為正，其互轉(zhuǎn)公式如下：
x = x'y = -y'
- 電子海圖可以隨意先縮放scale（比例尺形如：1:10000，為方便計(jì)算與表示，使scale=10000），再平移(dx, dy)，因此還需要進(jìn)一步將原始屏幕坐標(biāo)(x', y')轉(zhuǎn)換為顯示坐標(biāo)(x“, y”):
x“ = x' / scale + dxy” = y' / scale - dy
- 電子海圖系統(tǒng)指定屏幕范圍內(nèi)的中心位置，或用戶跳轉(zhuǎn)到指定位置時(shí)，此時(shí)可通過視窗中心位置(λ, φ)、當(dāng)前比例尺scale、視窗寬度cw及視窗高度ch，反求出平移量(dx, dy)，最后才能完成地理坐標(biāo)到屏幕坐標(biāo)的轉(zhuǎn)換：
dx = cw / 2 - λ * 60 * MinPP / scaledy = - ch / 2 - MP(φ) * MinPP / scale
- 綜合如上公式可知：
x“ = λ * 60 * MinPP / scale + dxy“ = -MP(φ) * MinPP / scale - dyλ = (x“ - dx) * scale / 60 / MinPPφ = AMP((-y - dy) * scale / MinPP)
地理相關(guān)的類
新建類S57Pos2D存儲(chǔ)經(jīng)緯度，新建類S57Pos3D存儲(chǔ)經(jīng)緯度及水深點(diǎn)信息。
public class S57Pos2D
{
public S57Pos2D(double lat, double lon)
{
Longitude = lon;
Latitude = lat;
}
public double Longitude;
public double Latitude;
public double MercatorLongitude
{
get
{
return Longitude * 60 * GeoTools.MinPP;
}
}
public double MercatorLatitude
{
get
{
var phi = GeoTools.MP(Latitude);
return phi * GeoTools.MinPP;
}
}
public override string ToString()
{
return $"{GeoTools.LatLonToString(Latitude, true)} {GeoTools.LatLonToString(Longitude, false)}";
}
}
public class S57Pos3D : S57Pos2D
{
public S57Pos3D(double lat, double lon, double depth) : base(lat, lon)
{
Depth = depth;
}
public double Depth { get; set; }
}
了解了屏幕坐標(biāo)與地理坐標(biāo)的轉(zhuǎn)換公式后，下面新建工具類GeoTools，并編寫代碼：
public static class GeoTools
{
//每海里像素?cái)?shù) MinPP 每英尺像素?cái)?shù)Dpi 一般電腦值為 96 一英尺 = 25.4 毫米
public const double MinPP = 1852 * 100 * 96 / 2.54;
//地球圓球體的緯度漸長(zhǎng)率
public static double MP(double latitude)
{
return 3437.746771 * Math.Log(Math.Tan(Math.PI / 4 + latitude * Math.PI / 360));
}
//反函數(shù) 南北自行判斷
public static double AMP(double mp)
{
return 360 * (Math.Atan(Math.Pow(Math.E, mp / 3437.746771)) - Math.PI / 4) / Math.PI;
}
//地理坐標(biāo)轉(zhuǎn)成屏幕坐標(biāo)
//dx dy 為偏移量 由左上角或中心點(diǎn)經(jīng)緯度確定
public static Point GeoPositionToScreenPoint(S57Pos2D pos, uint scale, double dx, double dy)
{
var pmy = pos.MercatorLatitude / scale;
var pmx = pos.MercatorLongitude / scale;
return new Point((int)(pmx + dx), (int)(-(pmy + dy)));
}
//屏幕坐標(biāo)轉(zhuǎn)成地理坐標(biāo)
public static S57Pos2D ScreenPointToGeoPosition(int x, int y, uint scale, double dx, double dy)
{
var lon = (double)(x - dx) * scale / (60 * MinPP);
var lat = AMP((double)(-y - dy) * scale / MinPP);
return new S57Pos2D(lat, lon);
}
//地理緯度轉(zhuǎn)成屏幕縱坐標(biāo)
public static int LatitudeToScreenPoint(double lat, uint scale, double dx, double dy)
{
var pmy = MP(lat) * MinPP / scale;
return (int)(-(pmy + dy));
}
//地理經(jīng)度轉(zhuǎn)成屏幕橫坐標(biāo)
public static int LongitudeToScreenPoint(double lon, uint scale, double dx, double dy)
{
var pmx = lon * 60 * MinPP / scale;
return (int)(pmx + dx);
}
//根據(jù)視窗中心經(jīng)度求橫向平移量
public static int GetOffsetX(double centerLongitude, uint scale, int cw)
{
return cw / 2 - (int)(centerLongitude * 60 * MinPP / scale);
}
//根據(jù)視窗中心緯度求縱向平移量
public static int GetOffsetY(double centerLatitude, uint scale, int ch)
{
return -ch / 2 - (int)(MP(centerLatitude) * MinPP / scale);
}
//緯度 轉(zhuǎn)化成 字符串 如：+38 = 38°0.000'N
public static string LatitudeToString(double latitude, int round = 3)
{
return LatLonToString(latitude, true, 3);
}
//經(jīng)度 轉(zhuǎn)化成 字符串 如：+120 = 120°0.000'E
public static string LongitudeToString(double longitude, int round = 3)
{
return LatLonToString(longitude, false, 3);
}
//經(jīng)緯度 轉(zhuǎn)化成 字符串
public static string LatLonToString(double latlon, bool isLat, int round = 3)
{
var d = Math.Abs(latlon);
var m = (d - (int)d) * 60;
if (m >= 59.999)
{
m = 0;
d += 1;
}
var suffix = "";
if (isLat) suffix = latlon >= 0 ? "N" : "S";
else suffix = latlon >= 0 ? "E" : "W";
if (round == 3)
return $"{(int)d}°{m:0.000}'{suffix}";
else
return $"{(int)d}°{m:0.00}'{suffix}";
}
}
