gdal栅格切图

今天药忘吃喽~ 2022-10-05 06:50 173阅读 0赞

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