Leaflet.js/Proj4.js自定义地图投影

Leaflet官网上就能下载到，国内几家做Web端GIS的底层技术大部分基于OpenLayer与Leaflet，由于Leaflet基础组件的轻便，广受青睐。

Leaflet自定义投影方法

最近在做海图的web版本，网络上提供了很多地图加载瓦片的方法。
1、主流是引入Proj4作为插件使用，如果想做循环显示，会发现这个插件对于L.Map中worldcopyjump此类的方法支持并不好。
2、不使用Proj4的方法也可以加载自定义的投影方法，当然必须是3857、3395此类基本转换的变种。太复杂的变换还是使用Proj4比较靠谱，也许效果不好，最少还是可以显示的。

Proj4

在leaflet中使用proj4定义的crs

 merc.x = 20037508.34279;
 merc.y = 15496570.73972;
 // 加载投影
 var r = new Array(19);
 for (var i = 0; i < r.length; i++) {
 	r[i] = merc.x * 2 / (Math.pow(2, i) * 256);
 }
 var crs = new L.Proj.CRS(
 	"EPSG:3395",
 	"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs",
 	{
 		resolutions: r,
 		origin: [-merc.x, merc.y],
 		bounds: L.bounds([-merc.x, -merc.y], [merc.x, merc.y])
 	});
 var map = L.map('map', {
		crs: crs  // 定义的坐标系
	});

标准的3395->https://epsg.io/3395

我定义的海图格式不是标准的3395，但是由于瓦片是公司自定义的，所以切图的经纬度边界
[Lng,Lat]->[-180,-80],[180,80]->[-20037508.34279,-15496570.73972],[20037508.34279,15496570.73972]

定义海图各层resolutions分辨率、原点、边界，加载到L.Map上就OK了。

注意：Leaflet CRS并不定义坐标中点，坐标中点默认为0.5*(Math.PI * Mercator.R)

Proj4解释一下

看一下L.Proj.CRS这个类

