基于GDAL的影像切图工具

1 功能需求

对一副3857坐标系的卫星影像按照谷歌TMS进行切片,并将切片文件存储在符合MBTiles规范的SQLite数据库中

2 依赖库

GDAL ,GeoTools

3 环境搭建

由于个人对Java语言比较熟悉，因此使用Java语言开发环境，Java工程复用原有的SpringBoot工程，然后在pom文件中添加对GeoTools库的依赖即可实现GeoTools的引入。

3.1 Windows环境下搭建GDAL Java Binding开发环境

GDAL是当前最流行的遥感影像处理库，提供了一系列对遥感影像数据的操作函数，支持常见的影像数据格式，并提供了相应的抽象结构方便进行算法操作。GDAL的源代码是用C和C++编写的，并通过SWIG发布其他语言版本的bindings。
首先访问http://www.gisinternals.com/下载对应版本的编译好的GDAL库，下载并解压到本地文件夹下，注意路径不能有中文。然后在计算机-环境变量Path中添加对GDAL DLL的引用。E:\soft\GIS\gdal\bin;E:\soft\GIS\gdal\bin\gdal\java
在E:\soft\GIS\gdal\bin\gdal\java即可找到GDAL对应的JAR包，在工程中引用后即完成了GDAL的引入。完成后运行测试代码，成功显示影像信息后即代表GDAL库引入成功。

public class GDALTest {
    public static void main(String[] args) {
        gdal.AllRegister();
        String fileName = "D:\\images\\2011_output_cai.tif";
        Dataset dataset = gdal.Open(fileName, gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
            System.err.println(gdal.GetLastErrorMsg());

            System.exit(1);
        }
        double[] ori_transform = dataset.GetGeoTransform();
        System.out.println(String.format("Origin = (%s, %s)", ori_transform[0], ori_transform[3]));
        System.out.println(String.format("Pixel Size = (%s, %s)", ori_transform[1], ori_transform[5]));
    }
}

4 切片处理流程

4.1 首先获取原始影像的地理坐标范围

double[] ori_transform = dataset.GetGeoTransform();
double latMin = ori_transform[3];
double latMax = ori_transform[3] + (yCount * ori_transform[1]);
double lonMin = ori_transform[0];
double lonMax = ori_transform[0] + (xCount * ori_transform[1]);
ReferencedEnvelope imageBound = new ReferencedEnvelope(lonMin, lonMax, latMin, latMax, crs);

4.2 获取原始影像的像素分辨率

int xCount = dataset.getRasterXSize();
int yCount = dataset.getRasterYSize();
// 原始图像东西方向像素分辨率
double src_w_e_pixel_resolution = (lonMax - lonMin) / xCount;
// 原始图像南北方向像素分辨率
double src_n_s_pixel_resolution = (latMax - latMin) / yCount;

4.3 根据原始影像地理范围求解切片行列号

int tileRowMax = Utils.getOSMTileYFromLatitude(latMin, zoom);
int tileRowMin = Utils.getOSMTileYFromLatitude(latMax, zoom);
int tileColMin = Utils.getOSMTileXFromLongitude(lonMin, zoom);
int tileColMax = Utils.getOSMTileXFromLongitude(lonMax, zoom);

4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集

double tempLatMin = tile2lat(row + 1, zoom);
double tempLatMax = tile2lat(row, zoom);
double tempLonMin = tile2lon(col, zoom);
double tempLonMax = tile2lon(col + 1, zoom);
ReferencedEnvelope tileBound = new ReferencedEnvelope(tempLonMin, tempLonMax, tempLatMin,tempLatMax, crs);
ReferencedEnvelope intersect = tileBound.intersection(imageBound);

4.5 求解当前切片的像素分辨率(默认切片大小为256*256)

// 切片东西方向像素分辨率
double dst_w_e_pixel_resolution = (tempLonMax - tempLonMin) / 256;
// 切片南北方向像素分辨率
double dst_n_s_pixel_resolution = (tempLatMax - tempLatMin) / 256;

4.6 计算交集的像素信息

// 求切图范围和原始图像交集的起始点像素坐标
int offset_x = (int) ((intersect.getMinX() - lonMin) / src_w_e_pixel_resolution);
int offset_y = (int) Math.abs((intersect.getMaxY() - latMax) / src_n_s_pixel_resolution);

// 求在切图地理范围内的原始图像的像素大小
int block_xsize = (int) ((intersect.getMaxX() - intersect.getMinX()) / src_w_e_pixel_resolution);
int block_ysize = (int) ((intersect.getMaxY() - intersect.getMinY()) / src_n_s_pixel_resolution);

// 求原始图像在切片内的像素大小
int image_Xbuf = (int) Math.ceil((intersect.getMaxX() - intersect.getMinX()) / dst_w_e_pixel_resolution);
int image_Ybuf = (int) Math.ceil(Math.abs((intersect.getMaxY() - intersect.getMinY()) / dst_n_s_pixel_resolution));

// 求原始图像在切片中的偏移坐标
int imageOffsetX = (int) ((intersect.getMinX() - tempLonMin) / dst_w_e_pixel_resolution);
int imageOffsetY = (int) Math.abs((intersect.getMaxY() - tempLatMax) / dst_n_s_pixel_resolution);
imageOffsetX = imageOffsetX > 0 ? imageOffsetX : 0;
imageOffsetY = imageOffsetY > 0 ? imageOffsetY : 0;

4.7 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。

4.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, image_Xbuf, image_Ybuf, band1BuffData);
dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData);
dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData);
alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData);

String tileFolder = System.getProperty("java.io.tmpdir");
String pngPath = tileFolder + File.separator + zoom + "c" + col + "r" + row +".png";
// 使用PNG驱动将内存中的图像数组写入文件
Driver pngDriver = gdal.GetDriverByName("PNG");
Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS);

// 释放内存
msmDS.FlushCache();
pngDs.delete();

4.8 读取临时文件，并转换成二进制存储到SQLite

try {
    InputStream fis = new FileInputStream(tileFile);
    byte[] imgBytes = toByteArray(fis);

    String uuid = UUID.randomUUID().toString();
    PreparedStatement mapPs = conn.prepareStatement("insert into map(zoom_level,tile_column,tile_row,tile_id) values (?,?,?,?)");
    PreparedStatement imagesPs = conn.prepareStatement("insert into images(tile_data,tile_id) values (?,?)");

         imagesPs.setBytes(1, imgBytes);
         imagesPs.setString(2, uuid);

         mapPs.setInt(1, zoom);
         mapPs.setInt(2, col);
         mapPs.setInt(3, row);
         mapPs.setString(4, uuid);

         imagesPs.executeUpdate();
         mapPs.executeUpdate();

         imagesPs.close();
         mapPs.close();
         fis.close();
         tileFile.delete();
} catch (Exception e) {
    e.printStackTrace();
}