cesium修改源码支持月球坐标系

随着各国探月任务的不断推进，月球地图可视化的需求日益增长。Cesium 作为一款功能强大的三维地球可视化框架，默认只支持地球坐标系（主要是 WGS84），无法直接用于月球数据的展示与分析。本文将详细介绍如何通过修改 Cesium 源码，为其添加月球坐标系的支持，实现月球数据的准确可视化。

月球与地球坐标系的差异

在进行修改前，我们需要理解月球与地球坐标系的主要区别：

1. 形状差异：地球是扁球体，而月球更接近于标准球体
2. 尺寸差异：月球半径约为 1737.4 公里，远小于地球的 6378.137 公里
3. 扁率差异：地球的扁率约为 1/298.257，而月球的扁率可以忽略不计
4. 坐标系统：月球通常使用基于它自身参数的坐标系，如 MOON2000 这些差异导致我们不能直接使用地球的椭球体参数来处理月球数据，否则会产生严重的位置偏差和变形。

const ExternalParameter = {
  ellipsoidCoordinates: null,
  changeApproximateTerrainHeights: false,
  minTerrainHeight: -20000,
  maxTerrainHeight: 20000,
  hasVertexNormals: false, //mpt坡度图
};
ExternalParameter.isMoon = function () {
  return (
    ExternalParameter.ellipsoidCoordinates === "MOON" ||
    ExternalParameter.ellipsoidCoordinates === "MOON2000"
  );
};
ExternalParameter.isVertexNormals = function () {
  return ExternalParameter.hasVertexNormals;
};
export default ExternalParameter;

const MOONRadiiSquared = new Cartesian3(
  1737400.0 * 1737400.0,
  1737400.0 * 1737400.0,
  // eslint-disable-next-line no-loss-of-precision
  1737400.0 * 1737400.0
);

Cartesian3.fromRadians = function (longitude, latitude, height, ellipsoid, result) {
  //>>includeStart('debug', pragmas.debug);
  Check.typeOf.number("longitude", longitude);
  Check.typeOf.number("latitude", latitude);
  //>>includeEnd('debug');

  height = defaultValue(height, 0.0);
  /** Geoworld Add **/
  let radiiSquared = null;
  if (ExternalParameter.ellipsoidCoordinates === "CGCS2000") {
    radiiSquared = CGCS2000RadiiSquared;
  } else if (
    ExternalParameter.ellipsoidCoordinates === "MOON" ||
    ExternalParameter.ellipsoidCoordinates === "MOON2000"
  ) {
    radiiSquared = MOONRadiiSquared;
  } else {
    radiiSquared = !defined(ellipsoid) ? Cartesian3._ellipsoidRadiiSquared : ellipsoid.radiiSquared;
  }
  const cosLatitude = Math.cos(latitude);
  scratchN.x = cosLatitude * Math.cos(longitude);
  scratchN.y = cosLatitude * Math.sin(longitude);
  scratchN.z = Math.sin(latitude);
  scratchN = Cartesian3.normalize(scratchN, scratchN);

  Cartesian3.multiplyComponents(radiiSquared, scratchN, scratchK);
  const gamma = Math.sqrt(Cartesian3.dot(scratchN, scratchK));
  scratchK = Cartesian3.divideByScalar(scratchK, gamma, scratchK);
  scratchN = Cartesian3.multiplyByScalar(scratchN, height, scratchN);

  if (!defined(result)) {
    result = new Cartesian3();
  }
  return Cartesian3.add(scratchK, scratchN, result);
};

const moonOneOverRadii = new Cartesian3(1.0 / 1737400.0, 1.0 / 1737400.0, 1.0 / 1737400.0);
const moonOneOverRadiiSquared = new Cartesian3(
  1.0 / (1737400.0 * 1737400.0),
  1.0 / (1737400.0 * 1737400.0),
  1.0 / (1737400.0 * 1737400.0)
);

Cartographic.fromCartesian = function (cartesian, ellipsoid, result) {
  // const oneOverRadii = defined(ellipsoid)
  //   ? ellipsoid.oneOverRadii
  //   : wgs84OneOverRadii;
  // const oneOverRadiiSquared = defined(ellipsoid)
  //   ? ellipsoid.oneOverRadiiSquared
  //   : wgs84OneOverRadiiSquared;
  // const centerToleranceSquared = defined(ellipsoid)
  //   ? ellipsoid._centerToleranceSquared
  //   : wgs84CenterToleranceSquared;
  //** Geoworld Add **/
  let oneOverRadii = null;
  let oneOverRadiiSquared = null;
  const centerToleranceSquared = defined(ellipsoid)
    ? ellipsoid._centerToleranceSquared
    : Cartographic._ellipsoidCenterToleranceSquared;
  if (ExternalParameter.ellipsoidCoordinates === "CGCS2000") {
    oneOverRadii = CGCS2000OneOverRadii;
    oneOverRadiiSquared = CGCS2000OneOverRadiiSquared;
  } else if (
    ExternalParameter.ellipsoidCoordinates === "MOON" ||
    ExternalParameter.ellipsoidCoordinates === "MOON2000"
  ) {
    oneOverRadii = moonOneOverRadii;
    oneOverRadiiSquared = moonOneOverRadiiSquared;
  } else {
    oneOverRadii = defined(ellipsoid)
      ? ellipsoid.oneOverRadii
      : Cartographic._ellipsoidOneOverRadii;
    oneOverRadiiSquared = defined(ellipsoid)
      ? ellipsoid.oneOverRadiiSquared
      : Cartographic._ellipsoidOneOverRadiiSquared;
  }
  //`cartesian is required.` is thrown from scaleToGeodeticSurface
  const p = scaleToGeodeticSurface(
    cartesian,
    oneOverRadii,
    oneOverRadiiSquared,
    centerToleranceSquared,
    cartesianToCartographicP
  );

  if (!defined(p)) {
    return undefined;
  }

  let n = Cartesian3.multiplyComponents(p, oneOverRadiiSquared, cartesianToCartographicN);
  n = Cartesian3.normalize(n, n);

  const h = Cartesian3.subtract(cartesian, p, cartesianToCartographicH);

  const longitude = Math.atan2(n.y, n.x);
  const latitude = Math.asin(n.z);
  const height = CesiumMath.sign(Cartesian3.dot(h, cartesian)) * Cartesian3.magnitude(h);

  if (!defined(result)) {
    return new Cartographic(longitude, latitude, height);
  }
  result.longitude = longitude;
  result.latitude = latitude;
  result.height = height;
  return result;
};

/** Geoworld Add **/
Ellipsoid.MOON2000 = Object.freeze(
  // eslint-disable-next-line no-loss-of-precision
  new Ellipsoid(1737400.0, 1737400.0, 1737400.0)
);

Ellipsoid.MOON = Object.freeze(
  new Ellipsoid(CesiumMath.LUNAR_RADIUS, CesiumMath.LUNAR_RADIUS, CesiumMath.LUNAR_RADIUS)
);