import { Cesium } from '../../namespace'
import { interPointsByDistance } from './Turf'
/**
* 空间位置相关方法集
* @module Coordinate
*/
/**
* 根据像素px值拾取位置点
* @method
* @param {Object} viewer 地图场景
* @param {Cartesian2} px 屏幕坐标
* @returns {Cartesian3 | null} 位置点笛卡尔坐标
*/
export const getCatesian3FromPX = (viewer, px) => {
const scene = viewer.scene
let picks = scene.drillPick(px)
let cartesian = undefined
let isOn3dtiles = false
for (let i in picks) {
let pick = picks[i]
const primitive = pick && pick.primitive
if (
primitive instanceof Cesium.Cesium3DTileFeature ||
primitive instanceof Cesium.Cesium3DTileset ||
primitive instanceof Cesium.Model
) {
isOn3dtiles = true
}
if (isOn3dtiles) {
scene.pick(px)
cartesian = scene.pickPosition(px)
if (cartesian) {
let cartographic = Cesium.Cartographic.fromCartesian(cartesian)
let x = Cesium.Math.toDegrees(cartographic.longitude),
y = Cesium.Math.toDegrees(cartographic.latitude),
z = cartographic.height
cartesian = transformWGS84ToCartesian({ x, y, z })
}
}
}
let isOnTerrain = false // 地形
let boolTerrain =
viewer.terrainProvider instanceof Cesium.EllipsoidTerrainProvider
if (!isOn3dtiles && !boolTerrain) {
let ray = scene.camera.getPickRay(px)
if (!ray) return null
cartesian = scene.globe.pick(ray, scene)
isOnTerrain = true
}
// 地球
if (!isOn3dtiles && !isOnTerrain && boolTerrain) {
cartesian = scene.camera.pickEllipsoid(px, scene.globe.ellipsoid)
}
if (cartesian) {
let position = transformCartesianToWGS84(cartesian)
if (position.z < 0) {
position.z = 0.1
cartesian = transformWGS84ToCartesian(position)
}
return cartesian
}
return null
}
/**
* WGS84坐标转笛卡尔坐标
* @method
* @param {DegreePosZ} position WGS84坐标
* @returns {Cartesian3} 笛卡尔坐标
*/
export const transformWGS84ToCartesian = (position) => {
return position
? Cesium.Cartesian3.fromDegrees(
position.x,
position.y,
position.z,
Cesium.Ellipsoid.WGS84
)
: Cesium.Cartesian3.ZERO
}
/**
* 笛卡尔坐标转WGS84
* @method
* @param {Cartesian3} cartesian3 笛卡尔坐标
* @return {DegreePosZ} WGS84坐标
*/
export const transformCartesianToWGS84 = (cartesian3) => {
let cartographic = Cesium.Ellipsoid.WGS84.cartesianToCartographic(cartesian3)
return {
x: Cesium.Math.toDegrees(cartographic.longitude),
y: Cesium.Math.toDegrees(cartographic.latitude),
z: cartographic.height,
}
}
/**
* 笛卡尔坐标点集合转WGS84点集合
* @method
* @param {Array<Cartesian3>} cartesianArr 笛卡尔坐标点集合
* @returns {Array<DegreePosZ>} WGS84坐标点集合
*/
export const transformCartesianArrayToWGS84Array = (cartesianArr) => {
return cartesianArr
? cartesianArr.map((item) => {
return transformCartesianToWGS84(item)
})
: []
}
/**
* WGS84坐标转弧度坐标
* @method
* @param {DegreePosZ} position WGS84坐标点
* @return {Cartographic} 弧度坐标点
*/
export const transformWGS84ToCartographic = (position) => {
return position
? Cesium.Cartographic.fromDegrees(position.x, position.y, position.z)
: Cesium.Cartographic.ZERO
}
/**
* 线段插值点（地表）
* @method
* @param {Object} [viewer] 地图场景
* @param {String} [type] 插值类型，project-投影点，高度为0；onground-地表点
* @param {Array<Cartesian3>} [positions] 待插值点集合，笛卡尔坐标数组
* @param {Array<Object>} [objectsToExclude=null] 高度采集时排除的对象集合
* @returns {Array<DegreePosZ>} 经纬度点集合
*/
export const interPoints = (
viewer,
type,
positions,
objectsToExclude = null
) => {
let positionsCartographic = []
let terrainSamplePositions = []
for (let index = 0; index < positions.length; index++) {
const element = positions[index]
let ellipsoid = viewer.scene.globe.ellipsoid
let cartographic = ellipsoid.cartesianToCartographic(element)
positionsCartographic.push(cartographic)
}
for (let i = 0; i < positionsCartographic.length; i++) {
const m_Cartographic0 = positionsCartographic[i]
const m_Cartographic1 = positionsCartographic[i + 1]
let lng0 = m_Cartographic0.longitude,
lat0 = m_Cartographic0.latitude
if (m_Cartographic1) {
let lng1 = m_Cartographic1.longitude,
lat1 = m_Cartographic1.latitude
let a = Math.abs(lng0 - lng1) * 10000000
let b = Math.abs(lat0 - lat1) * 10000000
//等距采样
if (a > b) b = a
let length = parseInt((b / 2).toString())
if (length > 150) length = 150
if (length < 2) length = 2
for (let j = 0; j < length; j++) {
terrainSamplePositions.push(
new Cesium.Cartographic(
Cesium.Math.lerp(lng0, lng1, j / (length - 1)),
Cesium.Math.lerp(lat0, lat1, j / (length - 1))
)
)
}
terrainSamplePositions.pop()
} else {
terrainSamplePositions.push(m_Cartographic0)
}
}
let _result = []
for (let n = 0; n < terrainSamplePositions.length; n++) {
//地理坐标（弧度）转经纬度坐标
let curCartographic = terrainSamplePositions[n]
const z =
type === 'project'
? 0
: viewer.scene.sampleHeight(curCartographic, objectsToExclude)
const x = (curCartographic.longitude / Math.PI) * 180
const y = (curCartographic.latitude / Math.PI) * 180
_result.push({ x, y, z })
}
return _result
}
/**
* 计算多边形质点
* @method
* @param {Array<Cartesian3>} positions 多边形节点，笛卡尔坐标集合
* @returns {Array<Number>} 多边形质点坐标值 => [x,y]
*/
export const getPolygonCentroid = (positions) => {
let centroidX = 0,
centroidY = 0,
signedArea = 0
for (let index = 0; index < positions.length; index++) {
const pos1 = transformCartesianToWGS84(positions[index])
const pos2 = transformCartesianToWGS84(
positions[(index + 1) % positions.length]
)
const crossProduct = pos1.x * pos2.y - pos2.x * pos1.y
signedArea += crossProduct // 计算多边形的有向面积
centroidX += (pos1.x + pos2.x) * crossProduct
centroidY += (pos1.y + pos2.y) * crossProduct
}
signedArea *= 0.5
centroidX /= 6 * signedArea
centroidY /= 6 * signedArea
return [centroidX, centroidY]
}
/**
* 弧度转换为角度
* @method
* @param {Number} radians 弧度值
* @returns {Number} 角度值
*/
export const radiansToDegrees = (radians) => {
const degrees = radians % (2 * Math.PI)
return (degrees * 180) / Math.PI < 0
? 360 + (degrees * 180) / Math.PI
: (degrees * 180) / Math.PI
}
/**
* 角度转弧度
* @method
* @param {Number} degrees 角度值
* @returns {Number} 弧度值
*/
export const degreesToRadians = (degrees) => {
const radians = Cesium.Math.toRadians(degrees)
return radians
}
/**
* 已知点根据角度和距离求取另一点坐标（二维），注：WGS84坐标系
* @method
* @param {Number} lng 经度X
* @param {Number} lat 纬度Y
* @param {Number} angle 角度 (0~360度)
* @param {Number} distance 距离 (米)
* @returns {DegreePos} 目标点坐标
*/
export const getPointByAngleDistance = (lng, lat, angle, distance) => {
let a = 6378137 // 赤道半径
let b = 6356752.3142 // 短半径
let f = 1 / 298.257223563 // 扁率
let alpha1 = angle * (Math.PI / 180)
let sinAlpha1 = Math.sin(alpha1)
let cosAlpha1 = Math.cos(alpha1)
let tanU1 = (1 - f) * Math.tan(lat * (Math.PI / 180))
let cosU1 = 1 / Math.sqrt(1 + tanU1 * tanU1),
sinU1 = tanU1 * cosU1
let sigma1 = Math.atan2(tanU1, cosAlpha1)
let sinAlpha = cosU1 * sinAlpha1
let cosSqAlpha = 1 - sinAlpha * sinAlpha
let uSq = (cosSqAlpha * (a * a - b * b)) / (b * b)
let A = 1 + (uSq / 16384) * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
let B = (uSq / 1024) * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
let sigma = distance / (b * A),
sigmaP = 2 * Math.PI
let cos2SigmaM = 0,
sinSigma = 0,
cosSigma = 0
while (Math.abs(sigma - sigmaP) > 1e-12) {
cos2SigmaM = Math.cos(2 * sigma1 + sigma)
sinSigma = Math.sin(sigma)
cosSigma = Math.cos(sigma)
let deltaSigma =
B *
sinSigma *
(cos2SigmaM +
(B / 4) *
(cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
(B / 6) *
cos2SigmaM *
(-3 + 4 * sinSigma * sinSigma) *
(-3 + 4 * cos2SigmaM * cos2SigmaM)))
sigmaP = sigma
sigma = distance / (b * A) + deltaSigma
}
let tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
let lat2 = Math.atan2(
sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
(1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp)
)
let lambda = Math.atan2(
sinSigma * sinAlpha1,
cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1
)
let C = (f / 16) * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
let L =
lambda -
(1 - C) *
f *
sinAlpha *
(sigma +
C *
sinSigma *
(cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
return {
x: lng + L * (180 / Math.PI),
y: lat2 * (180 / Math.PI),
}
}
/**
* 区域内生成三角网
* @method
* @param {Array<Cartesian3>} positions 三角网生成的范围点集合，笛卡尔坐标
* @param {Number} s - 精度，一般为 4/经度or纬度差
* @returns {Array<DegreePosZ[]>} 三角网集合，WGS84坐标
*/
export const subTriangle = (positions, s) => {
let granularity = Math.PI / Math.pow(2, 11)
granularity = granularity / s
const polygonGeometry = Cesium.PolygonGeometry.fromPositions({
positions: positions,
vertexFormat: Cesium.PerInstanceColorAppearance.FLAT_VERTEX_FORMAT,
granularity: granularity,
})
const geom = Cesium.PolygonGeometry.createGeometry(polygonGeometry)
const triangles = []
for (let i = 0; i < geom.indices.length; i += 3) {
const i0 = geom?.indices[i]
const i1 = geom?.indices[i + 1]
const i2 = geom?.indices[i + 2]
const subTrianglePositions = geom?.attributes.position.values
const pos0 = {
x: subTrianglePositions[i0 * 3],
y: subTrianglePositions[i0 * 3 + 1],
z: subTrianglePositions[i0 * 3 + 2],
}
const pos1 = {
x: subTrianglePositions[i1 * 3],
y: subTrianglePositions[i1 * 3 + 1],
z: subTrianglePositions[i1 * 3 + 2],
}
const pos2 = {
x: subTrianglePositions[i2 * 3],
y: subTrianglePositions[i2 * 3 + 1],
z: subTrianglePositions[i2 * 3 + 2],
}
let points = [pos0, pos1, pos2]
triangles.push(points)
}
return triangles
}
/**
* 自动进行插值并计算获取区域内的极值点（最高点和最低点）
* @method
* @param {Object} [viewer] 地图场景
* @param {Array<Cartesian3>} [positions] 坐标点集合，笛卡尔坐标
* @param {Number} [num=25] 精度按照经纬度1°大约100公里，默认按25等分来计算
* @returns {Object} 最高和最低点的高度值和位置 => {minH:0,maxH:0,minPos:{x,y}, maxPos:{x,y}}
*/
export const getRegionExtreme = (viewer, positions, num = 25) => {
let minX = 1000,
maxX = -1000,
minY = 1000,
maxY = -1000
const points = positions.map((pos) => {
const { x, y } = transformCartesianToWGS84(pos)
minX = Math.min(x, minX)
maxX = Math.max(x, maxX)
minY = Math.min(y, minY)
maxY = Math.max(y, maxY)
return [x, y]
})
points.push(points[0])
const box = [minX, minY, maxX, maxY]
const _d = Math.min(maxX - minX, maxY - minY)
// 精度按照经纬度1°大约100公里，默认按25等分来计算
const grid = interPointsByDistance(box, points, (_d * 100) / num)
const features = grid.features
const posList = []
const cartographics = []
for (let index = 0; index < features.length; index++) {
const feature = features[index]
const coord = feature.geometry.coordinates
posList.push({ x: coord[0], y: coord[1] })
cartographics.push(Cesium.Cartographic.fromDegrees(coord[0], coord[1]))
}
let minH = 1000000
let maxH = -1000000
let minPos = posList[0]
let maxPos = posList[0]
for (let i = 0; i < cartographics.length; i++) {
const graphic = cartographics[i]
let h = viewer.scene.sampleHeightSupported
? viewer.scene.sampleHeight(graphic)
: viewer.scene.globe.getHeight(graphic)
if (h < minH) {
minH = h
minPos = posList[i]
}
if (h > maxH) {
maxH = h
maxPos = posList[i]
}
}
return { minH, maxH, minPos, maxPos }
}
/**
* 二维状态下计算两点对于正北方向的朝向角度，范围[0,360]
* @method
* @param {DegreePos} start 起始点，经纬度点
* @param {DegreePos} end 目标点，经纬度点
* @returns {Number} 角度值
*/
export const calculateAngle = (start, end) => {
let rad = Math.PI / 180,
lat1 = start.y * rad,
lat2 = end.y * rad,
lon1 = start.x * rad,
lon2 = end.x * rad
const a = Math.sin(lon2 - lon1) * Math.cos(lat2)
const b =
Math.cos(lat1) * Math.sin(lat2) -
Math.sin(lat1) * Math.cos(lat2) * Math.cos(lon2 - lon1)
return radiansToDegrees(Math.atan2(a, b))
}
/**
* 经纬度坐标数组转笛卡尔坐标数组
* @method
* @param {Array<Number[]>} positions 经纬度坐标数组，如：[[120,34],[121,35]]
* @returns {Array<Cartesian3>} 笛卡尔坐标数组
*/
export const PosFromDegreeArray = (positions) => {
let result = []
for (let index = 0; index < positions.length; index++) {
const pos = positions[index]
result.push(pos[0], pos[1])
}
const points = Cesium.Cartesian3.fromDegreesArray(result)
return points
}
/**
* 根据经纬度数组获取地表位置点数组
* @method
* @param {Object} terrain 当前场景地形对象
* @param {Array<Number[]>} points 经纬度数组，如：[[120,34],[121,35]]
* @returns {Promise<Array<Cartesian3>>} 异步返回，地表点集合，笛卡尔坐标数组
*/
export const PosOnTerrainFromDegreeArray = async (terrain, points) => {
const positions = points.map((point) => {
return Cesium.Cartographic.fromDegrees(point[0], point[1])
})
const updatedPositions = await Cesium.sampleTerrainMostDetailed(
terrain,
positions
)
const result = updatedPositions.map((pos) => {
return Cesium.Cartographic.toCartesian(pos)
})
return result
}
/**
* 根据经纬度获取地表位置点
* @method
* @param {Object} terrain 当前场景地形对象
* @param {Array<Number>} position 经纬度坐标，如[120,35]
* @returns {Promise<Cartesian3>} 异步返回：地表坐标点，笛卡尔坐标
*/
export const PosOnTerrainFromDegree = async (terrain, position) => {
const cartographic = Cesium.Cartographic.fromDegrees(position[0], position[1])
const updatedPositions = await Cesium.sampleTerrainMostDetailed(terrain, [
cartographic,
])
const result = Cesium.Cartographic.toCartesian(updatedPositions[0])
return result
}
/**
* 获取指定经纬度点地表高度
* @method
* @param {Object} terrain 当前场景地形对象
* @param {Array<Number>} position 经纬度坐标，如[120,35]
* @returns {Promise<Number>} 异步返回：当前点的地表高度值
*/
export const getHeight = async (terrain, position) => {
const cartographic = Cesium.Cartographic.fromDegrees(position[0], position[1])
const updatedPositions = await Cesium.sampleTerrainMostDetailed(terrain, [
cartographic,
])
const height = updatedPositions[0].height
return height
}
/**
* 获取经纬度点集高度
* @method
* @param {Object} terrain 当前场景地形对象
* @param {Array<Number[]>} points 经纬度坐标数组，如[[120,35],[121,36]]
* @returns {Promise<Array<Number>>} 异步返回：高度值数组
*/
export const getHeightArray = async (terrain, points) => {
const positions = points.map((point) => {
return Cesium.Cartographic.fromDegrees(point[0], point[1])
})
const updatedPositions = await Cesium.sampleTerrainMostDetailed(
terrain,
positions
)
const result = updatedPositions.map((pos) => {
return pos.height
})
return result
}
/**
* 按倍数调整笛卡尔坐标点的高度，并返回更新后的笛卡尔坐标
* @method
* @param {Cartesian3} [pos] 笛卡尔坐标
* @param {Number} [scale] 高度放大倍数
* @param {Number} [min=0] 限定最小高度值，默认值为0，设置后，如果计算后的高度值小于min值，则高度直接取min值
* @returns {Cartesian3} 更新后的笛卡尔坐标
*/
export const scaleHeight = (pos, scale, min = 0) => {
const cartographic = Cesium.Cartographic.fromCartesian(pos)
let { longitude, latitude, height } = cartographic
const h = Math.max(height * scale, min)
const cartographic1 = new Cesium.Cartographic(longitude, latitude, h)
return Cesium.Cartographic.toCartesian(cartographic1)
}
/**
* 调整点高度为指定高度，并返回更新后的笛卡尔坐标
* @method
* @param {Cartesian3} pos 笛卡尔坐标
* @param {Number} height 高度值
* @returns {Cartesian3} 更新后的笛卡尔坐标
*/
export const changeHeight = (pos, height) => {
const cartographic = Cesium.Cartographic.fromCartesian(pos)
let { longitude, latitude } = cartographic
const cartographic1 = new Cesium.Cartographic(longitude, latitude, height)
return Cesium.Cartographic.toCartesian(cartographic1)
}
/**
* 经纬度四至转弧度四至
* @method
* @param {Array<Number>} extent 经纬度四至，如：[112,23,120,30]
* @returns {Array<Number>} 转换后的弧度四至
*/
export const extentToRadians = (extent) => {
return extent.map((val) => degreesToRadians(val))
}
/**
* 弧度四至转经纬度四至
* @method
* @param {Array<Number>} extent 弧度四至
* @returns {Array<Number>} 转换后的经纬度四至
*/
export const extentToDegrees = (extent) => {
return extent.map((val) => radiansToDegrees(val))
}
/**
* 判断点在地球的正面还是反面
* @method
* @param {Cartesian3} position 点位，笛卡尔坐标
* @param {Cartesian3} cameraPosition 相机位置，笛卡尔坐标
* @returns true-正面，false-反面
*/
export const isVisible = (position, cameraPosition) => {
let cameraOccluder = new Cesium.EllipsoidalOccluder(
Cesium.Ellipsoid.WGS84,
cameraPosition
)
let viewerVisible = cameraOccluder.isPointVisible(position)
return viewerVisible
}