基于dotspatial的拓扑检查

华坚成

2023-12-01

基于dotspatial的拓扑检查

前言

主要针对于不使用AE的用户（如果想使用最新版的AE至少要花几万块钱的授权费或者安装一大堆的破解版ArcGIS环境，因此我打算基于OpenGIS组件开发，解除对AE的依赖）。其中在进行对.shp文件图斑叠置检查也就是**Overlaps（）**函数执行时，发现OpenGIS组件库中无论是GEOS、Dotspatial、MapWinGIS等，结果都与ArcGIS的叠置结果不同。经过研究和测试，发现了这个问题是由坐标系和投影坐标系引起的。长话短说，接下来开始介绍我的解决方案。

技术标准

我这边是通过两个图斑相交部分的面积s>100平方米或者s>min(图斑1面积*10%，图斑2面积*10%)来判断两个图斑是否叠置。（当然两个图斑进行好投影坐标系的转换也可以进行叠置分析）

逻辑

创建图斑集合并判断shp文件的坐标系类型是大地2000还是西安80

           DotSpatial.Data.IFeatureSet fs =DotSpatial.Data.FeatureSet.OpenFile(shpfilepath);//shpfilepath是.shp文件所在路径
           string type = fs.Projection.ToEsriString();
           if (type.Contains("CGCS2000"))
                type = "CGCS2000";
            else
                type = "Xian1980";

获取相交的两个图斑对象以及相交部分的图斑对象

 DotSpatial.Topology.Geometry ge1 = (DotSpatial.Topology.Geometry)fs.Features[0].BasicGeometry;
 DotSpatial.Topology.Geometry ge2 = (DotSpatial.Topology.Geometry)fs.Features[1].BasicGeometry;            
 DotSpatial.Topology.Geometry ge3 = (DotSpatial.Topology.Geometry)ge1.Intersection(ge2);

计算三者的面积，使用getArea（）函数 一定要看注释

            IFeatureSet fs =new FeatureSet();
            fs.AddFeature(ge);
            fs.Projection = KnownCoordinateSystems.Geographic.World.WGS1984;//设置初始的平面投影坐标系
            double longitude = ge.Centroid.X;//获取经度 为了通过不同的度带来转化不同的坐标系
            int count = getProjectCodeByLongitude(longitude);//获取不同度带对应的count号 方便下述坐标转换
            string path = System.Environment.CurrentDirectory
            +"\\CommonClass\\CGCS2000projected\\CGCS2000_3_Degree_GK_CM_";//dotspatial中没有大地2000的投影坐标系 自己引入大地2000的路径
            double area=-1;
            if (type== "CGCS2000" && count != -1)
                count += 100;
            switch (count)
            {
                case -1:
                    return -1;                 
                case 0:
                   fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM75E); //转换坐标系为西安80的 75E的度带 该坐标系为dotspatial自带的。
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));//计算面积公式
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);//！！！重点 必须把坐标系重新设置为平面坐标系 不然后面相交部分面积计算会出错 怀疑是Reproject()执行的是全局坐标系转换。 
                    break;
                case 1:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM78E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 2:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM81E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 3:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM84E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 4:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM87E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 5:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM90E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 6:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM93E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 7:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM96E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 8:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM99E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 9:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM102E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 10:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM105E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 11:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM108E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 12:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM111E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 13:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM114E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 14:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM117E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 15:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM120E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 16:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM123E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 17:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM126E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 18:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM129E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 19:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM132E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 20:
                    fs.Reproject(KnownCoordinateSystems.Projected.KrugerXian1980.Xian19803DegreeGKCM135E);
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 100:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path+"75E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 101:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "78E.prj")));//这个是转换大地2000坐标 通过自己定义坐标系的文件 如果转为其他坐标系也可以
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));//计算面积
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);//重点！！必须设置回来
                    break;
                case 102:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "81E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 103:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "83.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 104:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "87E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 105:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "90E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 106:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "93E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 107:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "96E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 108:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "99E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 109:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "102E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 110:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "105E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 111:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "108E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 112:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "111E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 113:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "114E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 114:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "117E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 115:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "120E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 116:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "123E.prj")));
                    break;
                case 117:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "126E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 118:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "129E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 119:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "132E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
                case 120:
                    fs.Reproject(ProjectionInfo.FromEsriString(System.IO.File.ReadAllText(path + "135E.prj")));
                    area = Math.Abs(CgAlgorithms.SignedArea(fs.Features[0].Coordinates));
                    fs.Reproject(KnownCoordinateSystems.Geographic.World.WGS1984);
                    break;
            }
return area;

补充将度带转换为count的函数方便计算面积公式中swich语句的使用

  public static int getProjectCodeByLongitude(double longitude)
        {
            int count = 0;
            if (longitude > 136.5 || longitude < 73.5)
            {
                return -1;
            }
            else if (longitude >= 73.5 && longitude < 76.5)
            {
                count = 0;
            }
            else if (longitude >= 76.5 && longitude < 79.5)
            {
                count = 1;
            }
            else if (longitude >= 79.5 && longitude < 82.5)
            {
                count = 2;
            }
            else if (longitude >= 82.5 && longitude < 85.5)
            {
                count = 3;
            }
            else if (longitude >= 85.5 && longitude < 88.5)
            {
                count = 4;
            }
            else if (longitude >= 88.5 && longitude < 91.5)
            {
                count = 5;
            }
            else if (longitude >= 91.5 && longitude < 94.5)
            {
                count = 6;
            }
            else if (longitude >= 94.5 && longitude < 97.5)
            {
                count = 7;
            }
            else if (longitude >= 97.5 && longitude < 100.5)
            {
                count = 8;
            }
            else if (longitude >= 100.5 && longitude < 103.5)
            {
                count = 9;
            }
            else if (longitude >= 103.5 && longitude < 106.5)
            {
                count = 10;
            }
            else if (longitude >= 106.5 && longitude < 109.5)
            {
                count = 11;
            }
            else if (longitude >= 109.5 && longitude < 112.5)
            {
                count = 12;
            }
            else if (longitude >= 112.5 && longitude < 115.5)
            {
                count = 13;
            }
            else if (longitude >= 115.5 && longitude < 118.5)
            {
                count = 14;
            }
            else if (longitude >= 118.5 && longitude < 121.5)
            {
                count = 15;
            }
            else if (longitude >= 121.5 && longitude < 124.5)
            {
                count = 16;
            }
            else if (longitude >= 124.5 && longitude < 127.5)
            {
                count = 17;
            }
            else if (longitude >= 127.5 && longitude < 130.5)
            {
                count = 18;
            }
            else if (longitude >= 130.5 && longitude < 133.5)
            {
                count = 19;
            }
            else if (longitude >= 133.5 && longitude < 136.5)
            {
                count = 20;
            }
            return count;
        }

结尾

将各个图斑的面积计算出来然后按照技术标准来计算是否叠置。

最后经测试，结果与ArcGIS保持一致（面积计算出的单位均是平方米）

猜测：转换坐标系后直接用dotspatial里面的overlaps()函数应该也是可以的由于我有了判断叠置的技术标准因此没有继续测试猜测的内容。

以上分割线

欢迎大家来交流如何使用Dotspatial组件完成更多的GIS操作。

基于dotspatial的拓扑检查

基于dotspatial的拓扑检查

前言

技术标准

逻辑

结尾

相关阅读

相关文章

相关问答

相关文档