NDBGJK5.0/Kingo.Plugins/单图斑建库/Kingo.RuleCheck-陕西/CheckHelper/GeometryOperations.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geodatabase;
using OSGeo.OGR;

namespace Kingo.RuleCheck.CheckHelper
{
    /// <summary>
    /// 几何操作类
    /// </summary>
    public class GeometryOperations
    {
        /// <summary>
        /// 判断面是否自相交
        /// </summary>
        /// <param name="pgeo">图形</param>
        /// <returns>true 不相交，false 相交</returns>
        public static bool RuleTopoSingleCheckerIPolygon(IGeometry pgeo)
        {
            ITopologicalOperator3 topoOp3 = pgeo as ITopologicalOperator3;//Must be polygon
            topoOp3.IsKnownSimple_2 = false;
            //pTopoErrors = topoOp3.GetTopologyErrors();
            bool eanpBool = topoOp3.IsSimple;
            return eanpBool;
        }
        /// <summary>
        /// 合并满足指定条件查询的所有图形，返回合并后图形
        /// </summary>
        /// <param name="featureClass">要素类</param>
        /// <param name="filter">查询过滤器，IQueryFilter对象</param>
        /// <returns>返回合并后的几何图形</returns>
        public static IGeometry CreateGeometryFromEnumerator(IFeatureClass featureClass, IQueryFilter filter)
        {
            IEnumGeometry pEnumGeom = new EnumFeatureGeometry();
            IEnumGeometryBind pEnumGeomBind = pEnumGeom as IEnumGeometryBind;
            IGeometryFactory pGeomFactory = new GeometryEnvironment() as IGeometryFactory;
            pEnumGeomBind.BindGeometrySource(filter, featureClass);

            if (pEnumGeom == null) return null;

            IGeometry pGeo = pGeomFactory.CreateGeometryFromEnumerator(pEnumGeom);
            return pGeo;
        }

