地球表面是一个不规则的曲面，而海图却是平面的，从不规则曲面到平面的过程被称为投影，这种投影必然在不同的地点会有不同的变形。传统航海确定方向的主要工具是指南针，如果一个人照着指南针给定的固定角度航行，肯定被认为是走的直线，这种”直线“的性质也应该体现在海图上；同时，现实世界中前后左右的物体被投影到海图其相对方位也应该保持不变，即“等角”性质。为了满足航海人员对海图的“直线”和“等角”的性质，荷兰地图学家墨卡托在1569年拟定了一种投影规则 -- “等角正圆柱投影”。因此，利用等角正圆柱投影的海图被称作墨卡托海图，世界上95%以上的海图都是墨卡托海图。
由上图可知墨卡托海图上：等经差在海图上等间距；等纬差在海图上不等间距，且随纬度增加，间距逐渐加大。
恒向线
比如，设定角度60度，一直沿着这个角度航行，真实航行轨迹是下图的曲线，而非直线。该曲线被称之为恒向线，是一条与所有经线相交成恒定角度的、具有双重曲率的球面螺旋线，它逐渐趋向地极，但永远达不到地极。
以地球为圆球体简化计算可知恒向线的方程如下：
对恒向线上的微分三角形DM1M2进行积分，即可得恒向线公式：
式中：起点(φ1, λ1)，终点(φ2, λ2)，航向角C
恒向线应用场景：
- 已知起点、终点，求解两点间恒向线航向，及线上所有点的坐标；
- 已知起点、航向，求解线上其他点的坐标
赤道里
赤道上经度1'所对应的弧长，各国采用的地球椭圆体参数不同而略有不同，约为1855.36米。简化计算中，赤道里长度 = 海里长度 = 1852米。地球赤道周长 = 360 * 60 * 赤道里
等角性质 -- 纬度渐长率
等角意味着每一投影点在任一方向上的变形是相同的。根据理论，若一点在经线方向与纬线方向变形一致，那该点在任何方向的变形都一致，即满足等角。因此，需要了解一个新的概念：纬度渐长率。纬度渐长率是指在墨卡托海图上任一纬度线到赤道的子午线的长度与图上1赤道里的比值。
对于地球的微分区域ABCD，B点在经线方向与纬线方向投影变形一致，即可推出纬度渐长率公式。
式中：a为地球椭圆体长半轴长度。纬度渐长率大小只与纬度有关，与墨卡托海图比例尺无关。可以证明，恒向线在墨卡托海图上是直线。
海图绘制中，需先确定1赤道里的图长，即海图绘制单位。然后计算任意两纬度间的纬度渐长率之差，乘以海图绘制单位，即可得两纬线间的图长。绘制出各经纬线图网后，将地理物标填充进图网即可。
电子海图投影步骤
- 一英尺=25.4毫米，一般电脑屏幕每英尺像素数Dpi=96，每毫米像素数=96/25.4；
- 赤道里的长度为1852米，每赤道里像素数
MinPP=1852*100*10*96/25.4； - 以经度0度、纬度0度为原点，东北方向为正，可将地球按照墨卡托投影成平面；
- 地理坐标(λ, φ)与平面坐标(x, y)间可由下列公式互转：
x = λ * 60 * MinPPy = MP(φ) * MinPPφ = AMP(y / MinPP)
MP为纬度渐长率函数；AMP为纬度渐长率的反函数。
- 平面坐标(x, y)右上为正，而屏幕坐标(x', y')原点在左上角，左下为正，其互转公式如下：
x = x'y = -y'
- 电子海图可以随意先缩放scale（比例尺形如：1:10000，为方便计算与表示，使scale=10000），再平移(dx, dy)，因此还需要进一步将原始屏幕坐标(x', y')转换为显示坐标(x“, y”):
x“ = x' / scale + dxy” = y' / scale - dy
- 电子海图系统指定屏幕范围内的中心位置，或用户跳转到指定位置时，此时可通过视窗中心位置(λ, φ)、当前比例尺scale、视窗宽度cw及视窗高度ch，反求出平移量(dx, dy)，最后才能完成地理坐标到屏幕坐标的转换：
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)
地理相关的类
新建类S57Pos2D存储经纬度，新建类S57Pos3D存储经纬度及水深点信息。
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; }
}
了解了屏幕坐标与地理坐标的转换公式后，下面新建工具类GeoTools，并编写代码：
public static class GeoTools
{
//每海里像素数 MinPP 每英尺像素数Dpi 一般电脑值为 96 一英尺 = 25.4 毫米
public const double MinPP = 1852 * 100 * 96 / 2.54;
//地球圆球体的纬度渐长率
public static double MP(double latitude)
{
return 3437.746771 * Math.Log(Math.Tan(Math.PI / 4 + latitude * Math.PI / 360));
}
//反函数 南北自行判断
public static double AMP(double mp)
{
return 360 * (Math.Atan(Math.Pow(Math.E, mp / 3437.746771)) - Math.PI / 4) / Math.PI;
}
//地理坐标转成屏幕坐标
//dx dy 为偏移量 由左上角或中心点经纬度确定
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)));
}
//屏幕坐标转成地理坐标
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);
}
//地理纬度转成屏幕纵坐标
public static int LatitudeToScreenPoint(double lat, uint scale, double dx, double dy)
{
var pmy = MP(lat) * MinPP / scale;
return (int)(-(pmy + dy));
}
//地理经度转成屏幕横坐标
public static int LongitudeToScreenPoint(double lon, uint scale, double dx, double dy)
{
var pmx = lon * 60 * MinPP / scale;
return (int)(pmx + dx);
}
//根据视窗中心经度求横向平移量
public static int GetOffsetX(double centerLongitude, uint scale, int cw)
{
return cw / 2 - (int)(centerLongitude * 60 * MinPP / scale);
}
//根据视窗中心纬度求纵向平移量
public static int GetOffsetY(double centerLatitude, uint scale, int ch)
{
return -ch / 2 - (int)(MP(centerLatitude) * MinPP / scale);
}
//纬度 转化成 字符串 如：+38 = 38°0.000'N
public static string LatitudeToString(double latitude, int round = 3)
{
return LatLonToString(latitude, true, 3);
}
//经度 转化成 字符串 如：+120 = 120°0.000'E
public static string LongitudeToString(double longitude, int round = 3)
{
return LatLonToString(longitude, false, 3);
}
//经纬度 转化成 字符串
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}";
}
}
