添加读取栅格信息的接口

13693261870

2023-08-23 a955e13d33579ed03ffdc8ce2f09671425a89920

添加读取栅格信息的接口

 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,6 +79,9 @@
            }
        }

        /// <summary>
        /// 注册GDAL
        /// </summary>
        public void RegisterGDal()
        {
            string gdalData = Path.Combine(Tools.BaseDir, "gdal-data");
@@ -85,57 +98,147 @@
        }
        #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) 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();
                if (sr != null)
                {
                    vd.Meta.coor_sys = sr.GetName(); // 坐标系统
                    if (MOON200 == vd.Meta.coor_sys)
                    {
                        vd.Meta.epsg = "ESRI:104903"; // EPSG编码
                    }
                    else
                    {
                        string code = sr.GetAuthorityCode(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); // 数据类型
                vd.Meta.ct = ds.GetRasterBand(1).GetRasterColorTable().ToString(); // 数据颜色表
                vd.Meta.h_datum = null; // 高程基准

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

                if (!EPSGS.Contains(vd.Meta.epsg)) 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 wkt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{0} {1}))", xmin, ymax, xmax, ymax, xmax, ymin, xmin, ymin);

                vd.Meta.geom = wkt; // 四至范围
            }
            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)
        /// <summary>
        /// 获取Dataset的最小点
        /// </summary>
        private Geometry GetMinPoint(Dataset ds)
        {
            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];
            string epsg = ds.GetSpatialRef().GetAuthorityCode(null);
            double xMin = transform[0];
            double yMin = transform[3];

            return new double[] { x, y };
            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(xMin, yMin, 0);

            return Transform(ds, point, epsg);
        }

        /// <summary>
        /// 获取Dataset的最大点
        /// </summary>
        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);

            return Transform(ds, point, epsg);
        }

        /// <summary>
        /// 坐标转换
        /// </summary>
        private Geometry Transform(Dataset ds, Geometry point, string epsg)
        {
            if (string.IsNullOrEmpty(epsg))
            {
                point.AssignSpatialReference(sr104903);
                return point;
            }
            if ("4326" == epsg)
            {
                point.AssignSpatialReference(sr4326);
                return point;
            }
            if ("4490" == epsg)
            {
                point.AssignSpatialReference(sr4490);
                return point;
            }

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

            return point;
        }

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

 DataLoader/CS/Importor.cs

@@ -239,7 +239,7 @@
            {
                try
                {
                    vd.Status = "生成MD5...";
                    vd.Status = "生成MD5码...";
                    vd.Meta.guid = GetFilesMD5(vd);

                    if (!string.IsNullOrEmpty(vd.Meta.guid) && DBHelper.IsFileExists(vd.Meta.guid))
@@ -248,10 +248,13 @@
                        return;
                    }

                    vd.Status = "读取栅格信息...";
                    ReadRasterInfo(vd);

                    vd.Status = "复制文件...";
                    CopyFiles(vd);

                    vd.Status = "数据入库...";
                    vd.Status = "元数据入库...";
                    int id = DBHelper.InsertMeta(vd.Meta);
                    vd.Status = id > 0 ? "成功。" : "失败！";
                }
@@ -315,6 +318,16 @@
        }

        /// <summary>
        /// 读取栅格信息
        /// </summary>
        private static void ReadRasterInfo(ViewData vd)
        {
            if (!StaticData.RASTER_EXT.Contains(vd.Ext) || !File.Exists(vd.FilePath)) return;

            GdalHelper.Instance.ReadRasterInfo(vd);
        }

        /// <summary>
        /// 复制文件
        /// </summary>
        private static void CopyFiles(ViewData vd)