        /// <summary>
        /// 获取要素合并区域
        /// </summary>
        /// <param name="listFeature"></param>
        /// <returns></returns>
        /// <remarks></remarks>
        public static IGeometry GetFeatureUnionZoom(IList<IFeature> listFeature)
        {
            IGeometry m_unionGeo = null;
            if (m_unionGeo == null && listFeature != null && listFeature.Count != 0)
            {
                if (listFeature[0].ShapeCopy == null || listFeature[0].ShapeCopy.IsEmpty)
                {
                    return m_unionGeo;
                }
                IGeometry pGeo = listFeature[0].ShapeCopy;
                ITopologicalOperator pTopo = pGeo as ITopologicalOperator;

                for (int i = 1; i <= listFeature.Count - 1; i++)
                {
                    pGeo = pTopo.Union(listFeature[i].ShapeCopy);
                    pTopo = pGeo as ITopologicalOperator;
                }
                pGeo.SpatialReference = listFeature[0].Shape.SpatialReference;
                m_unionGeo = pGeo;
            }
            return m_unionGeo;
        }
        /// <summary>
        /// 获取面状要素边界缓冲区融合区域
        /// </summary>
        /// <param name="listFeature">要素集</param>
        /// <param name="pBuffer">缓冲设置</param>
        /// <returns></returns>
        /// <remarks></remarks>
        public static IGeometry GetFeatureBoundaryUnionZoom(IList<IFeature> listFeature, double pBuffer)
        {
            IGeometry m_unionGeo = null;
            if (m_unionGeo == null && listFeature != null && listFeature.Count != 0)
            {
                IGeometry pGeo = GetFeatureBoundary(listFeature[0], pBuffer);
                if (pGeo == null)
                    return m_unionGeo;
                IGeometry pGeoCurrent = null;
                ITopologicalOperator pTopo = pGeo as ITopologicalOperator;

                for (int i = 1; i <= listFeature.Count - 1; i++)
                {
                    pGeoCurrent = GetFeatureBoundary(listFeature[i], pBuffer);
                    if (pGeoCurrent == null)
                        return m_unionGeo;
                    pGeo = pTopo.Union(pGeoCurrent);
                    pTopo = pGeo as ITopologicalOperator;
                }
                pGeo.SpatialReference = listFeature[0].Shape.SpatialReference;
                m_unionGeo = pGeo;
            }
            return m_unionGeo;
        }
        /// <summary>
        /// 获取面状要素边界缓冲区域
        /// </summary>
        /// <param name="pFeature"></param>
        /// <param name="pBuffer"></param>
        /// <returns></returns>
        /// <remarks></remarks>
        public static IGeometry GetFeatureBoundary(IFeature pFeature, double pBuffer)
        {
            IGeometry pGeo = pFeature.ShapeCopy;
            if (pGeo == null)
                return null;
            if (pGeo.GeometryType == esriGeometryType.esriGeometryPolygon)
            {
                ITopologicalOperator pTopo = pGeo as ITopologicalOperator;
                return GetBuffer(pTopo.Boundary, pBuffer);
            }
            return null;
        }
        /// <summary>
        /// 创建面要素(单个面要素)
        /// </summary>
        /// <param name="pPointArray"></param>
        /// <param name="pSpatialReference"></param>
        /// <returns></returns>
        /// <remarks></remarks>
        public static IGeometry CreatePolygonFeatures(IPointArray pPointArray, ISpatialReference pSpatialReference = null)
        {
            IPointArray pValidPointArray = new PointArrayClass();
            //判断构造面的点是否超出数据集空间域
            for (int iIndex = 0; iIndex <= pPointArray.Count - 1; iIndex++)
            {
                if (pSpatialReference == null)
                {
                    IPoint pTempPoint = pPointArray.Element[iIndex];
                    pValidPointArray.Add(pTempPoint);
                }
                else
                {
                    if (ValidateXYDomain(pPointArray.Element[iIndex], pSpatialReference))
                    {
                        IPoint pTempPoint = pPointArray.Element[iIndex];
                        pValidPointArray.Add(pTempPoint);
                    }
                }
            }
            //若有超出数据集空间域的点，则不执行导入
            if (pPointArray.Count != pValidPointArray.Count)
            {
                return null;
            }
            else
            {
                //由文本文件中的点构造出单个面要素并导入
                //此处构造的面要素只考虑单个ring组成的polygon，还存在多个ring和岛环的情况
                IPointCollection pPolygon = new PolygonClass();
                for (int iIndex = 0; iIndex <= pValidPointArray.Count - 1; iIndex++)
                {
                    object Missing1 = Type.Missing;
                    object Missing2 = Type.Missing;
                    pPolygon.AddPoint(pValidPointArray.Element[iIndex], ref Missing1, ref Missing2);
                }
                pPolygon.AddPoint(pValidPointArray.Element[0]);

                return pPolygon as IPolygon;
            }
        }
        /// <summary>
        /// 创建线要素(单个线要素)
        /// </summary>
        /// <param name="pPointArray"></param>
        /// <param name="pSpatialReference"></param>
        /// <remarks></remarks>
        public static IGeometry CreatePolyLineFeatures(IPointArray pPointArray, ISpatialReference pSpatialReference = null)
        {
            IPointArray pValidPointArray = new PointArrayClass();
            //判断构造线的点是否超出数据集空间域
            for (int iIndex = 0; iIndex <= pPointArray.Count - 1; iIndex++)
            {
                if (pSpatialReference == null)
                {
                    IPoint pTempPoint = pPointArray.Element[iIndex];
                    pValidPointArray.Add(pTempPoint);
                }
                else
                {
                    if (ValidateXYDomain(pPointArray.Element[iIndex], pSpatialReference))
                    {
                        IPoint pTempPoint = pPointArray.Element[iIndex];
                        pValidPointArray.Add(pTempPoint);
                    }
                }
            }
            //若有超出数据集空间域的点，则不执行导入
            if (pPointArray.Count != pValidPointArray.Count)
            {
                return null;
            }
            else
            {
                //由文本文件中的点构造出单个线要素并导入
                IPointCollection pPolyline = new PolylineClass();
                for (int iIndex = 0; iIndex <= pValidPointArray.Count - 1; iIndex++)
                {
                    object Missing1 = Type.Missing;
                    object Missing2 = Type.Missing;
                    pPolyline.AddPoint(pValidPointArray.Element[iIndex], ref Missing1, ref Missing2);
                }
                return pPolyline as IGeometry;
            }
        }
        /// <summary>
        /// 验证几何图形坐标是否超出数据集空间域
        /// </summary>
        /// <param name="pGeo">几何图形</param>
        /// <param name="pSpatialReference"></param>
        /// <remarks></remarks>
        private static bool ValidateXYDomain(IGeometry pGeo, ISpatialReference pSpatialReference)
        {
            double dXMin, dXMax, dYMin, dYMax;
            pSpatialReference.GetDomain(out dXMin, out dXMax, out dYMin, out dYMax);
            if (pGeo.Envelope.XMin >= dXMin && pGeo.Envelope.YMin >= dYMin &&
                pGeo.Envelope.XMax <= dXMax && pGeo.Envelope.YMax <= dYMax)
            {
                return true;
            }
            return false;
        }

