上篇文章介绍了地图瓦片加载计算的相关的原理,学习原理性的东西,到底有没有用,怎么用? 日常工作中一般有成熟的第三方工具或者库可以实现。虽然很多时候是有工具的,大部分工具需要手动操作不能做到批量自化、动化操作,或者工具能力有限不能满足需求。作为技术同学不能只局限于使用工具,要学会创造工具,提供最优的解决方案。 本系列教程教程语言实现均为 JS语言,,读者需要具有 JS 开发基础并熟悉 Node。
本文主要介绍如何实现地图瓦片切片工具,地图切片除了需要了解瓦片计算逻辑,还会涉及到数据投影、分辨率,重采样相关能力这里我们使用第三方 库 Gdal 的node 版本 node-gdal。工具功能介绍可以看这篇文章
需求
遥感数据一般比较大,几百MB 或者几个 GB,如何在前端浏览器可视化确实是个问题,数据量不仅数据加载缓慢、而且地图渲染性能有很高的要求。瓦片是解决数据量大最好的方案,那问题怎么把一个 GeoTiff切成瓦片,目前没有合适的开源工具。
切片流程
读取文件
Gdal node 版本
naturalatlas.github.io/node-gdal/c…
读取文件为GeoTiff 文件,这里使用 Gdal 进行读取。
// 引入node-gdal 库
var gdal = require("gdal-next")
var dataset = gdal.open("../landcover_JiangXi.tif");
读取到文件 Gdal 提供了一些方法获取栅格影像的相关信息如坐标系、大小、分辨率,数据等等
重采样生成影像金字塔
为什么需要重采样,瓦片切片需要根据地图缩放层级、生成不同分辨率的影像也就生成瓦片金字塔模型。 不同分辨率
分辨率计算
地球周长除以地图像素大小既为分辨率,地球半径为固定值,地图像素大小跟地图缩放等级相关
- 初始分辨率第 0 层级 一张 256 * 256 的瓦片即可包含整个世界。
var initialResolution =2 * Math.PI * 6378137 / 256
- 各层级分辨率计算
var resolution = initialResolution / Math.pow(2, N)
重采样大小
- width、height 原栅格数据长、宽
- xSize、ySize 原栅格数据分辨率
- resolution 目标分辨率
const x = Math.round(width * xSize / resolution),
const y = Math.abs(Math.round(height * ySize / resolution))
相关参数计算完成之后,我们就可以进行重采样了,这里可以使用 reprojectImage 方法。
reprojectRaster(dataset, zoom, path) {
const driver = dataset.driver;
const tileHelp = this.tileHelper;
const resolution = tileHelp.resolution(zoom);
const rasterSize = dataset.rasterSize;
const newSize = this.scaleZoomSize(rasterSize.x, rasterSize.y, dataset.geoTransform[1], dataset.geoTransform[5], resolution);
const outDs = driver.create(path, newSize.x, newSize.y, 1, this.getDataType());
const gcp = dataset.getGCPs();
const gcpProject = dataset.srs.toWKT();
const [x, xtr, xr, yx, yr, ytr] = dataset.geoTransform;
outDs.geoTransform = [x, resolution, xr, yx, yr, -resolution];
outDs.setGCPs(gcp, gcpProject);
outDs.bands.get(1).colorInterpretation = dataset.bands.get(1).colorInterpretation;
const option = {
src: dataset,
dst: outDs,
s_srs: dataset.srs,
t_srs: dataset.srs,
resampling: gdal.GRA_NearestNeighbor,
};
gdal.reprojectImage(option)
outDs.flush();
outDs.close();
}
瓦片切片-平面分块
重采样之后,我们得到了不同分辨率的栅格影像,接下来就可以根据规则进行切片了。
栅格数据范围
栅格数据的经纬度范围直接 Gdal 可以读取,范围主要用于计算数据包含的瓦片范围进而进行切片。原始数据坐标系一般不是是经纬度坐标系,需要进行坐标转换。 原始坐标范围计算,ds.geoTransform 具体参数如下。 geoTransform坐标转换
getRasterExtent() {
const ds = this.dataset
const size = ds.rasterSize
const minCorner = { x: 0, y: size.y };
const maxCorner = { x: size.x, y: 0 };
const geotransform = ds.geoTransform;
const wgs84 = gdal.SpatialReference.fromEPSG(4326);
const coord_transform = new gdal.CoordinateTransformation(ds.srs, wgs84);
const min = {
x: geotransform[0] + minCorner.x * geotransform[1] + minCorner.y * geotransform[2],
y: geotransform[3] + minCorner.x * geotransform[4] + minCorner.y * geotransform[5]
}
const max = {
x: geotransform[0] + maxCorner.x * geotransform[1] + maxCorner.y * geotransform[2],
y: geotransform[3] + maxCorner.x * geotransform[4] + maxCorner.y * geotransform[5]
}
const minWgs84 = coord_transform.transformPoint(min);
const maxWgs84 = coord_transform.transformPoint(max);
return [minWgs84.y, minWgs84.x, maxWgs84.y, maxWgs84.x]
}
瓦片范围计算
瓦片范围计算就需要了解瓦片技术规则了,经纬度坐标如何转瓦片坐标,如何计算参照上篇文章
boundsToTileExtent(minLon, minLat, maxLon, maxLat, zoom) {
const [minTx, minTy] = this.lonLatToTile(minLon, maxLat, zoom);
const [maxTx, maxTy] = this.lonLatToTile(maxLon, minLat, zoom);
return [[minTx, minTy], [maxTx, maxTy]]
}
读取数据
瓦片行列号范围计算好之后我们就可以根据行列号逐个读取瓦片的数据了。读取数据需要处理为边缘的瓦片,数据超出区域的情况。
// 读取指定起点的瓦片数据
readRaster(origin) {
const tileHelp = this.tileHelper;
const band = this.dataset.bands.get(1);
const rasterSize = this.dataset.rasterSize;
const start = [Math.min(Math.max(0, origin[0]), rasterSize.x - 1), Math.min(Math.max(0, origin[1]), rasterSize.y - 1)];
const end = [Math.min(origin[0] + 256, rasterSize.x), Math.min(origin[1] + 256, rasterSize.y)]
const x = origin[0] < 0 ? - origin[0] : 0;
const y = origin[1] < 0 ? - origin[1] : 0;
const height= end[1] -start[1]
const width= end[0] -start[0]
const data = Array.from(band.pixels.read(start[0], start[1], width, height))
return {
data,
x,
y,
height,
width
}
}
写入数据
写入数据也需要注意,在边缘的瓦片数据,数据并不是完整的,需要计算数据写入起始坐标。
writeTiffTile(x, y, z, data) {
const dataset = this.dataset;
const driver = dataset.driver;
const tileHelp = this.tileHelper;
if (!fs.existsSync(`${this.outPath}/${z}`)) {
fs.mkdirSync(`${this.outPath}/${z}`)
}
if (!fs.existsSync(`${this.outPath}/${z}/${x}`)) {
fs.mkdirSync(`${this.outPath}/${z}/${x}`)
}
const outDs = driver.create(`${this.outPath}/${z}/${x}/${y}.tiff`, tileHelp.tileSize, tileHelp.tileSize, 1, this.getDataType());
const gcp = dataset.getGCPs();
const gcpProject = dataset.srs.toWKT();
const resolution = tileHelp.resolution(z);
const originMeters = tileHelp.pixelsToMeters(x * 256, y * 256, z)
const geoTransform = [originMeters[0], resolution, 0, originMeters[1], 0, -resolution];
outDs.geoTransform = geoTransform;
outDs.setGCPs(gcp, gcpProject)
const outBand = outDs.bands.get(1);
outBand.colorInterpretation = dataset.bands.get(1).colorInterpretation;
outBand.noDataValue = this.noDataValue;
outDs.srs = dataset.srs;
outBand.fill(this.noDataValue);
// TODO 数据类型
outBand.pixels.write(data.x, data.y, data.width, data.height, new DataTypeEnum[this.getDataType()](data.data))
outBand.computeStatistics(true)
outDs.flush()
outDs.close();
}
数据可视化
完成上面的布局,我们就可以实现栅格数据切片了,不再需要担心数据量大,无法显示了。
附录
- 在线demo 查看 l7.antv.vision/zh/examples…
- 切片工具源码:github.com/lzxue/raste…
- 可视化引擎L7: github.com/antvis/l7