L.Proj.CRS = L.Class.extend({
		includes: L.CRS,

		options: {
			transformation: new L.Transformation(1, 0, -1, 0)
		},

L.Proj.CRS继承L.CRS，重新定义transformation这个option。

贴一下Leaflet L.CRS.EPSG3395 源码

var EPSG3395 = extend({}, Earth, {
	code: 'EPSG:3395',
	projection: Mercator,

	transformation: (function () {
		//上面说过了，定义原点也是不生效的
		var scale = 0.5 / (Math.PI * Mercator.R);
		return toTransformation(scale, 0.5, -scale, 0.5);
	}())
});

var Earth = extend({}, CRS, {
	wrapLng: [-180, 180],

	// Mean Earth Radius, as recommended for use by
	// the International Union of Geodesy and Geophysics,
	// see http://rosettacode.org/wiki/Haversine_formula
	R: 6371000,

	// distance between two geographical points using spherical law of cosines approximation
	distance: function (latlng1, latlng2) {
		var rad = Math.PI / 180,
		    lat1 = latlng1.lat * rad,
		    lat2 = latlng2.lat * rad,
		    a = Math.sin(lat1) * Math.sin(lat2) +
		        Math.cos(lat1) * Math.cos(lat2) * Math.cos((latlng2.lng - latlng1.lng) * rad);

		return this.R * Math.acos(Math.min(a, 1));
	}
});

好了，到Earth这层已经可以了，3395定义的球体是正球基础上定义椭球。之后将中点坐标的墨卡托坐标转换成像素坐标，
如下面的代码：

if (this.options.origin) {
	this.transformation =
		new L.Transformation(1, -this.options.origin[0], -1, this.options.origin[1]);
}

transformation这个函数。看下这个函数做什么用的

L.Transformation.prototype = {
	transform: function (point, scale) { // (Point, Number) -> Point
		return this._transform(point.clone(), scale);
	},

	// destructive transform (faster)
	_transform: function (point, scale) {
		scale = scale || 1;
		point.x = scale * (this._a * point.x + this._b);
		point.y = scale * (this._c * point.y + this._d);
		return point;
	},

	untransform: function (point, scale) {
		scale = scale || 1;
		return new L.Point(
		        (point.x / scale - this._b) / this._a,
		        (point.y / scale - this._d) / this._c);
	}
};

简单说，transform就是将Mercator的坐标转化成像素坐标，untransform就是反转。其中的Zoom参数也很好理解，不同的比例尺下同一经纬度的像素坐标一定是不一样的。

---

自定义一个CRS

我之前写了很多Proj与Leaflet中CRS的转换解释，有点过于复杂了。

当墨卡托投影出现边界不一致的情况。只需要在transform这个经纬度转换函数乘以一个系数就行了。这个系数就是墨卡托边界坐标的比例。

其实不用一定要使用Proj4，可以直接在leaflet里面在定义一个投影方式，我们能用到的投影方式无外乎Mercator与Web Mercator这两种，在leaflet基础方法中就已经定义的很完善了，我们只需要在墨卡托坐标转像素坐标的时候转换正确就可以了。

var EPSG3396 = extend({}, Earth, {
	code: 'EPSG:3396',
	projection: Mercator,

	transformation: (function () {
		var scale = 0.5 / (Math.PI * Mercator.R);
		return toTransformation(scale, 0.5, -scale, 0.5*15496570.73972/20037508.34279);
	}())
});

---

L.CRS与L.Proj.CRS

那L.CRS.EPSG3395 使用的transformation是什么呢

 toTransformation(scale, 0.5, -scale, 0.5);

而L.Proj.CRS的transformation

L.Transformation(1, -this.options.origin[0],-1, this.options.origin[1]);

这两段代码参数不一样，计算方法就不一样么？往下看就知道根本就是一回事。
先定义一下长短轴，等会方便解释，我们来算一下：
Ma = 20037508.34279;
Mb = 15496570.73972;

从L.CRS.EPSG3395转L.Proj.CRS的transform

point.x = scale * (this._a * point.x + this._b);

point.x = 256×2 ^m (1/2Ma * point.x + 1/2)

point.x = 256×2 ^m/2Ma(1* point.x + Ma)

//不同比例尺下 _scales = 1/分辨率
this._scales[i] = 1 / this.options.resolutions[i];

L.Proj.CRS的resolutions是256×2 ^m /2*Ma

point.x = 1 / this.options.resolutions[i] × (1* point.x + Ma)

看出来一样了吧。

因为两种投影方法 X轴是一样的，Y轴就有差异了
L.CRS.EPSG3395–>point.y = 256×2 ^m (1/2Ma * point.y + 1/2) 依然如此
L.Proj.CRS–>point.y = 256×2 ^m /2Ma × (1 point.y + Mb)；
所以在L.Proj.CRS是比较完美的自定义地图的方法。

---

贴一个leaflet中的Mercator计算
网上的Mercator精度太差了，贴一个leaflet定义的方法。

var Mercator = {
	R: 6378137,
	R_MINOR: 6356752.314245179,

	bounds: new Bounds([-20037508.34279, -15496570.73972], [20037508.34279, 18764656.23138]),

	project: function (latlng) {
			var d = Math.PI / 180,
				r = this.R,
				y = latlng.lat * d,
				tmp = this.R_MINOR / r,
				e = Math.sqrt(1 - tmp * tmp),
				con = e * Math.sin(y);

			var ts = Math.tan(Math.PI / 4 - y / 2) / Math.pow((1 - con) / (1 + con), e / 2);
			y = -r * Math.log(Math.max(ts, 1E-10));

			return new Point(latlng.lng * d * r, y);
	},

	unproject: function (point) {
		var d = 180 / Math.PI,
		    r = this.R,
		    tmp = this.R_MINOR / r,
		    e = Math.sqrt(1 - tmp * tmp),
		    ts = Math.exp(-point.y / r),
		    phi = Math.PI / 2 - 2 * Math.atan(ts);

		for (var i = 0, dphi = 0.1, con; i < 15 && Math.abs(dphi) > 1e-7; i++) {
			con = e * Math.sin(phi);
			con = Math.pow((1 - con) / (1 + con), e / 2);
			dphi = Math.PI / 2 - 2 * Math.atan(ts * con) - phi;
			phi += dphi;
		}

		return new LatLng(phi * d, point.x * d / r);
	}
};

---

墨卡托投影
https://blog.csdn.net/mygisforum/article/details/13295223

比例尺，分辨率，级别，及其之间的转换
https://gaoyi2009.iteye.com/blog/1271014