基于GDAL的影像切图工具

男娘i 2022-11-07 01:53 221阅读 0赞

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/][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方法对影像指定范围进行读取与压缩。

推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。

Band in_band1 = dataset.GetRasterBand(1);
    Band in_band2 = dataset.GetRasterBand(2);
    Band in_band3 = dataset.GetRasterBand(3);
    
    int[] band1BuffData = new int[256 * 256 * gdalconst.GDT_Int32];
    int[] band2BuffData = new int[256 * 256 * gdalconst.GDT_Int32];
    int[] band3BuffData = new int[256 * 256 * gdalconst.GDT_Int32];
    
    in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf,gdalconst.GDT_Int32, band1BuffData, 0, 0);
    in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf,gdalconst.GDT_Int32, band2BuffData, 0, 0);
    in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, image_Xbuf, image_Ybuf,gdalconst.GDT_Int32, band3BuffData, 0, 0);

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();
    }

[http_www.gisinternals.com]: http://www.gisinternals.com/