Cesium載入ArcGIS Server4490且orgin -400 400的切片服務

孤草之魂發表於2023-04-25
Server

Cesium在使用載入Cesium.ArcGisMapServerImageryProvider載入切片服務時，預設只支援wgs84的4326座標系，不支援CGCS2000的4490座標系。

如果是ArcGIS釋出的4490座標系的切片服務，如果原點在orgin X: -180.0Y: 90.0的情況下，我們可以透過WebMapTileServiceImageryProvider按照WMTS的方式載入（需符合OGC標準的WMTS型別）。

但是對於ArcGIS釋出4490座標系的切片服務，如果原點在orgin X: -400.0Y: 400.0的情況下，我們無法實現載入，本文透過示例演示實現Cesium載入ArcGIS Server4490且orgin -400 400的切片服務。

本文使用：

Cesium原始碼版本：1.94

原始碼打包測試參考另一篇文件：cesium原始碼編譯除錯及呼叫全過程

另外，本文的一些解釋需要對切片原理有一定了解（想進一步瞭解的話可以去看看相關文件說明，不想了解的話按步驟修改就行了）。

為了能夠除錯原始碼，打包的時候使用命令：npm run combine

一、透過修改原始碼實現ArcGIS的切片服務，需要修改的原始碼檔案包括：

  • ArcGisMapServerImageryProvider
  • GeographicTilingScheme
  • Ellipsoid

1、修改ArcGisMapServerImageryProvider類

透過檢視ArcGisMapServerImageryProvider（\Source\Scene\ArcGisMapServerImageryProvider.js）原始碼，我們發現它不支援CGCS2000的4490座標系（僅支援wgs84的4326座標系）：

找到metadataSuccess方法，進行以下修改：

（1）讀取切片後設資料時增加支援wkid 4490座標系的判斷，同時將切片資訊也傳入，目的是為了後面在獲取行列號xy時，可以透過讀取切片資訊，使用自定義方法改寫行列號的獲取方式。

else if (data.tileInfo.spatialReference.wkid === 4490) {
        that._tilingScheme = new GeographicTilingScheme({
          ellipsoid: options.ellipsoid,
          tileInfo: data.tileInfo,
          rectangle: that._rectangle,
          numberOfLevelZeroTilesX: options.numberOfLevelZeroTilesX,
          numberOfLevelZeroTilesY: options.numberOfLevelZeroTilesY
        });
        
        that._tilingScheme._tileInfo = data.tileInfo;//附加自定義屬性
      }

 具體位置如圖所示：

（2）fullExtent範圍增加wkid 4490座標系判斷。

else if (data.fullExtent.spatialReference.wkid === 4326 || data.fullExtent.spatialReference.wkid === 4490) {
            that._rectangle = Rectangle.fromDegrees(
              data.fullExtent.xmin,
              data.fullExtent.ymin,
              data.fullExtent.xmax,
              data.fullExtent.ymax
            );
          }

 程式碼位置：

 2、修改GeographicTilingScheme類

GeographicTilingScheme類的位置是：\Source\Core\GeographicTilingScheme.js

透過增加4490座標系的橢球、矩陣範圍等定義，4490座標系預設橢球為CGCS2000，矩陣範圍為（-180，-90，180，90），開放矩陣範圍的目的就是為了支援自定義的origin原點。

