Cesium使用DEM地形数据(.tif)创建自定义Primitive

最新推荐文章于 2024-10-09 09:25:56 发布

微笑面对bug

最新推荐文章于 2024-10-09 09:25:56 发布

阅读量1.9k

点赞数 1

分类专栏： WebGL实战文章标签：前端

本文链接：https://blog.csdn.net/qq_25169689/article/details/126509479

版权

WebGL实战专栏收录该内容

8 篇文章 10 订阅

订阅专栏

DEM数据源：地理空间数据云 GDEMV3 30M 分辨率数字高程数据

第一步：用python把tif文件解析成高程数据，生成一个js文件，里面是高程数组，长度3601*3601

from osgeo import gdal, osr
import numpy as np

fname = 'ASTGTMV003_N30E104_dem';
url = './'+fname+'.tif'
ds = gdal.Open(url)
gt = ds.GetGeoTransform()
print(gt)

dem = ds.ReadAsArray()
print(dem.max())

f = open(fname+'.js', "a")

f.write('var dataDem = [')
for item in dem:
    for item2 in item:
        f.write(str(item2)+',')

f.write(']')
f.close()

第二步：把地形数据分成小块创建geomotry(一个3601*3601顶点的模型浏览器会崩溃，所以要切成小块)。

贴图使用mapbox的web墨卡托影像图，uv计算涉及到web墨卡托与经纬度的相互转换。

影像图的级别的选择：按照当前地形块的四个角在同一张影像图、且级别最大为原则。

源码地址：Apps/dem_tif_2.html · changjiuxiong/Cesium-1.62Test - Gitee.com

完整代码：

<!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>Hello World!</title>
  <script src="../Build/Cesium/Cesium.js"></script>
  <script src="./file/test.js"></script>


  <style>
    @import url(../Build/Cesium/Widgets/widgets.css);
    html, body, #cesiumContainer {
      width: 100%; height: 100%; margin: 0; padding: 0; overflow: hidden;
    }
  </style>
