Cesium加载多坐标系OGC服务

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的属性。

image.png

我们在看下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);
image.png
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,271评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,275评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,151评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,550评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,553评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,559评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,924评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,580评论 0赞 257
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,826评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2赞 320
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,363评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,940评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2赞 342

推荐阅读更多精彩内容