标签：case 基于 fs Features 拓扑 dotspatial longitude KnownCoordinateSystems Reproject
基于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操作。
标签：case,基于,fs,Features,拓扑,dotspatial,longitude,KnownCoordinateSystems,Reproject 来源： https://blog.csdn.net/bc6666/article/details/113704518
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。