</head>
<body>
<div id="cesiumContainer">
</div>
<script>
    Cesium.Ion.defaultAccessToken='eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJqdGkiOiJmNjJjMzY0OS1hZGQxLTRiZmYtYWYwNS03NmIyM2MwMDgwZDAiLCJpZCI6MTIzMTgsInNjb3BlcyI6WyJhc3IiLCJnYyJdLCJpYXQiOjE1NjA4NDQ3Mjd9.OLTL_rs2gAi2R9zoztBHcJPDHnVl2Q7OZxRtZhoCeZE';
    let viewer = new Cesium.Viewer("cesiumContainer");

    let data = dataDem;
    let dataMinL = 86;
    let dataMinB = 27;
    let dataMaxL = 87;
    let dataMaxB = 28;

    let maxLength = 20037508.3427892;
    let maxLength2 = maxLength*2;

    let destination = Cesium.Cartesian3.fromDegrees(dataMinL, dataMaxB, 10000);

    let cutWidth = 36;
    let cutWidthL = cutWidth/3600;
    let cutWidthB = cutWidth/3600;
    for(let startY=0; startY<cutWidth*30; startY+=cutWidth){
      for(let startX=0; startX<cutWidth*30; startX+=cutWidth){

        let posList = [];
        let indexList = [];
        let uvList = [];
        let nList = [];

        let ld = [dataMinL+startX/3600, dataMaxB-(startY+cutWidth)/3600];
        let lu = [dataMinL+startX/3600, dataMaxB-startY/3600];
        let rd = [dataMinL+(startX+cutWidth)/3600, dataMaxB-(startY+cutWidth)/3600];
        let ru = [dataMinL+(startX+cutWidth)/3600, dataMaxB-startY/3600];
        let findXY;
        let z1 = 18;
        for(;z1>=0;z1--){
          let ldxy = getXYbyBLZ(ld[0],ld[1],z1);
          let luxy = getXYbyBLZ(lu[0],lu[1],z1);
          let rdxy = getXYbyBLZ(rd[0],rd[1],z1);
          let ruxy = getXYbyBLZ(ru[0],ru[1],z1);

          if(ldxy[0]===luxy[0]&&ldxy[0]===rdxy[0]&&ldxy[0]===ruxy[0] && ldxy[1]===luxy[1]&&ldxy[1]===rdxy[1]&&ldxy[1]===ruxy[1]){
            findXY = ldxy;
            break;
          }
        }
        let url1 = 'https://api.mapbox.com/v4/mapbox.satellite/'+z1+'/'+findXY[0]+'/'+findXY[1]+'.webp?sku=101iGH34Pnrym&access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4M29iazA2Z2gycXA4N2pmbDZmangifQ.-g_vE53SD2WrJ6tFX7QHmA';
        // console.log(url1);

        let count1 = Math.pow(2,z1);
        let ceil1 = maxLength2/count1;

        let minMx = findXY[0]*ceil1 - maxLength;
        let maxMx = minMx + ceil1;

        let maxMy = maxLength - (findXY[1]*ceil1);
        let minMy = maxMy - ceil1;

        let minLB = webMercator2LonLat(minMx, minMy);
        let maxLB = webMercator2LonLat(maxMx, maxMy);

        let lLen = maxLB[0] - minLB[0];
        let bLen = maxLB[1] - minLB[1];

        for(let y=startY; y<=startY+cutWidth; y++){
          for(let x=startX; x<=startX+cutWidth; x++){
            let height = data[y*3601+x];
            let B = dataMaxB - y/3600;
            let L = dataMinL + x/3600;

            let s = (L - minLB[0])/lLen;
            let t = (B - minLB[1])/bLen;
            uvList.push(s,t);

            let v3 = Cesium.Cartesian3.fromDegrees(L, B, height);
            posList.push(v3.x, v3.y, v3.z);

            nList.push(v3.x, v3.y, v3.z);
          }
        }

        let count = cutWidth + 1;
        for (let i = 0; i < count*count - count; i++) {
          if (i % count === count - 1) {
            continue;
          }

          indexList.push(i + 1);
          indexList.push(i);
          indexList.push(i + count);

          indexList.push(i + 1 + count);
          indexList.push(i + 1);
          indexList.push(i + count);
        }

        let pos = new Float64Array(posList);
        let normal = new Float32Array(nList);
        let index = new Uint32Array(indexList);
        let st = new Float32Array(uvList);
        let boundingSphere = Cesium.BoundingSphere.fromVertices(pos);

        //创建geometry
        let geometry = new Cesium.Geometry({
          attributes: {
            position: new Cesium.GeometryAttribute({
              componentDatatype: Cesium.ComponentDatatype.DOUBLE,
              componentsPerAttribute: 3,
              values: pos
            }),

            st: new Cesium.GeometryAttribute({
              componentDatatype: Cesium.ComponentDatatype.FLOAT,
              componentsPerAttribute: 2,
              values: st
            }),

            normal: new Cesium.GeometryAttribute({
              componentDatatype: Cesium.ComponentDatatype.FLOAT,
              componentsPerAttribute: 3,
              values: normal
            }),

          },
          indices: index,
          primitiveType: Cesium.PrimitiveType.TRIANGLES,
          boundingSphere: boundingSphere
        });

        //绘制三角带
        let myInstance = new Cesium.GeometryInstance({
          geometry: geometry,
          // geometry: Cesium.GeometryPipeline.toWireframe(geometry),
        });

        let primitive = new Cesium.Primitive({
          geometryInstances: myInstance,
          appearance: new Cesium.MaterialAppearance({
            material: Cesium.Material.fromType('Image', {
              image: url1,
            }),
            flat: false,
            faceForward: false,
          }),
        });

        viewer.scene.primitives.add(primitive);

      }
    }

    viewer.camera.flyTo({
      destination: destination,
      duration:1,
    });

    //经纬度转Web墨卡托
    function lonLat2WebMercator(l, b) {
      let x = l * 20037508.34/180;
      let y = Math.log(Math.tan((90+b)*Math.PI/360))/(Math.PI/180);
      y = y *20037508.34/180;
      return [x,y];
    }
    //Web墨卡托转经纬度
    function webMercator2LonLat(x, y) {
      let l = x/20037508.34*180;
      let b = y/20037508.34*180;
      b = 180/Math.PI*(2*Math.atan(Math.exp(b*Math.PI/180))-Math.PI/2);
      return [l,b];
    }

    function getXYbyBLZ(l,b,z0){

      let count0 = Math.pow(2,z0);
      let BL = [l, b];
      let ceil = maxLength2/count0;

      let me = lonLat2WebMercator(BL[0],BL[1]);
      let mex = me[0];
      let mey = me[1];

      let x = Math.floor((mex+maxLength)/ceil);
      let y = Math.floor((maxLength2-(mey+maxLength))/ceil);

      return [x,y];
    }

</script>
</body>
</html>