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面积计算问题的主要内容,如果未能解决你的问题,请参考以下文章

GIS面积计算问题

GIS计算多个地块面积

GIS面积计算

如何用GIS算指定区域的面积和体积

gis怎么算两个面相交的占比

计算 shapefile 的面积