当前位置： 首页>编程语言>正文

Cesium加载多坐标系OGC服务

编程语言2024-10-17 22:13:07

Cesium对坐标系支持情况

Cesium可以加载WMTS、WMS、ArcGIS已经谷歌、必应、天地图等在线服务，但是这些服务都是Web墨卡托或者WGS坐标系，并且Cesium也只支持WGS84（EPSG：43226）、Web墨卡托（EPSG：3857）坐标系的WMTS、ArcGIS MapServer服务，当你加载其他服务的时候，Cesium会提示An error occurred in "WebMapTileServiceImageryProvider": Failed to obtain image tile X: 0 Y: 0 Level: 1.或者An error occurred in "St": Tile spatial reference WKID 4526 is not supported.

原理

WMTS服务

Cesium提供了统一的接口ImageryProvider来加载多种地图服务，对于WMTS服务，Cesium提供了WebMapTileServiceImageryProvider，支持WGS84（EPSG：43226）、Web墨卡托（EPSG：3857）坐标系。同时，WebMapTileServiceImageryProvider构造函数中还提供了tilingScheme: TilingScheme字段，允许我们自定义WMTS瓦片计算方式。

WMTS和TilingScheme

首先我们来了解下WMTS瓦片计算方式，然后看下如何自定义一个自己的TilingScheme对象，来扩充Cesium的多坐标系支持。OGC官网OpenGIS Web Map Tile Service Implementation Standard - Open Geospatial Consortium (ogc.org)讲的很清楚，最重要的就是下面这张图，图中标出来的属性就是WMTS的WMTS的Capabilities.xml的属性。

Cesium加载多坐标系OGC服务,第1张

我们在看下WMTS瓦片计算公式，其实OGC提供的这个规范也给出来了：

pixelSpan = scaleDenominator × 0.28 10^-3 / metersPerUnit(crs)
tileSpanX = tileWidth × pixelSpan
tileSpanY = tileHeight × pixelSpan
tileMatrixMaxX = tileMatrixMinX + tileSpanX × matrixWidth
tileMatrixMinY = tileMatrixMaxY - tileSpanY × matrixHeight

scaleDenominator大家尤其要注意下这个字段，它代表了实际尺寸和像素尺寸的比例，计算的时候不能忽略。
我们来结合一份实际的WMTS服务，来看看具体这些公式是什么意思。

<ows:BoundingBox crs="urn:ogc:def:crs:EPSG::4526">
<ows:LowerCorner>2705867.630193442 3.841100167045193E7</ows:LowerCorner>
<ows:UpperCorner>2776804.1782085816 3.849748247676485E7</ows:UpperCorner>
</ows:BoundingBox>

<TileMatrixSet>
<ows:Title>TileMatrix using 0.28mm</ows:Title>
<ows:Abstract>The tile matrix set that has scale values calculated based on the dpi defined by OGC specification (dpi assumes 0.28mm as the physical distance of a pixel).</ows:Abstract>
<ows:Identifier>default028mm</ows:Identifier>
<ows:SupportedCRS>urn:ogc:def:crs:EPSG::4526</ows:SupportedCRS>
<TileMatrix>
<ows:Identifier>0</ows:Identifier>
<ScaleDenominator>944942.3660750897</ScaleDenominator>
<TopLeftCorner>1.00021E7 3.28768E7</TopLeftCorner>
<TileWidth>256</TileWidth>
<TileHeight>256</TileHeight>
<MatrixWidth>83</MatrixWidth>
<MatrixHeight>108</MatrixHeight>
</TileMatrix>
</TileMatrixSet>
  1. 这是我用ArcGIS Server发布的4526坐标系的MWS部分内容，第一段是说明这个是4526坐标系，注意，在<ows:LowerCorner>2705867.630193442 3.841100167045193E7</ows:LowerCorner>这里，用的是高斯坐标系,也就是x是向北的，y是向东的，和我们的笛卡尔坐标系刚好反过来。
  2. 其中TopLeftCorner就是在你WMTS切片坐标系下的左上角顶点坐标，注意WMTS的坐标系，对于瓦片的tx ty来说，原点在左上角，向右是tx，向下是ty。
  3. 当我想要根据瓦片的tx ty来计算实际坐标x y，我只用根据公式算出tileSpanY tileSpanX，然后加上这个TopLeftCorner就是真实坐标了。同理，当我想要根据实际坐标x y来计算瓦片的tx ty时，首先和TopLeftCorner作差，然后除以tileSpanY tileSpanX，在向下取整，就得出瓦片的坐标了。

