1

13693261870

2023-08-19 59e7594e461c011d582bb6fb3ea6b150e8be6d54

1

 DataLoader/CS/GdalHelper.cs

@@ -137,474 +137,5 @@

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

        #region GDAL切图
        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];

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

        private Geometry GetMaxPoint(Dataset ds)
        {
            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]);

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

            point.AssignSpatialReference(ds.GetSpatialRef());
            point.TransformTo(sr4326);
            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>
        /// 创建多边形
        /// </summary>
        public Geometry CreatePolygon(double xmin, double xmax, double ymin, double ymax, SpatialReference sr)
        {
            string kwt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin);
            Geometry geo = Geometry.CreateFromWkt(kwt);
            geo.AssignSpatialReference(sr);

            return geo;
        }

        /// <summary>
        /// 经度转瓦片编号
        /// </summary>
        public static int Lon2tile(double lon, int zoom)
        {
            return (int)(Math.Floor((lon + 180) / 360 * Math.Pow(2, zoom)));
        }

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

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

        /// <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;
        }
        #endregion
    }
}

 DataLoader/MainWindow.xaml.cs

@@ -88,9 +88,8 @@
            //this.Hide();
            //win.Show();

            //SetLoginInfo();
            SetLoginInfo();
            //GdalHelper.Instance.ReadTiff();
            GdalHelper.Instance.GenerateTiles();
        }

        public void SetLoginInfo()