基于GDAL的Java遥感影像瓦片切片生成方式
感谢https://www.jianshu.com/p/98389e3168e8 提供了基础代码实现方式，但由于作者只提供了代码片段，导致某些实现不完整，这里对其补充了下。同时我们常用web地图引擎都是基于Web墨卡托切片，所以在对方实现上增加了投影转换。
package com.itfsw.example;
import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconst;
import org.gdal.gdalconst.gdalconstConstants;
import org.gdal.osr.SpatialReference;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.locationtech.jts.geom.Coordinate;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import java.io.File;
import java.util.UUID;
public class Tif2Tiles {
/**
* 瓦片大小
*/
public final static int TILE_SIZE = 256;
/**
* web 墨卡托 EPSG 编码
*/
public final static int WEB_MERCATOR_EPSG_CODE = 3857;
/**
* Wgs84 EPSG 编码
*/
public final static int WGS_84_EGPS_CODE = 4326;
/**
* Tif 切片
*
* @param tifFile
* @param outputDir
* @param minZoom
* @param maxZoom
*/
public static void tif2Tiles(String tifFile, String outputDir, int minZoom, int maxZoom) throws FactoryException, TransformException {
Dataset srcDs = gdal.Open(tifFile, gdalconstConstants.GA_ReadOnly);
if (srcDs == null) {
System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
System.err.println(gdal.GetLastErrorMsg());
return;
}
// 1. 影像重投影到web墨卡托
String tmpFile = System.getProperty("java.io.tmpdir") + File.separator + UUID.randomUUID().toString();
SpatialReference spatialReference = new SpatialReference();
spatialReference.ImportFromEPSG(WEB_MERCATOR_EPSG_CODE);
Dataset webMercatorDs = gdal.AutoCreateWarpedVRT(srcDs, null, spatialReference.ExportToWkt(), gdalconst.GRA_Bilinear);
Dataset dataset = gdal.GetDriverByName("GTiff").CreateCopy(tmpFile, webMercatorDs);
srcDs.delete();
webMercatorDs.delete();
try {
// 2. 获取遥感影像经纬度范围，计算遥感影像像素分辨率
// 获取原始影像的地理坐标范围
double[] geoTransform = dataset.GetGeoTransform();
// 获取原始影像的像素分辨率
int xSize = dataset.getRasterXSize();
int ySize = dataset.getRasterYSize();
// 计算经纬度范围
double lngMin = geoTransform[0];
double latMax = geoTransform[3];
double lngMax = lngMin + (xSize * geoTransform[1]) + (ySize * geoTransform[2]);
double latMin = latMax + (xSize * geoTransform[4]) + (ySize * geoTransform[5]);
// EPSG:3857 坐标转Wgs84经纬度
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:" + WEB_MERCATOR_EPSG_CODE);
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:" + WGS_84_EGPS_CODE);
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, false);
Coordinate lngLatMax = JTS.transform(new Coordinate(lngMax, latMax), null, transform);
Coordinate lngLatMin = JTS.transform(new Coordinate(lngMin, latMin), null, transform);
lngMax = lngLatMax.getY();
latMax = lngLatMax.getX();
lngMin = lngLatMin.getY();
latMin = lngLatMin.getX();
// 原始图像东西方向像素分辨率
double srcWePixelResolution = (lngMax - lngMin) / xSize;
// 原始图像南北方向像素分辨率
double srcNsPixelResolution = (latMax - latMin) / ySize;
for (int zoom = minZoom; zoom <= maxZoom; zoom++) {
// 3. 根据原始影像地理范围求解切片行列号
int[] tilePointMax = coordinates2tile(lngMax, latMax, zoom);
int[] tilePointMin = coordinates2tile(lngMin, latMin, zoom);
int tileRowMax = tilePointMin[1];
int tileColMax = tilePointMax[0];
int tileRowMin = tilePointMax[1];
int tileColMin = tilePointMin[0];
for (int row = tileRowMin; row <= tileRowMax; row++) {
for (int col = tileColMin; col <= tileColMax; col++) {
// 4. 求原始影像地理范围与指定缩放级别指定行列号的切片交集
double tempLatMin = tile2lat(row + 1, zoom);
double tempLatMax = tile2lat(row, zoom);
double tempLonMin = tile2lng(col, zoom);
double tempLonMax = tile2lng(col + 1, zoom);
ReferencedEnvelope tileBound = new ReferencedEnvelope(tempLonMin, tempLonMax, tempLatMin, tempLatMax, DefaultGeographicCRS.WGS84);
ReferencedEnvelope imageBound = new ReferencedEnvelope(lngMin, lngMax, latMin, latMax, DefaultGeographicCRS.WGS84);
ReferencedEnvelope intersect = tileBound.intersection(imageBound);
// 5. 求解当前切片的像素分辨率(默认切片大小为256*256)
// 切片东西方向像素分辨率
double dstWePixelResolution = (tempLonMax - tempLonMin) / 256;
// 切片南北方向像素分辨率
double dstNsPixelResolution = (tempLatMax - tempLatMin) / 256;
// 6. 计算交集的像素信息
// 求切图范围和原始图像交集的起始点像素坐标
int offsetX = (int) ((intersect.getMinX() - lngMin) / srcWePixelResolution);
int offsetY = (int) Math.abs((intersect.getMaxY() - latMax) / srcNsPixelResolution);
// 求在切图地理范围内的原始图像的像素大小
int blockXSize = (int) ((intersect.getMaxX() - intersect.getMinX()) / srcWePixelResolution);
int blockYSize = (int) ((intersect.getMaxY() - intersect.getMinY()) / srcNsPixelResolution);
// 求原始图像在切片内的像素大小
int imageXBuf = (int) Math.ceil((intersect.getMaxX() - intersect.getMinX()) / dstWePixelResolution);
int imageYBuf = (int) Math.ceil(Math.abs((intersect.getMaxY() - intersect.getMinY()) / dstNsPixelResolution));
// 求原始图像在切片中的偏移坐标
int imageOffsetX = (int) ((intersect.getMinX() - tempLonMin) / dstWePixelResolution);
int imageOffsetY = (int) Math.abs((intersect.getMaxY() - tempLatMax) / dstNsPixelResolution);
imageOffsetX = imageOffsetX > 0 ? imageOffsetX : 0;
imageOffsetY = imageOffsetY > 0 ? imageOffsetY : 0;
// 7. 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
// 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
Band inBand1 = dataset.GetRasterBand(1);
Band inBand2 = dataset.GetRasterBand(2);
Band inBand3 = dataset.GetRasterBand(3);
int[] band1BuffData = new int[TILE_SIZE * TILE_SIZE * gdalconst.GDT_Int32];
int[] band2BuffData = new int[TILE_SIZE * TILE_SIZE * gdalconst.GDT_Int32];
int[] band3BuffData = new int[TILE_SIZE * TILE_SIZE * gdalconst.GDT_Int32];
inBand1.ReadRaster(offsetX, offsetY, blockXSize, blockYSize, imageXBuf, imageYBuf, gdalconst.GDT_Int32, band1BuffData, 0, 0);
inBand2.ReadRaster(offsetX, offsetY, blockXSize, blockYSize, imageXBuf, imageYBuf, gdalconst.GDT_Int32, band2BuffData, 0, 0);
inBand3.ReadRaster(offsetX, offsetY, blockXSize, blockYSize, imageXBuf, imageYBuf, gdalconst.GDT_Int32, band3BuffData, 0, 0);
// 8. 将切片数据写入文件
// 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
Driver memDriver = gdal.GetDriverByName("MEM");
Dataset msmDS = memDriver.Create("msmDS", 256, 256, 4);
Band dstBand1 = msmDS.GetRasterBand(1);
Band dstBand2 = msmDS.GetRasterBand(2);
Band dstBand3 = msmDS.GetRasterBand(3);
// 设置alpha波段数据,实现背景透明
Band alphaBand = msmDS.GetRasterBand(4);
int[] alphaData = new int[256 * 256 * gdalconst.GDT_Int32];
for (int index = 0; index < alphaData.length; index++) {
if (band1BuffData[index] > 0) {
alphaData[index] = 255;
}
}
// 写各个波段数据
dstBand1.WriteRaster(imageOffsetX, imageOffsetY, imageXBuf, imageYBuf, band1BuffData);
dstBand2.WriteRaster(imageOffsetX, imageOffsetY, imageXBuf, imageYBuf, band2BuffData);
dstBand3.WriteRaster(imageOffsetX, imageOffsetY, imageXBuf, imageYBuf, band3BuffData);
alphaBand.WriteRaster(imageOffsetX, imageOffsetY, imageXBuf, imageYBuf, alphaData);
String pngPath = outputDir + File.separator + zoom + "-" + col + "-" + row + ".png";
// 使用PNG驱动将内存中的图像数组写入文件
Driver pngDriver = gdal.GetDriverByName("PNG");
Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS);
// 释放内存
msmDS.FlushCache();
msmDS.delete();
pngDs.delete();
}
}
}
} finally {
// 释放并删除临时文件
dataset.delete();
new File(tmpFile).delete();
}
}
/**
* 经纬度坐标转瓦片坐标
*
* @param lng
* @param lat
* @param zoom
* @return
*/
public static int[] coordinates2tile(double lng, double lat, int zoom) {
int x = (int) Math.floor((lng + 180) / 360 * (1 << zoom));
int y = (int) Math.floor((1 - Math.log(Math.tan(Math.toRadians(lat)) + 1 / Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1 << zoom));
if (x < 0) {
x = 0;
}
if (x >= (1 << zoom)) {
x = ((1 << zoom) - 1);
}
if (y < 0) {
y = 0;
}
if (y >= (1 << zoom)) {
y = ((1 << zoom) - 1);
}
return new int[]{x, y};
}
/**
* 瓦片坐标转经度
*
* @param x
* @param z
* @return
*/
public static double tile2lng(double x, int z) {
return x / Math.pow(2.0, z) * 360.0 - 180;
}
/**
* 瓦片坐标转纬度
*
* @param y
* @param z
* @return
*/
public static double tile2lat(double y, int z) {
double n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2.0, z);
return Math.toDegrees(Math.atan(Math.sinh(n)));
}
/**
* @param args
* @throws FactoryException
* @throws TransformException
*/
public static void main(String[] args) throws FactoryException, TransformException {
gdal.AllRegister();
tif2Tiles("D:/test.tif", "D:/tiles", 3, 15);
}
}