        /// <summary>
        /// 获取缓冲值
        /// </summary>
        /// <param name="curGeo">缓冲要素</param>
        /// <param name="bufferValue">缓冲值(可选)</param>
        /// <returns></returns>
        public static IGeometry GetBuffer(IGeometry curGeo, Double bufferValue)
        {
            if (bufferValue == 0) return curGeo;
            IBufferConstruction pBufferConstrucion = new BufferConstruction();
            IGeometry pBuffer = null;
            for (Double dValue = bufferValue; dValue < bufferValue + 0.3; dValue = dValue + 0.01)
            {
                try
                {
                    pBuffer = pBufferConstrucion.Buffer(curGeo, dValue);
                    System.Runtime.InteropServices.Marshal.ReleaseComObject(pBufferConstrucion);
                    return pBuffer;
                }
                catch (Exception ex)
                {
                    Console.WriteLine(ex.ToString());
                }
            }
            return curGeo;
        }


        /// <summary>
        /// 根据宽度计算角平分线上的点（三点确定一个角）
        /// </summary>
        /// <param name="p0">第一个点</param>
        /// <param name="p1">第二个点</param>
        /// <param name="p2">第三个点</param>
        /// <param name="distance">宽度</param>
        /// <returns>角平分线上的点</returns>
        public static IPoint getParallelPoint(IPoint p0, IPoint p1, IPoint p2, double distance)
        {
            var y12 = p2.Y - p1.Y;
            var x12 = p2.X - p1.X;
            var y01 = p1.Y - p0.Y;
            var x01 = p1.X - p0.X;
            double a, b;
            if (x12 == 0)
            {
                a = Math.PI / 2;
                if (y12 < 0)
                    a = -a;
            }
            else
            {
                a = Math.Atan(y12 / x12);
            }
            if (x01 == 0)
            {
                b = Math.PI / 2;
                if (y01 < 0)
                    b = -b;
            }
            else
            {
                b = Math.Atan(y01 / x01);
            }
            // 关键核心处
            if (p2.X < p1.X)
            {
                a += Math.PI;
            }
            if (p1.X < p0.X)
            {
                b += Math.PI;
            }
            var k = (b - a - Math.PI) / 2;
            var r = a + k;
            var d = distance / Math.Sin(k);
            var sinr = Math.Sin(r);
            var cosr = Math.Cos(r);
            var d_lat = p1.Y + d * sinr;
            var d_lng = p1.X + d * cosr;
            return new PointClass() { X = d_lng, Y = d_lat };
        }

