1

13693261870

2023-09-20 124486dbffb46b6a80008b0ba591c254643de0ad

 DataLoader/CS/GdalHelper.cs

@@ -1,14 +1,11 @@
using OSGeo.GDAL;
using DataLoader.Model;
using OSGeo.GDAL;
using OSGeo.OGR;
using OSGeo.OSR;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Threading.Tasks;
using System.Windows;
using System.Windows.Markup;

namespace DataLoader.CS
{
@@ -23,10 +20,17 @@

        public static SpatialReference sr4326 = null;

        public static SpatialReference sr4490 = null;

        public static SpatialReference sr104903 = null;

        public static string MOON200 = "GCS_Moon_2000";

        public static List<string> EPSGS = new List<string>() { "EPSG:4326", "EPSG:4490", "ESRI:104903" };

        /// <summary>
        /// 构造函数
        /// </summary>
        private GdalHelper()
        {
            lock (obj)
@@ -37,6 +41,9 @@

                    sr4326 = new SpatialReference(null);
                    sr4326.ImportFromEPSG(4326);

                    sr4490 = new SpatialReference(null);
                    sr4490.ImportFromEPSG(4490);

                    string wkt = "GEOGCS[\"GCS_Moon_2000\",\r\n" +
                        "    DATUM[\"D_Moon_2000\",\r\n" +
@@ -56,6 +63,9 @@
            }
        }

        /// <summary>
        /// 实例
        /// </summary>
        public static GdalHelper Instance
        {
            get
@@ -69,12 +79,15 @@
            }
        }

        /// <summary>
        /// 注册GDAL
        /// </summary>
        public void RegisterGDal()
        {
            string gdalData = Path.Combine(StaticData.BasePath, "gdal-data");
            string gdalData = Path.Combine(Tools.BaseDir, "gdal-data");
            Environment.SetEnvironmentVariable("GDAL_DATA", gdalData);

            string proj7 = Path.Combine(StaticData.BasePath, "proj7\\share");
            string proj7 = Path.Combine(Tools.BaseDir, "proj7\\share");
            Environment.SetEnvironmentVariable("PROJ_LIB", proj7);

            Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // NO,YES
@@ -85,454 +98,131 @@
        }
        #endregion

        public void ReadTiff()
        #region 读取栅格信息
        /// <summary>
        /// 读取栅格信息
        /// </summary>
        public void ReadRasterInfo(ViewData vd)
        {
            string outPath = "D:\\xyz\\ce\\xyz";
            string tif = "D:\\Moon\\data\\dom_tif\\moon.tif";
            Dataset ds = null;
            try
            {
                ds = Gdal.Open(vd.FilePath, Access.GA_ReadOnly);
                if (null == ds || ds.RasterCount == 0 || null == ds.GetSpatialRef()) return;

            Dataset ds = Gdal.Open(tif, Access.GA_ReadOnly);
            int x = ds.RasterXSize;
            int y = ds.RasterYSize;
                vd.Meta.gridsize = string.Format("{0},{1}", ds.RasterXSize, ds.RasterYSize); // 行列数

            SpatialReference sr = ds.GetSpatialRef();
            string srName = sr.GetName(); // GCS_Moon_2000
            // string code = sr4326.GetAuthorityCode(null); // 4326
            if (MOON200 != srName) return;
                SpatialReference sr = ds.GetSpatialRef();
                vd.Meta.coor_sys = sr.GetName(); // 坐标系统
                if (MOON200 == vd.Meta.coor_sys)
                {
                    vd.Meta.epsg = "ESRI:104903"; // EPSG编码
                }
                else
                {
                    string code = sr.GetAuthorityCode(null); // PROJCS、GEOGCS、GEOGCS|UNIT 或 NULL
                    vd.Meta.epsg = string.IsNullOrEmpty(code) ? null : "EPSG:" + code; // EPSG编码
                }

            double[] lt = ImageToGeoSpace(ds, 0, 0); // -179.99998351102391, 89.999983511056882
            double[] rb = ImageToGeoSpace(ds, x, y); // 179.999788852709, -89.9999026708093
                vd.Meta.bands = ds.RasterCount.ToString(); // 波段数
                vd.Meta.band_type = Enum.GetName(typeof(DataType), ds.GetRasterBand(1).DataType); // 数据类型
                ColorTable ct = ds.GetRasterBand(1).GetRasterColorTable();
                vd.Meta.ct = null == ct ? null : ct.ToString(); // 数据颜色表
                vd.Meta.h_datum = null; // 高程基准
                setMinAndMax(ds, vd.Meta); // 设置最值

            string polygon = string.Format("", "");
                double[] tr = new double[6];
                ds.GetGeoTransform(tr);
                vd.Meta.resolution = string.Format("{0},{1}", tr[1], Math.Abs(tr[5])); // 分辨率

                if (tr[0] == 0.0 && tr[1] == 1.0 && tr[2] == 0.0 && tr[3] == 0.0 && tr[4] == 0.0 && tr[5] == 1.0) return;

                Geometry minPoint = GetMinPoint(ds);
                Geometry maxPoint = GetMaxPoint(ds);
                double xmin = minPoint.GetX(0);
                double ymin = minPoint.GetY(0);
                double xmax = maxPoint.GetX(0);
                double ymax = maxPoint.GetY(0);
                string geom = string.Format("ST_GeomFromText('POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))')", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin);

                vd.Meta.geom = geom; // 四至范围
            }
            catch (Exception ex)
            {
                LogOut.Error(ex.StackTrace);
            }
            finally
            {
                if (ds != null) ds.Dispose();
            }
        }

        // 从地理空间转换到像素空间
        private int[] Geo2ImageSpace(Dataset ds, double x, double y)
        {
            double[] transform = new double[6];
            ds.GetGeoTransform(transform); // 影像坐标变换参数

            int col = (int)((y * transform[1] - x * transform[4] + transform[0] * transform[4] - transform[3] * transform[1]) / (transform[5] * transform[1] - transform[2] * transform[4])); // 像素所在列
            int row = (int)((x - transform[0] - col * transform[2]) / transform[1]); // 像素所在行

            return new int[] { row, col };
        }

        // 从地理空间转换到像素空间
        private int[] Geo2ImageSpace(double[] transform, double x, double y)
        {
            int col = (int)((y * transform[1] - x * transform[4] + transform[0] * transform[4] - transform[3] * transform[1]) / (transform[5] * transform[1] - transform[2] * transform[4])); // 像素所在列
            int row = (int)((x - transform[0] - col * transform[2]) / transform[1]); // 像素所在行

            return new int[] { row, col };
        }

        // 从像素空间转换到地理空间
        private double[] ImageToGeoSpace(Dataset ds, int row, int col)
        {
            double[] transform = new double[6];
            ds.GetGeoTransform(transform);

            double x = transform[0] + row * transform[1] + col * transform[2];
            double y = transform[3] + row * transform[4] + col * transform[5];

            return new double[] { x, y };
        }

        #region GDAL切图
        /// <summary>
        /// 获取Dataset的最小点
        /// </summary>
        private Geometry GetMinPoint(Dataset ds)
        {
            double[] transform = new double[6];
            ds.GetGeoTransform(transform);

            string epsg = ds.GetSpatialRef().GetAuthorityCode(null);
            double xMin = transform[0];
            double yMin = transform[3];
            double yMin = transform[3] - ds.RasterYSize * transform[1];

            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(xMin, yMin, 0);
            if ("4326" == epsg)
            {
                point.AssignSpatialReference(sr4326);
                return point;
            }

            point.AssignSpatialReference(ds.GetSpatialRef());
            point.TransformTo(sr4326);
            return point;
            return Transform(ds, point);
        }

        /// <summary>
        /// 获取Dataset的最大点
        /// </summary>
        private Geometry GetMaxPoint(Dataset ds)
        {
            /**
             * transform[0] 左上角x坐标
             * transform[1] 东西方向分辨率
             * transform[2] 旋转角度, 0表示图像 "北方朝上"
             * 
             * transform[3] 左上角y坐标
             * transform[4] 旋转角度, 0表示图像 "北方朝上"
             * transform[5] 南北方向分辨率
             */
            double[] transform = new double[6];
            ds.GetGeoTransform(transform);

            string epsg = ds.GetSpatialRef().GetAuthorityCode(null);
            double xMax = transform[0] + (ds.RasterXSize * transform[1]);
            double yMax = transform[3] + (ds.RasterYSize * transform[1]);
            double yMax = transform[3];

            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(xMax, yMax, 0);
            if ("4326" == epsg)

            return Transform(ds, point);
        }

        /// <summary>
        /// 坐标转换
        /// </summary>
        private Geometry Transform(Dataset ds, Geometry point)
        {
            point.AssignSpatialReference(ds.GetSpatialRef());
            if (ds.GetSpatialRef().IsGeographic() > 0)
            {
                point.AssignSpatialReference(sr4326);
                return point;
            }

            point.AssignSpatialReference(ds.GetSpatialRef());
            point.TransformTo(sr4326);
            string srsName = ds.GetSpatialRef().GetName();
            if (srsName.Contains("CGCS2000"))
            {
                point.TransformTo(sr4490);
            }
            else
            {
                point.TransformTo(sr4326);
            }
            point.SwapXY();

            return point;
        }

        public void GenerateTiles()
        {
            string outPath = "D:\\xyz\\ce\\xyz2";
            string file = "D:\\xyz\\ce\\5_A1.tif";
            if (!File.Exists(file)) return;

            Stopwatch sw = new Stopwatch();
            sw.Start(); // 程序开始时间
            Dataset ds = null;
            try
            {
                ds = Gdal.Open(file, Access.GA_Update);
                if (null == ds) return;

                double[] transform = new double[6];
                ds.GetGeoTransform(transform);
                int rasterCount = ds.RasterCount;
                Console.WriteLine(string.Format("Origin = ({0}, {1})", transform[0], transform[3]));
                Console.WriteLine(string.Format("Pixel Size = ({0}, {1})", transform[1], transform[5]));

                // 4.1 首先获取原始影像的地理坐标范围
                Geometry minPoint = GetMinPoint(ds);
                Geometry maxPoint = GetMaxPoint(ds);
                double xmin = minPoint.GetX(0);
                double ymin = minPoint.GetY(0);
                double xmax = maxPoint.GetX(0);
                double ymax = maxPoint.GetY(0);
                Geometry imageBound = CreatePolygon(xmin, xmax, ymin, ymax, sr4326);

                // 4.2 获取原始影像的像素分辨率
                // 原始图像东西方向像素分辨率
                double src_w_e_pixel_resolution = (xmax - xmin) / ds.RasterXSize;
                // 原始图像南北方向像素分辨率
                double src_n_s_pixel_resolution = (ymax - ymin) / ds.RasterYSize;

                // 获取Band
                Band in_band1 = ds.GetRasterBand(1);
                Band in_band2 = ds.GetRasterBand(2);
                Band in_band3 = ds.GetRasterBand(3);
                //in_band1.Fill(0, 255); // GdalConst.GMF_NODATA
                in_band1.SetNoDataValue(0); // -9999
                in_band2.SetNoDataValue(0);
                in_band3.SetNoDataValue(0);

                int BufferSize = 256 * 256 * 5; // GdalConst.GDT_Int32
                for (int zoom = 10; zoom <= 18; zoom++)
                {
                    // 4.3 根据原始影像地理范围求解切片行列号  // 经纬度转瓦片编号
                    int tileRowMax = Lat2tile(ymin, zoom); // 纬度  -90 ->  90 lat
                    int tileRowMin = Lat2tile(ymax, zoom);
                    int tileColMin = Lon2tile(xmin, zoom); // 经度 -180 -> 180 lon
                    int tileColMax = Lon2tile(xmax, zoom);

                    Parallel.For(tileColMin, tileColMax + 1, (col) =>
                    {
                        Parallel.For(tileRowMin, tileRowMax + 1, (row) =>
                        {
                            // 4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集
                            double tempLatMin = Tile2lat(row + 1, zoom);
                            double tempLatMax = Tile2lat(row, zoom);
                            double tempLonMin = Tile2lon(col, zoom);
                            double tempLonMax = Tile2lon(col + 1, zoom);

                            Console.WriteLine(string.Format("{0}\\{1}\\{2}.png", zoom, col, row));
                            Geometry tileBound = CreatePolygon(tempLonMin, tempLonMax, tempLatMin, tempLatMax, sr4326);
                            Geometry intersect = tileBound.Intersection(imageBound);
                            if (null == intersect)
                            {
                                Console.WriteLine(string.Format("{0}\\{1}\\{2}.png，不存在", zoom, col, row));
                                return;
                            }
                            Envelope env = new Envelope();
                            intersect.GetEnvelope(env);

                            // 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)((env.MinX - xmin) / src_w_e_pixel_resolution);
                            int offset_y = (int)Math.Abs((env.MaxY - ymax) / src_n_s_pixel_resolution);

                            // 求在切图地理范围内的原始图像的像素大小
                            int block_xsize = (int)((env.MaxX - env.MinX) / src_w_e_pixel_resolution);
                            int block_ysize = (int)((env.MaxY - env.MinY) / src_n_s_pixel_resolution);

                            // 求原始图像在切片内的像素大小
                            int image_Xbuf = (int)Math.Ceiling((env.MaxX - env.MinX) / dst_w_e_pixel_resolution);
                            int image_Ybuf = (int)Math.Ceiling(Math.Abs((env.MaxY - env.MinY) / dst_n_s_pixel_resolution));

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

                            // 4.7 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
                            // 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
                            int[] band1BuffData = new int[BufferSize]; // 256 * 256 * GdalConst.GDT_Int32
                            int[] band2BuffData = new int[BufferSize];
                            int[] band3BuffData = new int[BufferSize];

                            try
                            {
                                in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);

                                // 4.8 将切片数据写入文件
                                // 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
                                OSGeo.GDAL.Driver memDriver = Gdal.GetDriverByName("MEM");
                                Dataset msmDS = memDriver.Create("msmDS", 256, 256, 4, DataType.GDT_Int32, null);
                                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[BufferSize];
                                for (int index = 0; index < alphaData.Length; index++)
                                {
                                    if (band1BuffData[index] > 0) alphaData[index] = 255;
                                }

                                // 写各个波段数据
                                dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData, image_Xbuf, image_Ybuf, 0, 0);

                                string path = Path.Combine(outPath, zoom.ToString(), col.ToString());
                                if (!Directory.Exists(path)) Directory.CreateDirectory(path);

                                // 使用PNG驱动将内存中的图像数组写入文件
                                string pngPath = path + "\\" + row + ".png";
                                OSGeo.GDAL.Driver pngDriver = Gdal.GetDriverByName("PNG");
                                Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS, 0, null, null, null);

                                msmDS.FlushCache();
                                pngDs.Dispose(); // 释放内存
                            }
                            catch (Exception ex)
                            {
                                Console.WriteLine(ex.Message);
                            }
                        });
                    });
                }
            }
            catch (Exception ex)
            {
                Console.WriteLine(ex.Message);
            }
            finally
            {
                if (null != ds) ds.Dispose();
                sw.Stop(); // 程序结束
                Console.WriteLine("耗时" + Math.Round(sw.ElapsedMilliseconds / 1000.0, 2) + "s"); // 秒
            }
        }

        public void GenerateTiles2()
        {
            string outPath = "D:\\xyz\\ce\\xyz";
            string file = "D:\\xyz\\ce\\5_A1.tif";
            //string file = "d:\\xyz\\dem\\dem\\33b.tif"; //string file = "D:\\xyz\\dq\\dq.vrt";
            if (!File.Exists(file)) return;

            Stopwatch sw = new Stopwatch();
            sw.Start(); // 程序开始时间
            Dataset ds = null;
            try
            {
                ds = Gdal.Open(file, Access.GA_Update);
                if (null == ds) return;

                double[] transform = new double[6];
                ds.GetGeoTransform(transform);
                int rasterCount = ds.RasterCount;
                Console.WriteLine(string.Format("Origin = ({0}, {1})", transform[0], transform[3]));
                Console.WriteLine(string.Format("Pixel Size = ({0}, {1})", transform[1], transform[5]));

                // 4.1 首先获取原始影像的地理坐标范围
                Geometry minPoint = GetMinPoint(ds);
                Geometry maxPoint = GetMaxPoint(ds);
                double xmin = minPoint.GetX(0);
                double ymin = minPoint.GetY(0);
                double xmax = maxPoint.GetX(0);
                double ymax = maxPoint.GetY(0);
                Geometry imageBound = CreatePolygon(xmin, xmax, ymin, ymax, sr4326);

                // 4.2 获取原始影像的像素分辨率
                // 原始图像东西方向像素分辨率
                double src_w_e_pixel_resolution = (xmax - xmin) / ds.RasterXSize;
                // 原始图像南北方向像素分辨率
                double src_n_s_pixel_resolution = (ymax - ymin) / ds.RasterYSize;

                // 获取Band
                Band in_band1 = ds.GetRasterBand(1);
                Band in_band2 = ds.GetRasterBand(2);
                Band in_band3 = ds.GetRasterBand(3);
                //in_band1.Fill(0, 255); // GdalConst.GMF_NODATA
                in_band1.SetNoDataValue(0); // -9999
                in_band2.SetNoDataValue(0);
                in_band3.SetNoDataValue(0);

                int BufferSize = 256 * 256 * 5; // GdalConst.GDT_Int32
                for (int zoom = 10; zoom <= 18; zoom++)
                {
                    // 4.3 根据原始影像地理范围求解切片行列号  // 经纬度转瓦片编号
                    int tileRowMax = Lat2tile(ymin, zoom); // 纬度  -90 ->  90 lat
                    int tileRowMin = Lat2tile(ymax, zoom);
                    int tileColMin = Lon2tile(xmin, zoom); // 经度 -180 -> 180 lon
                    int tileColMax = Lon2tile(xmax, zoom);
                    //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);

                    for (int col = tileColMin; col <= tileColMax; col++)
                    {
                        for (int row = tileRowMin; row <= tileRowMax; row++)
                        {
                            // 4.4 求原始影像地理范围与指定缩放级别指定行列号的切片交集
                            double tempLatMin = Tile2lat(row + 1, zoom);
                            double tempLatMax = Tile2lat(row, zoom);
                            double tempLonMin = Tile2lon(col, zoom);
                            double tempLonMax = Tile2lon(col + 1, zoom);
                            //double tempLatMin = tile2lat(row + 1, zoom);
                            //double tempLatMax = tile2lat(row, zoom);
                            //double tempLonMin = tile2lon(col, zoom);
                            //double tempLonMax = tile2lon(col + 1, zoom);

                            Console.WriteLine(string.Format("{0}\\{1}\\{2}.png", zoom, col, row));
                            Geometry tileBound = CreatePolygon(tempLonMin, tempLonMax, tempLatMin, tempLatMax, sr4326);
                            Geometry intersect = tileBound.Intersection(imageBound);
                            if (null == intersect)
                            {
                                Console.WriteLine(string.Format("{0}\\{1}\\{2}.png，不存在", zoom, col, row));
                                continue;
                            }
                            Envelope env = new Envelope();
                            intersect.GetEnvelope(env);

                            // 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)((env.MinX - xmin) / src_w_e_pixel_resolution);
                            int offset_y = (int)Math.Abs((env.MaxY - ymax) / src_n_s_pixel_resolution);

                            // 求在切图地理范围内的原始图像的像素大小
                            int block_xsize = (int)((env.MaxX - env.MinX) / src_w_e_pixel_resolution);
                            int block_ysize = (int)((env.MaxY - env.MinY) / src_n_s_pixel_resolution);

                            // 求原始图像在切片内的像素大小
                            int image_Xbuf = (int)Math.Ceiling((env.MaxX - env.MinX) / dst_w_e_pixel_resolution);
                            int image_Ybuf = (int)Math.Ceiling(Math.Abs((env.MaxY - env.MinY) / dst_n_s_pixel_resolution));

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

                            // 4.7 使用GDAL的ReadRaster方法对影像指定范围进行读取与压缩。
                            // 推荐在切片前建立原始影像的金字塔文件，ReadRaster在内部实现中可直接读取相应级别的金字塔文件，提高效率。
                            int[] band1BuffData = new int[BufferSize]; // 256 * 256 * GdalConst.GDT_Int32
                            int[] band2BuffData = new int[BufferSize];
                            int[] band3BuffData = new int[BufferSize];

                            try
                            {
                                //         ReadRaster(int xOff, int yOff, int xSize, int ySize, int[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace)
                                //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);
                                in_band1.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                in_band2.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                in_band3.ReadRaster(offset_x, offset_y, block_xsize, block_ysize, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);

                                // 4.8 将切片数据写入文件
                                // 使用gdal的MEM驱动在内存中创建一块区域存储图像数组
                                OSGeo.GDAL.Driver memDriver = Gdal.GetDriverByName("MEM");
                                Dataset msmDS = memDriver.Create("msmDS", 256, 256, 4, DataType.GDT_Int32, null);
                                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[BufferSize];
                                for (int index = 0; index < alphaData.Length; index++)
                                {
                                    if (band1BuffData[index] > 0)
                                    {
                                        alphaData[index] = 255;
                                    }
                                }

                                // 写各个波段数据
                                //         WriteRaster(int xOff, int yOff, int xSize, int ySize, int[] buffer, int buf_xSize, int buf_ySize, int pixelSpace, int lineSpace)
                                //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);
                                dstBand1.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band1BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                dstBand2.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band2BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                dstBand3.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, band3BuffData, image_Xbuf, image_Ybuf, 0, 0);
                                alphaBand.WriteRaster(imageOffsetX, imageOffsetY, image_Xbuf, image_Ybuf, alphaData, image_Xbuf, image_Ybuf, 0, 0);

                                //String pngPath = "D:\\xyz\\temp" + "\\" + zoom + "c" + col + "r" + row + ".png";
                                //Console.WriteLine("pngPath=" + pngPath);
                                string path = Path.Combine(outPath, zoom.ToString(), col.ToString());
                                if (!Directory.Exists(path)) Directory.CreateDirectory(path);

                                // 使用PNG驱动将内存中的图像数组写入文件
                                string pngPath = path + "\\" + row + ".png";
                                OSGeo.GDAL.Driver pngDriver = Gdal.GetDriverByName("PNG");
                                Dataset pngDs = pngDriver.CreateCopy(pngPath, msmDS, 0, null, null, null);

                                // 释放内存
                                msmDS.FlushCache();
                                pngDs.Dispose();
                            }
                            catch (Exception ex)
                            {
                                Console.WriteLine(ex.Message);
                            }
                        }
                    }
                }
            }
            catch (Exception ex)
            {
                Console.WriteLine(ex.Message);
            }
            finally
            {
                if (null != ds) ds.Dispose();
                sw.Stop(); // 程序结束
                Console.WriteLine("耗时" + Math.Round(sw.ElapsedMilliseconds / 1000.0, 2) + "s"); // 秒
            }
        }

        /// <summary>
