一、卫星轨迹
实现卫星轨迹，我这里采用的是加载czml的方式，采用cesium.CzmlDataSource API来加载卫星轨迹，话不多说，先看效果。
在这个项目中，除了卫星实时轨道，还有一个重要的功能过境预测，就是可以上传一个区域或者框选一个区域，亦或者选取一个行政边界，根据选的卫星和时间，通过计算卫星轨迹，来预测卫星何时经过此地，并且将经过此地的轨道画出来。大致的效果如下：
二、实现方法
首先需要czml文件，通过卫星轨道二根数生成卫星轨道的czml文件，文件的格式及案例如资源中czml所示。在初始化cesium地图后，初始化cesium及cesium相关参数在前文中已经讲解过了，这里就不过多赘述了，在页面初始化时，调用satelliteRoad即可加载卫星轨迹。
2.1卫星轨迹
function satelliteRoad(params, format) {
//去除卫星轨迹
viewer.dataSources._dataSources.forEach((czmlDataSource) => {
if (czmlDataSource instanceof Cesium.CzmlDataSource) {
viewer.dataSources.remove(czmlDataSource, true);
}
});
$.ajax({
type: "post",
url: build_config+ "/ems/address/addSatellite",
dataType: "json",
async: false,
contentType: "application/json;charset=utf-8",
data: JSON.stringify(params),
processData: false,
success: function(msg) {
var czml = msg.result;
console.log(czml);
var satelliteCZML = new Cesium.CzmlDataSource.load(msg.result);
viewer.dataSources.add(satelliteCZML).then(function(dataSource) {
for (let i = 0; i < dataSource.entities._entities.length; i++) {
addEntities(dataSource.entities._entities._array[i],czml[i+1].label.text);
}
});
},
error: function(jqXHR, textStatus, errorThrown) {
alert("添加轨道出错");
},
});
}
//添加卫星圆锥实体
function addEntities(satellite,satelliteName) {
var length = 625000.0;
if(satelliteName == "高分四号"){
length = 0.0;
}else{
length = 625000.0;
}
var color1 = new Cesium.Color(
satellite.path.material.color._value.red,
satellite.path.material.color._value.green,
satellite.path.material.color._value.blue
);
var property = new Cesium.SampledPositionProperty();
var projectionPositions = [];
var cylinderEntity1 = viewer.entities.add({
// position:property,
name: "cylinder",
cylinder: {
show: true,
length: length, // 圆柱体长度
topRadius: 0.0, // 圆柱体顶部半径
bottomRadius: 200000.0, // 圆柱体底部半径
heightReference: Cesium.HeightReference.CLAMP_TO_GROUND,
fill: true,
material: color1.withAlpha(0.5),
outline: true,
outlineColor: color1,
outlineWidth: 1.0,
numberOfVerticalLines: 0, // 沿轮廓的周长绘制的垂直线的数量
shadows: Cesium.ShadowMode.DISABLED,
slices: 23, // 圆柱周围的边缘数量
},
});
for (var ind = 0; ind < 292; ind++) {
var time = Cesium.JulianDate.addSeconds(
viewer.clock.startTime,
3 * ind,
new Cesium.JulianDate()
);
var position = satellite.position.getValue(time);
var cartographic = viewer.scene.globe.ellipsoid.cartesianToCartographic(
position
);
var temp = cartographic;
var lat = Cesium.Math.toDegrees(cartographic.latitude),
lng = Cesium.Math.toDegrees(cartographic.longitude),
hei = cartographic.height / 1;
property.addSample(time, Cesium.Cartesian3.fromDegrees(lng, lat, hei));
}
cylinderEntity1.position = property;
cylinderEntity1.position.setInterpolationOptions({
//设定位置的插值算法
interpolationDegree: 5,
interpolationAlgorithm: Cesium.LagrangePolynomialApproximation,
});
cylinderEntity1.position.setInterpolationOptions({
//设定位置的插值算法
interpolationDegree: 5,
interpolationAlgorithm: Cesium.LagrangePolynomialApproximation,
});
var projectionPositions = [];
viewer.clock.onTick.addEventListener(function() {
var satellitePosition = satellite.position.getValue(
viewer.clock.currentTime
);
if (satellitePosition != null) {
var projectionPosition = viewer.scene.globe.ellipsoid.cartesianToCartographic(
satellitePosition
);
projectionPosition.height = 0;
projectionPositions.push(
viewer.scene.globe.ellipsoid.cartographicToCartesian(projectionPosition)
);
}
});
}
2.2 过境预测
了解卫星轨道的各位朋友都知道卫星轨道是一个宽幅，他是扫过去的，所以他这个宽幅在扫到预测区域的时候，有可能是这个区域在这个幅宽内，有可能在这个幅宽外，也有可能与这个幅宽相交，还有可能是这个幅宽在区域内，因为区域有可能比这个幅宽还要大。那这样，在计算的时候就要把这些情况都要考虑进去。
因为卫星每时每刻都在运动，产生的面有无数个，为了减小计算量，首先需要利用预测区域的四至经纬度进行筛选一下，将区域外的都筛掉，这样会大大减少计算量。在实际项目中，仅这一步操作就筛选掉了三万多次计算，如果这三万多次计算加上，这个计算时间可想而知。
首先需要知道所选区域的四至经纬度，通过四至经纬度确定预测区域的外接矩形，也就是
//外接矩形
var tempArray = [];
for (var i = 0; i < data.features[0].geometry.coordinates.length; i++) {
tempArray.push(getRect(data.features[0].geometry.coordinates[i][0]));
}
let xmax = tempArray.reduce((a, b) => {
return a.xmax > b.xmax ? a.xmax : b.xmax;
});
let xmin = tempArray.reduce((a, b) => {
return a.xmin < b.xmin ? a.xmin : b.xmin;
});
let ymax = tempArray.reduce((a, b) => {
return a.ymax > b.ymax ? a.ymax : b.ymax;
});
let ymin = tempArray.reduce((a, b) => {
return a.ymin < b.ymin ? a.ymin : b.ymin;
});
var temp_polyGeo = [xmin, ymin, xmax, ymax];
然后利用外接矩形先进行一次筛选，最后在将预测区域与卫星宽幅进行裁剪，也就是过境分析。
//过境分析
function transitAnalysis() {
gjarr = [];
sessionStorage.setItem("wxgjlist", gjarr);
indexSatelliteInArea = 1;
htmlStr1a = "";
viewer.dataSources._dataSources.forEach((geoJsonDataSource) => {
if (geoJsonDataSource instanceof Cesium.GeoJsonDataSource) {
viewer.dataSources._dataSources.forEach((czmlDataSource) => {
if (czmlDataSource instanceof Cesium.CzmlDataSource) {
var polygon =
geoJsonDataSource.entities.values[0].polygon.hierarchy._value
.positions;
for (var i = 0; i < czmlDataSource.entities._entities.length; i++) {
var satellite = czmlDataSource.entities._entities._array[i];
SatelliteInArea(satellite, polygon);
}
}
});
}
});
}
function SatelliteInArea(satellite, polygon) {
var timeData = satellite.position._property._times;
var projectionPositions = [];
var color1 = new Cesium.Color(
satellite.path.material.color._value.red,
satellite.path.material.color._value.green,
satellite.path.material.color._value.blue
);
for (var timeinr = 0; true; timeinr++) {
var time = Cesium.JulianDate.addSeconds(
viewer.clock.startTime,
1 * timeinr,
new Cesium.JulianDate()
);
var time1 = Cesium.JulianDate.addSeconds(
viewer.clock.startTime,
1 * (timeinr + 1),
new Cesium.JulianDate()
);
if (time > timeData[timeData.length – 1]) {
break;
} else {
var position = satellite.position.getValue(time);
//下一秒卫星的位置
var position1 = satellite.position.getValue(time1);
if (position == null || position1 == null) {
break;
}
var cartographic = viewer.scene.globe.ellipsoid.cartesianToCartographic(
position
);
var cartographiclst = viewer.scene.globe.ellipsoid.cartesianToCartographic(
position1
);
var distance = mapSatileWidth.get(satellite.label.text._value) / 2;
var lat0 = Cesium.Math.toDegrees(cartographic.latitude);
var lng0 = Cesium.Math.toDegrees(cartographic.longitude);
var latlst = Cesium.Math.toDegrees(cartographiclst.latitude);
var lnglst = Cesium.Math.toDegrees(cartographiclst.longitude);
//当前时刻卫星上下点坐标
var lat1 = getDisLnglat(lng0, lat0, 90, distance).lat;
var lng1 = getDisLnglat(lng0, lat0, 90, distance).lng;
var lat2 = getDisLnglat(lng0, lat0, 270, distance).lat;
var lng2 = getDisLnglat(lng0, lat0, 270, distance).lng;
var line = turf.lineString([[lng2, lat2], [lng1, lat1]]);
if ((lng0 < temp_polyGeo[2] && lng0 > temp_polyGeo[0]) ||
(lng1 < temp_polyGeo[2] && lng1 > temp_polyGeo[0]) ||
(lng2 < temp_polyGeo[2] && lng2 > temp_polyGeo[0]) ||
(temp_polyGeo[0] > lng2 && temp_polyGeo[2] < lng1)) {
if (isPointInPolygon(lng0, lat0, polyGeo) || isPointInPolygon(lng1, lat1, polyGeo) || isPointInPolygon(lng2, lat2, polyGeo) || turf.intersect(line, polyGeo)) {
projectionPositions.push(Cesium.Cartesian3.fromDegrees(lng0, lat0));
} else {
if (projectionPositions.length >= 2) {
var cartographic1 = viewer.scene.globe.ellipsoid.cartesianToCartographic(
projectionPositions[0]
);
var cartographic2 = viewer.scene.globe.ellipsoid.cartesianToCartographic(
projectionPositions[1]
);
if (cartographic1.latitude > cartographic2.latitude) {
var nowTime = time
.toString()
.replace(/T/, " ")
.replace(/Z/, "")
.replace(/\\.\\d+/g, "");
let obj = {
index: indexSatelliteInArea,
time: nowTime,
address: satellite.label.text._value,
};
indexSatelliteInArea = indexSatelliteInArea + 1;
gjarr.push(obj);
}
}
projectionPositions = [];
}
}
//下一时刻卫星上下点坐标
var latlst1 = getDisLnglat(lnglst, latlst, 90, distance).lat;
var lnglst1 = getDisLnglat(lnglst, latlst, 90, distance).lng;
var latlst2 = getDisLnglat(lnglst, latlst, 270, distance).lat;
var lnglst2 = getDisLnglat(lnglst, latlst, 270, distance).lng;
//1假设卫星轨迹过行政区域，不算宽幅
//2假设卫星左侧轨迹过行政区域
//3假设卫星右侧轨迹过行政区域
//4假设行政区域在卫星左侧轨迹内
//5假设行政区域卫星右侧轨迹内
if (((lat0 > temp_polyGeo[3]) && (latlst < temp_polyGeo[3]) && (lng0 > temp_polyGeo[0]) && (lng0 < temp_polyGeo[2])) ||
((lat2 > temp_polyGeo[3]) && (latlst2 < temp_polyGeo[3]) && (lng2 > temp_polyGeo[0]) && (lng2 < temp_polyGeo[2])) ||
((lat1 > temp_polyGeo[3]) && (latlst1 < temp_polyGeo[3]) && (lng1 > temp_polyGeo[0]) && (lng1 < temp_polyGeo[2])) ||
((lat0 > temp_polyGeo[3]) && (latlst < temp_polyGeo[3]) && (lng0 > temp_polyGeo[3]) && (lng2 < temp_polyGeo[0])) ||
((lat0 > temp_polyGeo[3]) && (latlst < temp_polyGeo[3]) && (lng1 > temp_polyGeo[3]) && (lng0 < temp_polyGeo[0]))) {
var coor00 = [lng1, lat1];
var coor03 = [lng2, lat2];
var coor04 = [lng1, lat1];
}
if (((lat0 > temp_polyGeo[1]) && (latlst < temp_polyGeo[1]) && (lng0 > temp_polyGeo[0]) && (lng0 < temp_polyGeo[2])) ||
((lat2 > temp_polyGeo[1]) && (latlst2 < temp_polyGeo[1]) && (lng2 > temp_polyGeo[0]) && (lng2 < temp_polyGeo[2])) ||
((lat1 > temp_polyGeo[1]) && (latlst1 < temp_polyGeo[1]) && (lng1 > temp_polyGeo[0]) && (lng1 < temp_polyGeo[2])) ||
((lat0 > temp_polyGeo[1]) && (latlst < temp_polyGeo[1]) && (lng0 > temp_polyGeo[3]) && (lng2 < temp_polyGeo[0])) ||
((lat0 > temp_polyGeo[1]) && (latlst < temp_polyGeo[1]) && (lng1 > temp_polyGeo[3]) && (lng0 < temp_polyGeo[0]))) {
var coor01 = [lnglst1, latlst1];
var coor02 = [lnglst2, latlst2];
var poly01 = turf.polygon([[
coor00,
coor01,
coor02,
coor03,
coor04
]]);
//卫星过境预测部分外接矩形实例测试
// viewer.dataSources
// .add(
// Cesium.GeoJsonDataSource.load(poly1, {
// stroke: Cesium.Color.YELLOW, //设置多边形轮廓的默认颜色
// fill: Cesium.Color.YELLOW.withAlpha(0.5), //多边形的内部默认颜色
// strokeWidth: 10, //轮廓的宽度
// clamToGround: true, //让地图贴地
// })
// )
if (polyGeo.geometry.coordinates[0].length > 0) {
for (let i = 0; i < polyGeo.geometry.coordinates.length; i++) {
var poly2 = turf.polygon([polyGeo.geometry.coordinates[0][i]]);
var intersection = turf.intersect(poly01, poly2);
if (intersection != null) {
viewer.dataSources
.add(
Cesium.GeoJsonDataSource.load(intersection, {
stroke: color1, //设置多边形轮廓的默认颜色
fill: color1.withAlpha(0.3), //多边形的内部默认颜色
strokeWidth: 10, //轮廓的宽度
clamToGround: true, //让地图贴地
})
)
}
}
}
}
}
}
//将过境分析结果添加到页面列表中
for (var i = 0; i < gjarr.length; i++) {
for (var j = i + 1; j < gjarr.length; j++) {
if (
gjarr[i].address == gjarr[j].address &&
gjarr[i].time.slice(0, 10) == gjarr[j].time.slice(0, 10)
) {
gjarr.splice(j, 1);
j–;
}
}
gjarr[i].index = i + 1;
}
sessionStorage.setItem("wxgjlist", JSON.stringify(gjarr));
}
/**
*
* @param {*} lng 经度 122
* @param {*} lat 纬度 24
* @param {*} brng 方位角 0~360度
* @param {*} dist 90000距离(米)
*
*/
function getDisLnglat(lng, lat, brng, dist) {
var a = 6378137;
var b = 6356752.3142;
var f = 1 / 298.257223563;
var lon1 = lng * 1;
var lat1 = lat * 1;
var s = dist;
var alpha1 = brng * (Math.PI / 180)
var sinAlpha1 = Math.sin(alpha1);
var cosAlpha1 = Math.cos(alpha1);
var tanU1 = (1 – f) * Math.tan(lat1 * (Math.PI / 180));
var cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
var sigma1 = Math.atan2(tanU1, cosAlpha1);
var sinAlpha = cosU1 * sinAlpha1;
var cosSqAlpha = 1 – sinAlpha * sinAlpha;
var uSq = cosSqAlpha * (a * a – b * b) / (b * b);
var A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 – 175 * uSq)));
var B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 – 47 * uSq)));
var sigma = s / (b * A), sigmaP = 2 * Math.PI;
while (Math.abs(sigma – sigmaP) > 1e-12) {
var cos2SigmaM = Math.cos(2 * sigma1 + sigma);
var sinSigma = Math.sin(sigma);
var cosSigma = Math.cos(sigma);
var 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 = s / (b * A) + deltaSigma;
}
var tmp = sinU1 * sinSigma – cosU1 * cosSigma * cosAlpha1;
var lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
(1 – f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
var lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma – sinU1 * sinSigma * cosAlpha1);
var C = f / 16 * cosSqAlpha * (4 + f * (4 – 3 * cosSqAlpha));
var L = lambda – (1 – C) * f * sinAlpha *
(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
var revAz = Math.atan2(sinAlpha, -tmp); // final bearing
var lngLatObj = { lng: lon1 + L * (180 / Math.PI), lat: lat2 * (180 / Math.PI) }
return lngLatObj;
}
function isPointInPolygon(lng, lat, poly) {
var pt = turf.point([lng, lat]);
return turf.booleanPointInPolygon(pt, poly);
}
如果有用请多多关注，谢谢~