        /// <summary>
        /// 已知点A，根据B点据A点的距离，和方位，求B点的坐标
        /// </summary>
        /// <param name="A">已知点A</param>
        /// <param name="distance"> B点到A点的距离</param>
        /// <param name="angle">B点相对于A点的方位角度</param>
        /// <returns>B点的坐标</returns>
        public static IPoint GetPoint(IPoint A, double distance, double angle)
        {
            //int nCoordSystem = GetCoordSystem();
            IPoint point = null;
            //if (A.X <= 180 && A.Y <= 90)
            //白明雅 2018-12-19
            if (A.X <= distance && A.Y <= distance)
            {
                //LongLat ll1 = LongLatHelper.GetLongLat(new LongLat(A.X, A.Y), distance, angle);
                //point = new Vertex2D(ll1.m_Longitude, ll1.m_Latitude);
            }
            else
            {
                double dbθ = (180 - angle) * Math.PI / 180;
                point = new PointClass();
                point.X = A.X + distance * Math.Sin(dbθ);
                point.Y = A.Y - distance * Math.Cos(dbθ);
            }
            return point;
        }
        /// <summary>
        /// 计算两点之间的角度
        /// </summary>
        /// <param name="p1">点1</param>
        /// <param name="p2">点2</param>
        /// <returns></returns>
        public static double getAngle(IPoint p1, IPoint p2)
        {
            //两点的x、y值
            double x = p2.X - p1.X;
            double y = p2.Y - p1.Y;
            double hypotenuse = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2));
            //斜边长度
            double cos = x / hypotenuse;
            double radian = Math.Acos(cos);
            //求出弧度
            double angle = 180 / (Math.PI / radian);
            //用弧度算出角度        
            if (y < 0)
            {
                angle = -angle;
            }
            else if ((y == 0) && (x < 0))
            {
                angle = 180;
            }
            return angle;
        }
        /// <summary>
        /// 计算平行线上的点集合
        /// </summary>
        /// <param name="points">参考线上的点集合</param>
        /// <param name="pDistance">距离</param>
        /// <param name="pOutPoints">平行线上的点集合</param>
        public static void PointsFromParalleLine(IPointCollection points, double pDistance, ref IPointCollection pOutPoints)
        {
            if (pOutPoints == null)
                return;
            for (int i = 0; i < points.PointCount; i++)
            {
                double angle = 0;
                IPoint lp = points.get_Point(i);
                if (i == points.PointCount - 1)
                {
                    angle = GeometryOperations.getAngle(points.Point[i - 1], points.Point[i]);
                }
                else
                {
                    angle = GeometryOperations.getAngle(points.Point[i], points.Point[i + 1]);
                }
                lp = GeometryOperations.GetPoint(points.Point[i], pDistance, -angle);
                if (i > 1)
                {
                    lp = GeometryOperations.getParallelPoint(points.get_Point(i - 2), points.get_Point(i - 1), points.get_Point(i), pDistance);
                    pOutPoints.UpdatePoint(i - 1, lp);
                }
                if (i == points.PointCount - 1 || i <= 1)
                {
                    lp = GeometryOperations.GetPoint(points.Point[i], pDistance, -angle);
                    pOutPoints.UpdatePoint(i, lp);
                }
            }
        }
        /// <summary>
        /// 获取最小夹角
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public static double GetMinAngle2(ESRI.ArcGIS.Geometry.IGeometry geometry)
        {
            ESRI.ArcGIS.Geometry.IGeometryFactory3 factory = new ESRI.ArcGIS.Geometry.GeometryEnvironment() as ESRI.ArcGIS.Geometry.IGeometryFactory3;
            byte[] b = factory.CreateWkbVariantFromGeometry(geometry) as byte[];
            Geometry gdalGeometry = Geometry.CreateFromWkb(b);
            return GetAngle(gdalGeometry);
        }

        public static double GetAngle(Geometry geometry)
        {
            string wkt = "";
            geometry.ExportToWkt(out wkt);
            double rstAngle = -1;
            double sumAngle = 0;
            if (wkt.StartsWith("POLYGON"))
            {

                //StringBuilder sb = new StringBuilder("POLYGON (");
                string wktN = wkt.Replace("POLYGON ((", "").Replace("))", "");
                string[] strs = wktN.Replace("),(", "@").Split('@');
                for (int i = 0; i < strs.Length; i++)
                {

                    string[] strsP = strs[i].Split(',');
                    for (int j = 0; j < strsP.Length; j++)
                    {
                        string[] strF = strsP[j].ToString().Split(' ');
                        string[] strC = null;
                        string[] strS = null;
                        if (j == strsP.Length - 2)
                        {
                            strC = strsP[j + 1].ToString().Split(' ');
                            strS = strsP[0].ToString().Split(' ');
                            if (strC[0] == strS[0] && strC[1] == strS[1])
                            {
                                strS = strsP[1].ToString().Split(' ');
                            }
                        }
                        else if (j == strsP.Length - 1)
                        {
                            strC = strsP[0].ToString().Split(' ');
                            strS = strsP[1].ToString().Split(' ');
                            continue;
                        }
                        else
                        {
                            strC = strsP[j + 1].ToString().Split(' ');
                            strS = strsP[j + 2].ToString().Split(' ');
                        }
                        double angle = GetAngle(Convert.ToDouble(strC[0]), Convert.ToDouble(strC[1]), Convert.ToDouble(strF[0]), Convert.ToDouble(strF[1]), Convert.ToDouble(strS[0]), Convert.ToDouble(strS[1]));
                        sumAngle += angle;
                        if (rstAngle == -1)
                        {
                            rstAngle = angle;
                        }
                        else
                        {
                            if (angle < rstAngle)
                            {
                                rstAngle = angle;
                            }
                        }
                    }

                }

                return rstAngle;

            }
            return 0;
        }
        public static double GetAngle(double cenX, double cenY, double firstX, double firstY, double secondX, double secondY)
        {
            double ma_x = firstX - cenX;
            double ma_y = firstY - cenY;
            double mb_x = secondX - cenX;
            double mb_y = secondY - cenY;
            double v1 = (ma_x * mb_x) + (ma_y * mb_y);
            double ma_val = Math.Sqrt(ma_x * ma_x + ma_y * ma_y);
            double mb_val = Math.Sqrt(mb_x * mb_x + mb_y * mb_y);
            double cosM = v1 / (ma_val * mb_val);
            double angleAMB = Math.Acos(cosM) * 180 / Math.PI;
            return angleAMB;
        }
        public static double GetAngle(IPoint cenPoint, IPoint firstPoint, IPoint secondPoint)
        {
            double ma_x = firstPoint.X - cenPoint.X;
            double ma_y = firstPoint.Y - cenPoint.Y;
            double mb_x = secondPoint.X - cenPoint.X;
            double mb_y = secondPoint.Y - cenPoint.Y;
            double v1 = (ma_x * mb_x) + (ma_y * mb_y);
            double ma_val = Math.Sqrt(ma_x * ma_x + ma_y * ma_y);
            double mb_val = Math.Sqrt(mb_x * mb_x + mb_y * mb_y);
            if (ma_val * mb_val == 0)
            {
                return -1;
            }
            double cosM = v1 / (ma_val * mb_val);
            double angleAMB = Math.Acos(cosM) * 180 / Math.PI;
            return angleAMB;
        }
        public static double GetMinAngle(IGeometry pGeo)
        {
            double rstAngle = -1;
            try
            {
                if (pGeo == null || pGeo.IsEmpty)
                    return rstAngle;
                IPolygon4 poly4 = pGeo as IPolygon4;
                ITopologicalOperator topo = poly4 as ITopologicalOperator;
                if (topo != null)
                {
                    topo.Simplify();
                }
                GeometryBag geoBag = poly4.ExteriorRingBag as GeometryBag;
                if (geoBag == null) return rstAngle;
                IGeometryCollection geoCollection = geoBag as IGeometryCollection;
                List<IGeometry> rings = new List<IGeometry>();
                for (int j = 0; j < geoCollection.GeometryCount; j++)
                {
                    IGeometry geo = geoCollection.get_Geometry(j);
                    rings.Add(geo);
                    //内环图形
                    IGeometryBag InteriorBag = (pGeo as ESRI.ArcGIS.Geometry.IPolygon4).get_InteriorRingBag(geo as IRing);
                    if (InteriorBag != null)
                    {
                        IGeometryCollection InteriorRingGeometryCollection = InteriorBag as IGeometryCollection;
                        for (int IR = 0; IR < InteriorRingGeometryCollection.GeometryCount; IR++)
                        {
                            rings.Add(InteriorRingGeometryCollection.get_Geometry(IR));
                        }
                    }
                }
                foreach (IGeometry ring in rings)
                {
                    if (ring.IsEmpty) continue;
                    IPointCollection points = ring as IPointCollection;
                    int num = points.PointCount - 1;
                    for (int i = 0; i < num; i++)
                    {
                        IPoint p1 = null;
                        IPoint p2 = points.get_Point(i);
                        IPoint p3 = null;
                        if (i == 0)
                        {
                            p1 = points.get_Point(num - 1);
                            p3 = points.get_Point(i + 1);
                        }
                        else if (i == num - 1)
                        {
                            p1 = points.get_Point(i - 1);
                            p3 = points.get_Point(0);
                        }
                        else
                        {
                            p1 = points.get_Point(i - 1);
                            p3 = points.get_Point(i + 1);
                        }
                        double angle = GetAngle(p2, p1, p3);
                        if (rstAngle == -1)
                        {
                            rstAngle = angle;
                        }
                        else
                        {
                            if (rstAngle > angle)
                            {
                                rstAngle = angle;
                            }
                        }
                    }
                }
            }
            catch (Exception ex)
            {
                throw ex;
            }
            return rstAngle;
        }
    }
}