if (defined(options.tileInfo)
    && defined(options.tileInfo.spatialReference)
    && defined(options.tileInfo.spatialReference.wkid)
    && options.tileInfo.spatialReference.wkid == 4490) {
    this._tileInfo = options.tileInfo;
    this._ellipsoid = defaultValue(options.ellipsoid, Ellipsoid.CGCS2000);
    this._rectangle = defaultValue(options.rectangle, Rectangle.fromDegrees(-180, -90, 180, 90));
    this._numberOfLevelZeroTilesX = defaultValue(options.numberOfLevelZeroTilesX, 4);
    this._numberOfLevelZeroTilesY = defaultValue(options.numberOfLevelZeroTilesY, 2);
  }
  else {
    this._ellipsoid = defaultValue(options.ellipsoid, Ellipsoid.WGS84);
    this._rectangle = defaultValue(options.rectangle, Rectangle.MAX_VALUE);
    this._numberOfLevelZeroTilesX = defaultValue(options.numberOfLevelZeroTilesX, 2);
    this._numberOfLevelZeroTilesY = defaultValue(options.numberOfLevelZeroTilesY, 1);
  }
  this._projection = new GeographicProjection(this._ellipsoid);

程式碼位置： 

 （2）修改切片矩陣計算獲取行列號數量xy值的原型方法getNumberOfXTilesAtLevel和getNumberOfYTilesAtLeve

/**
 * Gets the total number of tiles in the X direction at a specified level-of-detail.
 *
 * @param {Number} level The level-of-detail.
 * @returns {Number} The number of tiles in the X direction at the given level.
 */
GeographicTilingScheme.prototype.getNumberOfXTilesAtLevel = function (level) {
  // return this._numberOfLevelZeroTilesX << level;
  if (!defined(this._tileInfo)) {
    return this._numberOfLevelZeroTilesX << level
  } else { // 使用切片矩陣計算
    var currentMatrix = this._tileInfo.lods.filter(function (item) {
      return item.level === level
    })
    var currentResolution = currentMatrix[0].resolution
    // return Math.round(360 / (this._tileInfo.rows * currentResolution))
    return Math.round(CesiumMath.toDegrees(CesiumMath.TWO_PI * 2) / (this._tileInfo.rows * currentResolution));
  }
};

/**
 * Gets the total number of tiles in the Y direction at a specified level-of-detail.
 *
 * @param {Number} level The level-of-detail.
 * @returns {Number} The number of tiles in the Y direction at the given level.
 */
GeographicTilingScheme.prototype.getNumberOfYTilesAtLevel = function (level) {
  // return this._numberOfLevelZeroTilesY << level;
  if (!defined(this._tileInfo)) {
    return this._numberOfLevelZeroTilesY << level
  } else { // 使用切片矩陣計算
    var currentMatrix = this._tileInfo.lods.filter(function (item) {
      return item.level === level
    })
    var currentResolution = currentMatrix[0].resolution
    // return Math.round(180 / (this._tileInfo.cols * currentResolution))
    return Math.round(CesiumMath.toDegrees(CesiumMath.TWO_PI * 2) / (this._tileInfo.cols * currentResolution));
  }
};

程式碼位置：

這段程式碼和參考文章的程式碼存在一定出入，在文末會做詳細說明。

3、修改Ellipsoid類，定義2000橢球引數

Ellipsoid類位置：\Source\Core\Ion.js

 定義2000橢球引數：

/**
 * An Ellipsoid instance initialized to the CGCS2000 standard.
 *
 * @type {Ellipsoid}
 * @constant
 */
Ellipsoid.CGCS2000 = Object.freeze(
  new Ellipsoid(6378137.0, 6378137.0, 6356752.31414035585)
);

 程式碼位置：

二、程式碼呼叫

原始碼修改後，為了能夠除錯原始碼，使用npm run combine打包下程式碼，並在呼叫地方修改下引用，

測試用html頁面（做測試的頁面，有一些其他程式碼，可以忽略，留意本文需要的程式碼部分）：

<!DOCTYPE html>
<html lang="en">
<head>
    <!-- Use correct character set. -->
    <meta charset="utf-8"/>
    <!-- Tell IE to use the latest, best version. -->
    <meta http-equiv="X-UA-Compatible" content="IE=edge"/>
    <!-- Make the application on mobile take up the full browser screen and disable user scaling. -->
    <meta
            name="viewport"
            content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no"
    />
    <title>cesium載入影像和向量資料</title>
    <script src="../Build/CesiumUnminified/Cesium.js"></script>
    <style>
        @import url(../Build/CesiumUnminified/Widgets/widgets.css);

        html,
        body,
        #cesiumContainer {
            width: 100%;
            height: 100%;
            margin: 0;
            padding: 0;
            overflow: hidden;
        }
    </style>