ArcGIS MapServer服务

对于ArcGis 的MapServer来说，有两种办法，第一种参考Cesium的ArcGisMapServerImageryProvider写一个自己的Provider，ArcGIS的瓦片计算规则和WMTS的类似，第二种就是用WMTS的方法来加，因为ArcGIS Server也可以发布WMTS服务。

代码

还以刚才那个4526坐标系的WMTS服务为例，我们看看如何写一个自己的TilingScheme。首先来看下TilingScheme的定义。首先是字段：

  1. ellipsoid : Ellipsoid，这个好说，WGS84或者CGCS2000
  2. projection : MapProjection 这个需要自己写，就是在WMTS坐标系和WGS84坐标系下的转换，有两个成员函数project unproject。
  3. rectangle : Rectangle，需要提供在WGS84坐标系下的切片范围，在WMTS的xml文件里可以找到。
    然后是成员函数：
  4. getNumberOfXTilesAtLevel(level) → Number，根据zoom级别返回x方向上的瓦片数量，WMTS的xm文件里也可以获取到
  5. getNumberOfYTilesAtLevel(level) → Number，根据zoom级别返回y方向上的瓦片数量，WMTS的xm文件里也可以获取到
  6. rectangleToNativeRectangle(rectangle, result) → Rectangle，将WMTS指定坐标系下的Rectangle转到WGS84坐标系下
  7. tileXYToNativeRectangle(x, y, level, result) → Rectangle，根据瓦片的坐标和层级，返回在WMTS指定坐标下的瓦片范围
  8. tileXYToRectangle(x, y, level, result) → Rectangle，根据瓦片的坐标和层级，返回在WGS84坐标下的瓦片范围

MapProjection

首先我们先来实现MapProjection，需要借助Proj4js这个库。

// 定义4526坐标系，可以在epsg.io上查
proj4.defs("EPSG:4526", "+proj=tmerc +lat_0=0 +lon_0=114 +k=1 +x_0=38500000 +y_0=0 +ellps=GRS80 +units=m +no_defs +type=crs");

// 定义自己的MapProjection
class MapProjection {
    static transform = proj4("EPSG:4326", "EPSG:4526");

    constructor() {
        this.ellipsoid = Cesium.Ellipsoid.WGS84;
    }

    project(cartographic, result) {
        let lon = Cesium.Math.toDegrees(cartographic.longitude);
        let lat = Cesium.Math.toDegrees(cartographic.latitude);
        // 这里要注意，EPSG规定是纬度在前，经度在后，加上true就表示强制按经纬度的顺序来
        const coord = MapProjection.transform.forward([lon, lat], true);

        if (!Cesium.defined(result)) {
            result = new Cesium.Cartesian3();
        }
        result.x = coord[0];
        result.y = coord[1];
        result.z = cartographic.height;
        return result;
    }

    unproject(cartesian, result) {
        const coord = MapProjection.transform.inverse([cartesian.x, cartesian.y]);
        const lon = Cesium.Math.toRadians(coord[0]);
        const lat = Cesium.Math.toRadians(coord[1]);
        if (!Cesium.defined(result)) {
            result = new Cesium.Cartographic();
        }

        console.log(coord);
        result.longitude = lon;
        result.latitude = lat;
        // result.longitude = coord[0];
        // result.latitude = coord[1];
        result.height = cartesian.z;
        return result;
    }
}

TilingScheme

