1

13693261870

2023-09-20 124486dbffb46b6a80008b0ba591c254643de0ad

 DataLoader/CS/GdalHelper.cs

@@ -108,23 +108,20 @@
            try
            {
                ds = Gdal.Open(vd.FilePath, Access.GA_ReadOnly);
                if (null == ds || ds.RasterCount == 0) return;
                if (null == ds || ds.RasterCount == 0 || null == ds.GetSpatialRef()) return;

                vd.Meta.gridsize = string.Format("{0},{1}", ds.RasterXSize, ds.RasterYSize); // 行列数

                SpatialReference sr = ds.GetSpatialRef();
                if (sr != null)
                vd.Meta.coor_sys = sr.GetName(); // 坐标系统
                if (MOON200 == vd.Meta.coor_sys)
                {
                    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编码
                    }
                    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编码
                }

                vd.Meta.bands = ds.RasterCount.ToString(); // 波段数
@@ -132,12 +129,13 @@
                ColorTable ct = ds.GetRasterBand(1).GetRasterColorTable();
                vd.Meta.ct = null == ct ? null : ct.ToString(); // 数据颜色表
                vd.Meta.h_datum = null; // 高程基准
                setMinAndMax(ds, vd.Meta); // 设置最值

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

                if (!EPSGS.Contains(vd.Meta.epsg)) return;
                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);
@@ -145,9 +143,9 @@
                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);
                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 = wkt; // 四至范围
                vd.Meta.geom = geom; // 四至范围
            }
            catch (Exception ex)
            {
@@ -167,14 +165,13 @@
            double[] transform = new double[6];
            ds.GetGeoTransform(transform);

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

            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(xMin, yMin, 0);

            return Transform(ds, point, epsg);
            return Transform(ds, point);
        }

        /// <summary>
@@ -182,42 +179,40 @@
        /// </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];

            Geometry point = new Geometry(wkbGeometryType.wkbPoint);
            point.AddPoint(xMax, yMax, 0);

            return Transform(ds, point, epsg);
            return Transform(ds, point);
        }

        /// <summary>
        /// 坐标转换
        /// </summary>
        private Geometry Transform(Dataset ds, Geometry point, string epsg)
        private Geometry Transform(Dataset ds, Geometry point)
        {
            if (string.IsNullOrEmpty(epsg))
            point.AssignSpatialReference(ds.GetSpatialRef());
            if (ds.GetSpatialRef().IsGeographic() > 0)
            {
                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"))
            string srsName = ds.GetSpatialRef().GetName();
            if (srsName.Contains("CGCS2000"))
            {
                point.TransformTo(sr4490);
            }
@@ -225,6 +220,7 @@
            {
                point.TransformTo(sr4326);
            }
            point.SwapXY();

            return point;
        }
@@ -240,6 +236,52 @@

            return geo;
        }

        /// <summary>
        /// 设置最值
        /// GDALRasterBand::GetHistogram ​​统计直方图
        /// </summary>
        private void setMinAndMax(Dataset ds, SysMeta meta)
        {
            List<double> minList = new List<double>();
            List<double> maxList = new List<double>();

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

                minList.Add(min);
                maxList.Add(max);
            }

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