基于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操作。