GIS面积计算问题
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了GIS面积计算问题相关的知识，希望对你有一定的参考价值。
好长时间不更新了，今天说点干货，项目用到的。
1、项目中要用到计算面积的，根据火星坐标；
2、百度找了各种面积计算，google了半天，也没发现那个比较准确；
直接说干货吧。咱也高大上一会，用 ArcGIS
需要的dll有 ESRI.ArcGIS.Client.Toolkit.dll 和 ESRI.ArcGIS.Client.dll 这两个dll即可，不需要那么多的dll，可到官网下载，也可以到dll之家下载（http://www.dllzj.com/ESRI.ArcGIS.Client.dll/）
计算面积的代码，这个代码比较简单
GeometryService geometryService = new GeometryService("http://tasks.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer");
Polygon polygon = new Polygon();
polygon.SpatialReference = new SpatialReference(4326);//球面坐标
PointCollection GPSpoints = new PointCollection();//保存84坐标的集合
for (int i = 0; i < arraytemp.Count; i++)
{
JArray arraystr = (JArray)JsonConvert.DeserializeObject(arraytemp[i].ToString());
GPSpoints.Add(new MapPoint(Convert.ToDouble(arraystr[0].ToString()), Convert.ToDouble(arraystr[1].ToString())));
}
polygon.Rings.Add(GPSpoints);
Graphic graphic = new Graphic();
graphic.Geometry = polygon;
IList<Graphic> graphicList = new List<Graphic>(); //构建参数
graphicList.Add(graphic);
graphicList = geometryService.Project(graphicList, new SpatialReference(102100));//执行转换
var args = geometryService.AreasAndLengths(graphicList, LinearUnit.Meter, LinearUnit.Meter, CalculationType.Geodesic);//Geodesic 采用球面面积计算方式 102113
area = args.Areas[0];
说下其中的注意点，其中4326指的84坐标，但是国内的地图一般都是火星坐标（腾讯和高德，百度更闹心，直接是自己的坐标），就是加了偏移量的，具体的84坐标和火星坐标的转化，这个可能比较常见，代码如下
public static string BAIDU_LBS_TYPE = "bd09ll";
public static double pi = 3.1415926535897932384626;
public static double a = 6378245.0;
public static double ee = 0.00669342162296594323;
/**
* 84 to 火星坐标系 (GCJ-02) World Geodetic System ==> Mars Geodetic System
*
* @param lat
* @param lon
* @return
*/
public static Gps gps84_To_Gcj02(double lat, double lon)
{
if (outOfChina(lat, lon))
{
return null;
}
double dLat = transformLat(lon - 105.0, lat - 35.0);
double dLon = transformLon(lon - 105.0, lat - 35.0);
double radLat = lat / 180.0 * pi;
double magic = Math.Sin(radLat);
magic = 1 - ee * magic * magic;
double sqrtMagic = Math.Sqrt(magic);
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
double mgLat = lat + dLat;
double mgLon = lon + dLon;
return new Gps(mgLat, mgLon);
}
/**
* * 火星坐标系 (GCJ-02) to 84 * * @param lon * @param lat * @return
* */
public static Gps gcj_To_Gps84(double lat, double lon)
{
Gps gps = transform(lat, lon);
double lontitude = lon * 2 - gps.getWgLon();
double latitude = lat * 2 - gps.getWgLat();
return new Gps(latitude, lontitude);
}
/**
* 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 将 GCJ-02 坐标转换成 BD-09 坐标
*
* @param gg_lat
* @param gg_lon
*/
public static Gps gcj02_To_Bd09(double gg_lat, double gg_lon)
{
double x = gg_lon, y = gg_lat;
double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * pi);
double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * pi);
double bd_lon = z * Math.Cos(theta) + 0.0065;
double bd_lat = z * Math.Sin(theta) + 0.006;
return new Gps(bd_lat, bd_lon);
}
/**
* * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 * * 将 BD-09 坐标转换成GCJ-02 坐标 * * @param
* bd_lat * @param bd_lon * @return
*/
public static Gps bd09_To_Gcj02(double bd_lat, double bd_lon)
{
double x = bd_lon - 0.0065, y = bd_lat - 0.006;
double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * pi);
double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * pi);
double gg_lon = z * Math.Cos(theta);
double gg_lat = z * Math.Sin(theta);
return new Gps(gg_lat, gg_lon);
}
/**
* (BD-09)-->84
* @param bd_lat
* @param bd_lon
* @return
*/
public static Gps bd09_To_Gps84(double bd_lat, double bd_lon)
{
Gps gcj02 = PositionUtil.bd09_To_Gcj02(bd_lat, bd_lon);
Gps map84 = PositionUtil.gcj_To_Gps84(gcj02.getWgLat(),
gcj02.getWgLon());
return map84;
}
/**
* 度分秒转84
* @param lat
* @param lon
* @return
*/
public static Gps dufenmiao_To_Gps84(double lat, double lon)
{
return new Gps(convertdufenmiao(lat), convertdufenmiao(lon));
}
private static double convertdufenmiao(double num)
{
var strlist = num.ToString().Split(‘.‘);
string m_fen = strlist[1].Substring(0, 2);//取分
string m_miao = strlist[1].Substring(2, 2);//取秒
double fen = Convert.ToDouble(m_fen) / 60;
double miao = Convert.ToDouble(int.Parse(m_miao)) / 3600;
double mn = fen + miao;
return Convert.ToDouble((Convert.ToDouble(strlist[0]) + mn).ToString("0.00000000"));
}
/**
* 84转度分秒
* @param lat
* @param lon
* @return
*/
public static Gps gps84_To_Dufenmiao(double lat, double lon)
{
return new Gps(convert84(lat), convert84(lon));
}
private static double convert84(double num)
{
var strlist = num.ToString().Split(‘.‘);
double fentemp = Convert.ToDouble("0." + strlist[1]) * 60;
string fen = Convert.ToDouble(fentemp).ToString();
if (fen.Contains("."))
{
var fenlist = fen.Split(‘.‘);//分
double miaotemp = Convert.ToDouble("0." + fenlist[1]) * 60;
string miao = Convert.ToInt32(miaotemp).ToString();
string temp = strlist[0] + "." + fenlist[0] + Convert.ToInt32(Convert.ToDouble(miao)).ToString();
return Convert.ToDouble(temp);
}
else
{
return Convert.ToDouble(strlist[0] + "." + fen);
}
}
public static Boolean outOfChina(double lat, double lon)
{
if (lon < 72.004 || lon > 137.8347)
return true;
if (lat < 0.8293 || lat > 55.8271)
return true;
return false;
}
public static Gps transform(double lat, double lon)
{
if (outOfChina(lat, lon))
{
return new Gps(lat, lon);
}
double dLat = transformLat(lon - 105.0, lat - 35.0);
double dLon = transformLon(lon - 105.0, lat - 35.0);
double radLat = lat / 180.0 * pi;
double magic = Math.Sin(radLat);
magic = 1 - ee * magic * magic;
double sqrtMagic = Math.Sqrt(magic);
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
double mgLat = lat + dLat;
double mgLon = lon + dLon;
return new Gps(mgLat, mgLon);
}
public static double transformLat(double x, double y)
{
double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y
+ 0.2 * Math.Sqrt(Math.Abs(x));
ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0;
ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0;
return ret;
}
public static double transformLon(double x, double y)
{
double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1
* Math.Sqrt(Math.Abs(x));
ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0;
ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0
* pi)) * 2.0 / 3.0;
return ret;
}
在这里，对于火星坐标，84坐标，还有平面坐标系、大地坐标系，已经高斯平面直角坐标系这类的问题，请自行百度 or google，必经我也不知道很清楚，我就知道他们区别挺大的。
以为到这里就搞定，那么我想着真的是侮辱宝宝的智商了，大爷的，走这个接口的时候，一直让我泪奔呀，10次里面5次的状态吗都是500，也是R了哈士奇了。
经过N次测试，终于确定了一个问题，那就是，原来是点太多了，就是利用gps采集器或者手持设备采集的点太多了，计算的过程中，因为返回的对象是 AreasAndLengths 也就是说返回的不仅仅是面积，还有周长等，那么这时候就需要注意了，如果点太多，就会 over 的，最终得到的结果是180个点以内，是最稳定的。
接下里就开始嘚瑟今天的重要内容。
对于N个点来讲，我们怎么让无限趋近于180这个数字，当然，如果是少于180个点，这个可以直接参与运算，如果是300，330，500这样的点数呢。也就是说我们需要一个算法让GPS的点的总数无线趋近于180。
1、那么最基本的肯定先取商，也就是说用N来除180，商是多少，当然，是取整数的，但是这个时候我们就要去想需要向上取整还是向下取整，OK，我们来回想一下，如果向下取整的话，那么就意味着这时候，商 × 180 小于 N。那么这时候就会出现这种情况，如果有359个点，我们设置的计算常数是180；359 ÷ 180 取整得到 1，那么这时候怎么取呢？
2、比如就这个359个点，我们采取向上取整，也就是每当余数是 0 的时候，我们把这个点加到我们的 计算的点的集合 里面，对于这种情况就已经搞定了。
3、那么，如果有240个点呢，怎么办？向上取整是去到的点数还是 2，这是按照 240 ÷ 2 =120，也就是每两个点取一个点的时候，这时候剩下的就有120个点，然而，我们才取了120个点，所以，我们这时候就要注意一个问题，排序的问题。
4、关于排序，我们现在想想，如果不排序的话，那么，就上一步从剩下的120个点里面取剩下的60个点的时候是不是和原来的顺序就有可能乱呢？比如本来 A点 采集的索引是 3 ，但是计算之后成了我们取了list[0]、list[2]之后，然后再取A点，这时候如果再把这个点加入到要参与计算的这个集合中，是不是顺序就乱了呢？这个点应该就成了 list[120] 了；
5、所以，我们这里需要排序，也就是说对于传输进来的点，参与计算之前，我们就需要对每个点进行排序，然后用每个点的索引参与计算，最终得到的是参与计算的点的索引的集合，从而取到对应的点，然后参与计算；
6、继续刚才的4的问题，如何取剩下的120个点呢，因为之前我们已经取了120个点，而且对于这240个点，这120个点已经参与了计算，那么，如果再次取点的时候，这时候是不能把已经取走的120个点参与计算的，也就是说我们接下来就是从剩下的120个点取出需要的60个点，这个问题就简单了，任然是取余数了。
综上所述，其实取点的过程就是一个循环取值的过程。然后需要注意每个点的排序的问题。
/// <summary>
///
/// </summary>
/// <param name="all">点的数量</param>
/// <param name="count">要取的点</param>
/// <returns></returns>
private static List<int> CalcPointIndex(int all, int count)
{
List<int> beingList = new List<int>();
List<int> selList = new List<int>();
count = count - 10;
//前后各去掉10个值
for (int i = 10; i < all - 10; i++)
{
beingList.Add(i);
}
//从前后各去掉10个值中均匀取5个
int[] nums = { 0, 2, 4, 6, 8 };//取开始的点
int[] nume = { all - 9, all - 7, all - 5, all - 3, all - 1 };//取结束的点
selList = ReturnListInt(nums, selList);
selList = ReturnListInt(nume, selList);
if (beingList.Count < (BASE_NUM - 10))//如果此时点数已经小于，则直接返回
{
foreach (var item in selList)
{
beingList.Add(item);
}
selList = null;
selList = SoftListInt(beingList);
}
else
{
List<int> sel = new List<int>();
List<int> temp = reSel(ref beingList, ref sel, count);
for (int i = 0; i < temp.Count; i++)
{
selList.Add(temp[i]);
}
selList = SoftListInt(selList);
}
return selList;
}
/// <summary>
/// 均匀采样
/// </summary>
/// <param name="beingList"></param>
/// <param name="sel">当前选择的点</param>
/// <param name="count">本次需要取的点的总数</param>
/// <returns></returns>
private static List<int> reSel(ref List<int> beingList, ref List<int> sel, int count)
{
List<int> thistimeAdd = new List<int>();//本次要添加的点
if (count > 0 && count != 1)
{
string t = Math.Ceiling(Convert.ToDouble(beingList.Count / double.Parse(count.ToString()))).ToString();
int sCount = int.Parse(t);//向上取整
int n = 0;//取的是临时变量，存储本次循环添加了多少新的点
for (int i = 0; i < beingList.Count; i++)
{
if (i % sCount == 0 && i != 0)
{
n++;
sel.Add(beingList[i]);
}
else
{
thistimeAdd.Add(beingList[i]);
}
}
beingList = null;
beingList = thistimeAdd;
count = count - n;
if (count > 0)
{
reSel(ref beingList, ref sel, count);
}
}
else
{
string t = Math.Ceiling(Convert.ToDouble(beingList.Count / 2)).ToString();
int sCount = int.Parse(t);//向上取整
sel.Add(beingList[sCount]);
beingList.RemoveAt(sCount);
}
return sel;
}
/// <summary>
/// 循环加载点
/// </summary>
/// <param name="nums"></param>
/// <param name="list"></param>
/// <returns></returns>
private static List<int> ReturnListInt(int[] nums, List<int> list)
{
for (int i = 0; i < nums.Count(); i++)
{
list.Add(nums[i]);
}
return SoftListInt(list);
}
/// <summary>
/// 排序
/// </summary>
/// <param name="sel"></param>
/// <returns></returns>
private static List<int> SoftListInt(List<int> sel)
{
sel.Sort();//Reverse(); 倒序
return sel;
}
这里，需要提醒一下 其中 int nums[] 和 int nume[] 分别取的是开始的点和结束的点，这个点是固定取值的，为了闭合整个图层；
还有对于其中的BASE_NUM这就是常数180
if (beingList.Count < (BASE_NUM - 10))
这句话的目的真的为了防止浪费资源，因为在我看来，如果有1000个点，去掉，880个点和870点区别不是很大。
我觉得，还是再来点实际的吧，必经都是开发项目的不是，你们难道没发现用高德地图的过程中，如果是卫星图，则会出现此区域无卫星图么？其实，人家是提供了解决方案的。直接上代码，这个最简单了，在js里面加入图层就可以了。
//谷歌地图
var googleLayerf = new AMap.TileLayer({
getTileUrl: ‘http://mt{1,2,3,0}.google.cn/vt/lyrs=s&hl=zh-CN&gl=cn&x=[x]&y=[y]&z=[z]&s=Galile‘
});//引用谷歌url
var roadNetLayerf = new AMap.TileLayer.RoadNet(); //
var map = new AMap.Map("container", {
resizeEnable: true,
zoom: 12,
layers: [googleLayerf, roadNetLayerf],
center: [longit, latitu]//地图中心点
});
好了，这事情就这么搞定了。虽然装逼挺累的，但是这玩意想出来还有得费点劲的，特别中间的排序和索引的问题。
以上是关于GIS面积计算问题的主要内容，如果未能解决你的问题，请参考以下文章