基于dotspatial的拓扑检查

基于dotspatial的拓扑检查

前言

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

技术标准

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

逻辑

  1. 创建图斑集合并判断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";
  1. 获取相交的两个图斑对象以及相交部分的图斑对象

     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);
    
  2. 计算三者的面积,使用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;   
            
  1. 补充将度带转换为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操作。

上一篇:将多个月的nc数据文件合成一个(月平均)


下一篇:微信小程序之高德地图多点路线规划