gdal栅格切图

package gdal.cutmap.raster;

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.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;

import java.io.File;
import java.io.FileInputStream;
import java.io.InputStream;
import java.sql.Connection;
import java.sql.PreparedStatement;
import java.util.UUID;

import static org.apache.commons.io.IOUtils.toByteArray;

/**
 * https://www.jianshu.com/p/98389e3168e8
 */
public class DemoRasterTest {

    public static void main(String[] args) throws FactoryException {
        gdal.AllRegister();
        String fileName = "C:\\Users\\h\\Desktop\\tif\\map.tif";
        Dataset dataset = gdal.Open(fileName, gdalconstConstants.GA_ReadOnly);
        double[] ori_transform = dataset.GetGeoTransform();
        int rasterCount = dataset.getRasterCount();

        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.1 首先获取原始影像的地理坐标范围
        int yCount = dataset.getRasterXSize();
        int xCount = dataset.getRasterYSize();
        System.err.println("yCount:"+yCount+"xCount:"+xCount);
        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]);
        CoordinateReferenceSystem crs = CRS.decode("EPSG:4326", true);
        ReferencedEnvelope imageBound = new ReferencedEnvelope(lonMin, lonMax, latMin, latMax, crs);
//4.2 获取原始影像的像素分辨率
        // 原始图像东西方向像素分辨率
        double src_w_e_pixel_resolution = (lonMax - lonMin) / xCount;
        // 原始图像南北方向像素分辨率
        double src_n_s_pixel_resolution = (latMax - latMin) / yCount;
// 4.3 根据原始影像地理范围求解切片行列号  // 经纬度转瓦片编号

        int zoom = 10;
        int tileRowMax = lat2tile(latMin, zoom); // 纬度  -90 ——90  lat
        int tileRowMin = lat2tile(latMax, zoom);  // 经度 -180 -- 180 lon
        int tileColMin = lon2tile(lonMin, zoom);
        int tileColMax = lon2tile(lonMax, zoom);
//4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集
        int row=44;
        int col=34;
        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";
        System.out.println("pngPath="+pngPath);
        // 使用PNG驱动将内存中的图像数组写入文件
        Driver pngDriver = gdal.GetDriverByName("PNG");
        Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS);

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

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

//        try {
//            File tileFile = new File(pngPath);
//            InputStream fis = new FileInputStream(tileFile);
//            byte[] imgBytes = toByteArray(fis);
//
//            String uuid = UUID.randomUUID().toString();
//            // 获取数据库链接
//            Connection conn=null;
//            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();
//        }
    }
    // 经纬度转瓦片编号
    public static int  lon2tile(double lon, int zoom) {
        return (int)(Math.floor((lon+180)/360*Math.pow(2,zoom)));
    }
    public static int  lat2tile(double lat, int zoom)  {
        return (int)(Math.floor((1-Math.log(Math.tan(lat*Math.PI/180) + 1/Math.cos(lat*Math.PI/180))/Math.PI)/2 *Math.pow(2,zoom)));
    }

    // 瓦片编号转经纬度
    private static double tile2lon(int col, int zoom) {
        return col / Math.pow(2.0, zoom) * 360.0 - 180;
    };

    private static double tile2lat(int row, int zoom) {
        double n = Math.PI - (2.0 * Math.PI * row) / Math.pow(2.0, zoom);
        return Math.toDegrees(Math.atan(Math.sinh(n)));
    };
}