</head>
<body>
<div id="cesiumContainer"></div>
<script>
    //天地圖token
    let TDT_tk = "b0df1f950b1fd6914abe9e17079c0345";
    //Cesium token
    let cesium_tk = "eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJqdGkiOiI1NThjYTk0MC03YjQwLTQ3YWYtOTY5Yy04NDk3OTJmMmI4NDciLCJpZCI6NzAxOTMsImlhdCI6MTYzNDA4ODE1N30.TB1v9XXATQLUGE5GNki_fYFMHddIyQ9arXPIx65e09s";
    //天地圖影像
    let TDT_IMG_C = "http://{s}.tianditu.gov.cn/img_c/wmts?service=wmts&request=GetTile&version=1.0.0" +
        "&LAYER=img&tileMatrixSet=c&TileMatrix={TileMatrix}&TileRow={TileRow}&TileCol={TileCol}" +
        "&style=default&format=tiles&tk=" + TDT_tk;

    //標註
    let TDT_CIA_C = "http://{s}.tianditu.gov.cn/cia_c/wmts?service=wmts&request=GetTile&version=1.0.0" +
        "&LAYER=cia&tileMatrixSet=c&TileMatrix={TileMatrix}&TileRow={TileRow}&TileCol={TileCol}" +
        "&style=default&format=tiles&tk=" + TDT_tk;
    
    //初始頁面載入
    //Cesium.Ion.defaultAccessToken = cesium_tk;
    var cgs2000Ellipsolid = Cesium.Ellipsoid.CGCS2000
    var cgs2000GeographicProj = new Cesium.GeographicProjection(cgs2000Ellipsolid)
    let viewer = new Cesium.Viewer('cesiumContainer', {
        // baseLayerPicker: false,
        timeline: true,
        homeButton: true,
        fullscreenButton: true,
        infoBox: true,
        animation: true,
        shouldAnimate: true,
        mapProjection: cgs2000GeographicProj
        //imageryProvider: layer, //設定預設底圖
    });
    let rightTilt = true;
    if (rightTilt) {
      viewer.scene.screenSpaceCameraController.tiltEventTypes = [
        Cesium.CameraEventType.RIGHT_DRAG,
        Cesium.CameraEventType.PINCH,
        {
          eventType: Cesium.CameraEventType.LEFT_DRAG,
          modifier: Cesium.KeyboardEventModifier.CTRL
        },
        {
          eventType: Cesium.CameraEventType.RIGHT_DRAG,
          modifier: Cesium.KeyboardEventModifier.CTRL
        }
      ]
      viewer.scene.screenSpaceCameraController.zoomEventTypes = [
        Cesium.CameraEventType.MIDDLE_DRAG,
        Cesium.CameraEventType.WHEEL,
        Cesium.CameraEventType.PINCH
      ]
    }
    var handler = new Cesium.ScreenSpaceEventHandler(viewer.scene.canvas);
    handler.setInputAction(function(evt) {
        var cartesian=viewer.camera.pickEllipsoid(evt.position,viewer.scene.globe.ellipsoid);
        var cartographic=Cesium.Cartographic.fromCartesian(cartesian);
        var lng=Cesium.Math.toDegrees(cartographic.longitude);//經度值
        var lat=Cesium.Math.toDegrees(cartographic.latitude);//緯度值
        var mapPosition={x:lng,y:lat,z:cartographic.height};//cartographic.height的值始終為零。
        alert("longitude:" + lng + ";latitude:" + lat );
    }, Cesium.ScreenSpaceEventType.LEFT_CLICK);
    
    viewer.imageryLayers.remove(viewer.imageryLayers.get(0))
    //新增tms
    let tms = {};
    tms.url =  "http://10.0.7.16:81/tms";
    if (tms) {
      const layerInfo = {
        url: tms.url,
        fileExtension: tms.fileExtension || 'jpg',
        maximumLevel: tms.maxZoom || 7,
        name: 'tms'
      }
      const tmsService = new Cesium.TileMapServiceImageryProvider(layerInfo)
      tmsService.layerInfo = layerInfo
    }
    //新增地形
    let terrain = {};
    terrain.url =  "http://data.marsgis.cn/terrain";
    if (terrain) {
        const terrainLayer = new Cesium.CesiumTerrainProvider({
            url: terrain.url
            })
        viewer.terrainProvider = terrainLayer
    }

    _matrixIds = ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18"]
    //呼叫影響中文註記服務
    /*viewer.imageryLayers.addImageryProvider(new Cesium.WebMapTileServiceImageryProvider({
        url: TDT_CIA_C,
        layer: "tdtImg_c",
        style: "default",
        format: "tiles",
        tileMatrixSetID: "c",
        subdomains: ["t0", "t1", "t2", "t3", "t4", "t5", "t6", "t7"],
        tilingScheme: new Cesium.GeographicTilingScheme(),
        tileMatrixLabels: ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"],
        maximumLevel: 50,
        show: false
    }))*/
    
    /*var myGeographicTilingScheme = new Cesium.GeographicTilingScheme({
        ellipsoid: cgs2000Ellipsolid,
        rectangle: Cesium.Rectangle.fromDegrees(-400, -399.9999999999998, 400, 399.9999999999998),
        numberOfLevelZeroTilesX: 4,
        numberOfLevelZeroTilesY: 4
    })*/
    var world4490 = new Cesium.ArcGisMapServerImageryProvider({
        url: 'http://10.1.88.200:6080/arcgis/rest/services/test/global4490ori400/MapServer',
        //tilingScheme: myGeographicTilingScheme,
        rectangle: Cesium.Rectangle.fromDegrees(-400, -320, 320, 400),
        ellipsoid:cgs2000Ellipsolid,
        numberOfLevelZeroTilesX: 4,
        numberOfLevelZeroTilesY: 4
    });
    viewer.imageryLayers.addImageryProvider(world4490);
    //viewer.imageryLayers.addImageryProvider(world);
    //使用ArcGisMapServerImageryProvider載入影像沒成功，改用WebMapServiceImageryProvider
    //var world = new Cesium.ArcGisMapServerImageryProvider({
        //url:'http://10.1.88.200:6080/arcgis/rest/services/test/globaltdt5/MapServer',
    //});
    //viewer.imageryLayers.addImageryProvider(world);
    var arcgisyx = new Cesium.WebMapServiceImageryProvider({
    url:'http://10.1.88.200:6080/arcgis/rest/services/test/globaltdt5/MapServer/tile/{z}/{y}/{x}',
    layers:[0]
    });
    // viewer.imageryLayers.addImageryProvider(arcgisyx);
    var china = new Cesium.ArcGisMapServerImageryProvider({
        url:'http://10.1.88.200:6080/arcgis/rest/services/test/china4490/MapServer'
    });
    viewer.imageryLayers.addImageryProvider(china);
