three-geo加载地形的实际应用

Q：threeJs与GIS结合，threeJs中已加载有真实地理位置的三维模型和点线数据，需加载地形和影像数据，与模型数据、矢量数据匹配。

A：

1. 地形影像加载

threeJs-geo已基于threeJS和mapbox封装了地形和影像的加载库，基本用法：

const tgeo = new ThreeGeo({
	tokenMapbox: '', // <---- set your Mapbox API token here
});

const terrain = await tgeo.getTerrainRgb(
	[lat, lng],   // [36, 118], 经纬度
	radius,       // 2， radius of bounding circle (km) , 
	zoom);   // 14,  影像级别
scene.add(terrain);

three-geo是基于mapbox封装的，因此需要mapbox地图的token(官网申请)，此处用到的方法是getTerrainRgb，地形作为模型加载而影像是地形模型的贴图，三个参数的含义在代码注释中已经写的很清楚，加载效果：

如果仅仅是看一下某个地方的地形影像，到此处完结！

但是要与真实地理数据匹配，就存在以下问题：

     1）地形位置问题：地形在threeJS加载的默认位置为（0,0,0）点， 需将地形定位到真实地理位置

     2）地形比例问题：地形单位与threejs场景单位不一致，此处获取到地形作为模型加载后，比例缩小(缩小规则未搞清)，如加载半径2km范围的地形，获取到的模型大小是长宽均为1.27，明显不对

 

2. 地形位置问题

threeJS使用的是世界坐标，在该问题中，获取地形的经纬度不是凭空捏造的

     1）拾取坐标点：在threeJS场景中拾取模型或矢量数据上的某一个点p1(a, b, c)，作为要加载和显示地形的地方

     2）p1为世界坐标，在GIS的概念中，也是投影坐标，通过proj4.js库将p1转为经纬坐标

let lonlat = proj4(poi1Projection, proj4('EPSG:4326'), p1);

p1为地形的加载位置, lonlat为获取地形&影像的经纬度参数

terrain.position.set(p1.x, p1.y, 0)

 

3. 地形比例问题

首先要认识一个问题，three-geo获取到的地形模型的真实大小是否与radius设置的值相符？并不。

假设参数radius设为2，获取到的地形&影像范围长宽并非为4km，而是做了相交处理，以2km为半径获取相交的影像瓦片及瓦片对应的地形范围，如：

radius=0.5， 瓦片数量：1x2，地形模型大小1.27

radius=2，瓦片梳理：2x2，地形模型大小1.27

所以地形初始化加载时的大小与radius设置的半径范围关联不大，加载的地形&影像真实范围屈居于相交的影像瓦片数量。

（以上问题可阅读threeJS-geo源码）

解决方式： 

    1）计算 R(纬度圈所在半径) = 地球半径 * cos(lat)

         (注：  纬度lat在前文算出的lonlat中已获得，JS中Math.cos方法用的是弧度而不是角度， 需把 lat 转为弧度 lat * PI /180)

    2）计算瓦片分辨率 resolution = 2*PI*R/(pixel*Math.pow(2,zoom)) 

         (注： R为所在纬度圈的半径，而非地球半径，dpi为影像瓦片像素，mapbox的瓦片像素为512， zoom为缩放等级)   

    3）计算xy瓦片数量 number： 用的办法比较土，观察terrain模型，它本质上是由瓦片组成，每个瓦片是模型的一个子节点，所以terrain.children.length 其实就是瓦片的数量, 而且存在一个特点： terrain.children.length = numberX * numberY (x方向瓦片数量和y方向瓦片数量， xy很多时候是相等的)。 然后初中数学就可以算出来numberX 和 numberY 了

            let xNum, yNum
			let sqrtSide = Math.sqrt( terrain.children.length )
			if(sqrtSide%1 === 0){
				xNum = yNum = sqrtSide
			}else{
				let n = boxSize.x/boxSize.y // terrain.children.length = (yNum * n) * yNum = n * Math.pow(yNum , 2)
				yNum = Math.round(Math.sqrt(terrain.children.length/n)) // 瓦片数量
				xNum = Math.round(yNum * n)
			}

    4）地形真实大小： width = numberX * resolution * pixel, 

                             height = numberY * resolution * pixel

    5）获取地形缩放比例尺： scale = width/modelWidth

 

以上，设置好地形模型的位置，修改比例尺，那么结果就是对的了，当然仍然存在一定的误差，radius越大，加载的面积越大误差越大。

 

关键代码：

var targetProjection = `PROJCS["CGCS2000 / 3-degree Gauss-Kruger CM 120E",
    GEOGCS["China Geodetic Coordinate System 2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101,
                AUTHORITY["EPSG","1024"]],
            AUTHORITY["EPSG","1043"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4490"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",120],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","4549"]]`;

// proj4的定义方式有两种
// proj4.defs([["EPSG:4549","+proj=tmerc +lat_0=0 +lon_0=120 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"]])
// targetProjection = proj4('EPSG:4549')

		let coord = proj4(targetProjection, proj4('EPSG:4326'), [p1.x, p1.y]);
		camera.position.set(0, 0, 1.5);

		(async () => { // main
			const tgeo = new ThreeGeo({
				tokenMapbox: 'pk.eyJ1IjoibGF2b25uZSIsImEiOiJja2xodXhwN20xeHJrMm9tanQzc2s3dzR1In0.M4qxPgDJDitR7JvmXrUXdQ', // <---- set your Mapbox API token here
			});

			if (tgeo.tokenMapbox === '********') {
				const warning = 'Please set your Mapbox API token in ThreeGeo constructor.';
				alert(warning);
				throw warning;
			}


			let zoomLevel = 16
			const terrain = await tgeo.getTerrainRgb(
					[coord[1], coord[0]], // [lat, lng]
					7,               // radius of bounding circle (km)
					zoomLevel);      // zoom resolution


			// resolution = 2*PI*R/(256*Math.pow(2,zoom)) 256为像素。根据不同地图瓦片像素进行修改
 
            let R = 20037508.34 * Math.cos(coord[0] * 3.1415926 /180)
			let resolution = R * 2 / (512 * Math.pow(2 , zoomLevel))
		
			let box = new THREE.Box3()
			box.setFromObject(terrain)
			let boxSize = box.size()
			let xNum, yNum
			let sqrtSide = Math.sqrt( terrain.children.length )
			if(sqrtSide%1 === 0){
				xNum = yNum = sqrtSide
			}else{
				let n = boxSize.x/boxSize.y // terrain.children.length = (yNum * n) * yNum = n * Math.pow(yNum , 2)
				yNum = Math.round(Math.sqrt(terrain.children.length/n)) // 瓦片数量
				xNum = Math.round(yNum * n)
			}
			let width = resolution * 512 * xNum
			let height = resolution * 512 * yNum
			terrain.scale.multiplyScalar(width/boxSize.x)
            terrain.position.set(p1.x, p1.y, 0);
			scene.add(terrain);

			renderer.render(scene, camera);
		})();