@@ -548,63 +238,50 @@
        }

        /// <summary>
        /// 经度转瓦片编号
        /// 设置最值
        /// GDALRasterBand::GetHistogram ​​统计直方图
        /// </summary>
        public static int Lon2tile(double lon, int zoom)
        private void setMinAndMax(Dataset ds, SysMeta meta)
        {
            return (int)(Math.Floor((lon + 180) / 360 * Math.Pow(2, zoom)));
        }
            List<double> minList = new List<double>();
            List<double> maxList = new List<double>();

        /// <summary>
        /// 纬度转瓦片编号
        /// </summary>
        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)));
        }
            for (int i = 1; i <= ds.RasterCount; i++)
            {
                int hasval;
                double min, max;
                ds.GetRasterBand(i).GetMinimum(out min, out hasval);
                ds.GetRasterBand(i).GetMaximum(out max, out hasval);

        /// <summary>
        /// 瓦片编号转经度
        /// </summary>
        public static double Tile2lon(int col, int zoom)
        {
            return col / Math.Pow(2.0, zoom) * 360.0 - 180;
        }
                minList.Add(min);
                maxList.Add(max);
            }

        /// <summary>
        /// 瓦片编号转纬度
        /// </summary>
        public static double Tile2lat(int row, int zoom)
        {
            double n = Math.PI - (2.0 * Math.PI * row) / Math.Pow(2.0, zoom);
            return ToDegrees(Math.Atan(Math.Sinh(n)));
        }

        /// <summary>
        /// 常数e
        /// </summary>
        public static readonly double E = 2.7182818284590452354;

        /// <summary>
        /// 常数Pi
        /// </summary>
        public static readonly double PI = 3.14159265358979323846;

        /// <summary>
        /// 度转弧度
        /// </summary>
        public static double ToRadians(double angdeg)
        {
            return angdeg / 180.0 * Math.PI;
        }

        /// <summary>
        /// 弧度转度
        /// </summary>
        public static double ToDegrees(double angrad)
        {
            return angrad * 180.0 / Math.PI;
            meta.min = string.Join(",", minList.ToArray());
            meta.max = string.Join(",", maxList.ToArray());
        }
        #endregion

        public void CsTransform(double x, double y, int epsg)
        {
            SpatialReference srs = new SpatialReference(null);
            srs.ImportFromEPSG(epsg);

            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(x, y, 0);
            point.AssignSpatialReference(srs);

            if (srs.GetName().Contains("CGCS2000"))
            {
                point.TransformTo(sr4490);
            }
            else
            {
                point.TransformTo(sr4326);
            }
            point.SwapXY();

            double[] xy = new double[] { point.GetX(0), point.GetY(0) };
        }
    }
}