</script>
</body>
</html>

1、定義橢球體部分：

 var cgs2000Ellipsolid = Cesium.Ellipsoid.CGCS2000
    var cgs2000GeographicProj = new Cesium.GeographicProjection(cgs2000Ellipsolid)

2、初始化Cesium.Viewer時，增加2000的mapProjection

let viewer = new Cesium.Viewer('cesiumContainer', {
        // baseLayerPicker: false,
        timeline: true,
        homeButton: true,
        fullscreenButton: true,
        infoBox: true,
        animation: true,
        shouldAnimate: true,
        mapProjection: cgs2000GeographicProj
        //imageryProvider: layer, //設定預設底圖
    });

3、呼叫關鍵程式碼，載入圖層

var world4490 = new Cesium.ArcGisMapServerImageryProvider({
        url: 'http://10.1.88.200:6080/arcgis/rest/services/test/global4490ori400/MapServer',
        //tilingScheme: myGeographicTilingScheme,
        rectangle: Cesium.Rectangle.fromDegrees(-400, -320, 320, 400),
        ellipsoid:cgs2000Ellipsolid,
        numberOfLevelZeroTilesX: 4,
        numberOfLevelZeroTilesY: 4
    });
    viewer.imageryLayers.addImageryProvider(world4490);

說明：

