地球表面是一個(gè)不規(guī)則的曲面，而海圖卻是平面的挺尿，從不規(guī)則曲面到平面的過(guò)程被稱(chēng)為投影奏黑，這種投影必然在不同的地點(diǎn)會(huì)有不同的變形。傳統(tǒng)航海確定方向的主要工具是指南針编矾，如果一個(gè)人照著指南針給定的固定角度航行熟史，肯定被認(rèn)為是走的直線(xiàn)，這種”直線(xiàn)“的性質(zhì)也應(yīng)該體現(xiàn)在海圖上窄俏；同時(shí)蹂匹，現(xiàn)實(shí)世界中前后左右的物體被投影到海圖其相對(duì)方位也應(yīng)該保持不變，即“等角”性質(zhì)裆操。為了滿(mǎn)足航海人員對(duì)海圖的“直線(xiàn)”和“等角”的性質(zhì)怒详，荷蘭地圖學(xué)家墨卡托在1569年擬定了一種投影規(guī)則 -- “等角正圓柱投影”。因此踪区，利用等角正圓柱投影的海圖被稱(chēng)作墨卡托海圖，世界上95%以上的海圖都是墨卡托海圖吊骤。
由上圖可知墨卡托海圖上：等經(jīng)差在海圖上等間距缎岗；等緯差在海圖上不等間距，且隨緯度增加白粉，間距逐漸加大传泊。
恒向線(xiàn)
比如，設(shè)定角度60度鸭巴，一直沿著這個(gè)角度航行眷细，真實(shí)航行軌跡是下圖的曲線(xiàn)，而非直線(xiàn)鹃祖。該曲線(xiàn)被稱(chēng)之為恒向線(xiàn)溪椎，是一條與所有經(jīng)線(xiàn)相交成恒定角度的、具有雙重曲率的球面螺旋線(xiàn)，它逐漸趨向地極校读，但永遠(yuǎn)達(dá)不到地極沼侣。
以地球?yàn)閳A球體簡(jiǎn)化計(jì)算可知恒向線(xiàn)的方程如下：
對(duì)恒向線(xiàn)上的微分三角形DM1M2進(jìn)行積分，即可得恒向線(xiàn)公式：
式中：起點(diǎn)(φ1, λ1)歉秫，終點(diǎn)(φ2, λ2)蛾洛，航向角C
恒向線(xiàn)應(yīng)用場(chǎng)景：
- 已知起點(diǎn)、終點(diǎn)雁芙，求解兩點(diǎn)間恒向線(xiàn)航向轧膘，及線(xiàn)上所有點(diǎn)的坐標(biāo)；
- 已知起點(diǎn)兔甘、航向谎碍，求解線(xià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)線(xiàn)方向與緯線(xiàn)方向變形一致仙蛉，那該點(diǎn)在任何方向的變形都一致笋敞，即滿(mǎn)足等角。因此荠瘪，需要了解一個(gè)新的概念：緯度漸長(zhǎng)率夯巷。緯度漸長(zhǎng)率是指在墨卡托海圖上任一緯度線(xiàn)到赤道的子午線(xiàn)的長(zhǎng)度與圖上1赤道里的比值。
對(duì)于地球的微分區(qū)域ABCD哀墓，B點(diǎn)在經(jīng)線(xiàn)方向與緯線(xiàn)方向投影變形一致趁餐，即可推出緯度漸長(zhǎng)率公式。
式中：a為地球橢圓體長(zhǎng)半軸長(zhǎng)度篮绰。緯度漸長(zhǎng)率大小只與緯度有關(guān)后雷，與墨卡托海圖比例尺無(wú)關(guān)。可以證明吠各，恒向線(xiàn)在墨卡托海圖上是直線(xiàn)臀突。
海圖繪制中，需先確定1赤道里的圖長(zhǎng)贾漏，即海圖繪制單位候学。然后計(jì)算任意兩緯度間的緯度漸長(zhǎng)率之差，乘以海圖繪制單位纵散，即可得兩緯線(xiàn)間的圖長(zhǎng)梳码。繪制出各經(jīng)緯線(xiàn)圖網(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)的中心位置，或用戶(hù)跳轉(zhuǎn)到指定位置時(shí)炬守，此時(shí)可通過(guò)視窗中心位置(λ, φ)牧嫉、當(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)的類(lèi)
新建類(lèi)S57Pos2D存儲(chǔ)經(jīng)緯度，新建類(lèi)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)換公式后辽剧，下面新建工具類(lèi)GeoTools，并編寫(xiě)代碼：
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}";
}
}