class ArcGISTilingScheme {
    static scales = [944942.3660750897, 472471.18303754483, 236235.59151877242, 118117.79575938621, 60476.31142880573];
    static xtiles = [83, 166, 332, 664, 1297];
    static ytiles = [108, 216, 431, 862, 1684];
    constructor() {
        this.TopCorner = 1.00021E7;
        this.LeftCorner = 3.28768E7;
        this.ellipsoid = Cesium.Ellipsoid.WGS84;
 
        this.projection = new MapProjection();
        // 这个信息可以在wmts服务那里找到
        this.rectangle = new Cesium.Rectangle(
            Cesium.Math.toRadians(113.11773979064854),
            Cesium.Math.toRadians(24.454099600244135),
            Cesium.Math.toRadians(113.97516978692987),
            Cesium.Math.toRadians(25.09704275817674));
    }

    getNumberOfXTilesAtLevel(level) {
        return ArcGISTilingScheme.xtiles[level];
    }

    getNumberOfYTilesAtLevel(level) {
        return ArcGISTilingScheme.ytiles[level];
    }

    positionToTileXY(position, level, result) {
        let cart = this.projection.project(position);

        const tileY = Math.floor((this.TopCorner - cart.y) / this.pixelSpan(level) / 256);
        const tileX = Math.floor((cart.x - this.LeftCorner) / this.pixelSpan(level) / 256);
        if (!Cesium.defined(result)) {
            result = new Cesium.Cartesian2();
        }

        result.x = tileX;
        result.y = tileY;
        return result;
    }

    rectangleToNativeRectangle(rectangle, result) {
        let coord1 = new Cesium.Cartographic(rectangle.west, rectangle.south, 0);
        let coord2 = new Cesium.Cartographic(rectangle.east, rectangle.north, 0);

        coord1 = this.projection.project(coord1);
        coord2 = this.projection.project(coord2);

        if (!Cesium.defined(result)) {
            result = new Cesium.Rectangle();
        }

        result.west = coord1.x;
        result.south = coord1.y;
        result.east = coord2.x;
        result.north = coord2.y;
        return result;
    }

    tileXYToNativeRectangle(x, y, level, result) {
        let xmin = this.pixelSpan(level) * 256 * x + this.LeftCorner;
        let xmax = this.pixelSpan(level) * (x + 1) * 256 + this.LeftCorner;

        let ymin = this.TopCorner - (y + 1) * 256 * this.pixelSpan(level);
        let ymax = this.TopCorner - (y) * 256 * this.pixelSpan(level);

        if (!Cesium.defined(result)) {
            result = new Cesium.Rectangle();
        }

        result.west = xmin;
        result.south = ymin;
        result.east = xmax;
        result.north = ymax;
        return result;
    }

    tileXYToRectangle(x, y, level, result) {
        const rectangle = this.tileXYToNativeRectangle(x, y, level);
        let coord1 = new Cesium.Cartesian3(rectangle.west, rectangle.south, 0);
        let coord2 = new Cesium.Cartesian3(rectangle.east, rectangle.north, 0);

        coord1 = this.projection.unproject(coord1);
        coord2 = this.projection.unproject(coord2);

        if (!Cesium.defined(result)) {
            result = new Cesium.Rectangle();
        }

        result.west = coord1.longitude;
        result.south = coord1.latitude;
        result.east = coord2.longitude;
        result.north = coord2.latitude;
        return result;
    }

    pixelSpan(level) {
        return ArcGISTilingScheme.scales[level] * 0.28 / 1000.0;
    }
}

应用

应用上就非常简单了，只用加上tilingScheme就好了。

let wmts = new Cesium.WebMapTileServiceImageryProvider({
    url: 'http://192.168.72.193:6080/arcgis/rest/services/sg/MapServer/WMTS',
    layer: 'sg',
    format: "image/png",
    style: 'default',
    tileMatrixSetID: "default028mm",
    maximumLevel: 4,
    tilingScheme: new ArcGISTilingScheme(),
    rectangle: new Cesium.Rectangle(
        Cesium.Math.toRadians(113.11773979064854),
        Cesium.Math.toRadians(24.454099600244135),
        Cesium.Math.toRadians(113.97516978692987),
        Cesium.Math.toRadians(25.09704275817674))
});
viewer.imageryLayers.addImageryProvider(wmts);
Cesium加载多坐标系OGC服务,第2张

https://www.xamrdz.com/lan/5xg1994396.html

相关文章：