（1）服務說明，釋出了一個全球影像4490座標系-400,400起點的資料（切片方案保證其他引數與-180,90一致，可從本文文末獲取切片方案xml檔案用來切片測試）

 （2）在新建ArcGisMapServerImageryProvider時，可以不設定tilingScheme。發現ArcGisMapServerImageryProvider裡有新建tilingScheme，如果設定了tilingScheme，發現這裡的行列數引數沒有起作用，故直接在新建ArcGisMapServerImageryProvider傳入行列數引數，並在類中新建切片方案的時候讀取。

rectangle切片範圍：發現用(-400, -400, 400, 400)帶入整個地圖偏移了：

對切片的原理進一步瞭解後：

針對-400,400起點切片，為了保證和-180,90按照0級兩列一行的大小，如上圖，按照格網一樣90°大小，劃分成4行4列，切片範圍應該是(-400, -320, 320, 400),請求的切片應該6塊，行列是（行在前列在後）：（1,1），（1,2），（1,3），（2,1），（2,2），（2,3）

按照這樣設定，此時再次執行後，能夠正常載入-400,400起點切片了：

以上是展開的效果，球體的效果：

 

注：

本文參考文章：https://blog.csdn.net/wokao253615105/article/details/123462643

關於getNumberOfXTilesAtLevel和getNumberOfYTilesAtLeve這段程式碼和上述參考的文件存在一定出入做進一步說明。

但是因為直接複製後發現不能正常載入，透過切片原理判斷得出，計算的行列式數量不對，改成了2π*2：

 

這裡0級的話，-400,400起點切片獲取的xy切片數量應該是4行4列，透過反推應該是4π，如果是-180,90起點的話，是2行1列，則x應該用2π，y應該用π，原文章應該是針對-180,90起點的計算方式。

這段程式碼的寫法只支援-400,400起點，因為只是為了測試能夠把-400,400起點的切片資料，偷懶直接這麼寫了。如果要同時支援兩種起點，這裡的寫法應該要改成公式：（右頂點X-左頂點X）/(256*解析度)；（上頂點Y-下頂點Y）/(256*解析度)

如：

-180,90起點：[180-（-180）]/(256*解析度)      範圍右頂點經度-左頂點經度，256是因為切片大小是256*256,當0級時候，解析度為0.7031250000026057

-400,400起點：[320-（-400）]/(256*解析度)      範圍右頂點經度-左頂點經度，256是因為切片大小是256*256,當0級時候，解析度為0.7031250000026057

 

本文使用的資料說明及測試下載：

影像資料：用來切4490起點-400,400的切片資料，全球影像，天地圖5級資料合成的。

連結：https://pan.baidu.com/s/14SrGonqHG9gL6ixixIFVUw

向量資料：用來驗證與影像資料疊加是否大致吻合，全國行政區資料，包括省會，國界、省界

連結：https://pan.baidu.com/s/1FdEaoD8rs5ogLZ2Pwa3Xxw

切片方案：起點-400,400的切片方案

連結：https://pan.baidu.com/s/1vTunkbMMaV3wFyXF1YwbGA

測試頁面：cesiumlayer.html

連結：https://pan.baidu.com/s/1uNHBg-dY2xIgyyTuzllaPw

<本文完>

 

相關文章