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]//地图中心点
    });

  好了,这事情就这么搞定了。虽然装逼挺累的,但是这玩意想出来还有得费点劲的,特别中间的排序和索引的问题。

时间: 2024-12-28 22:00:57

GIS面积计算问题的相关文章

网上一些过来人对GIS的看法

GIS学生编程不如计算机的,搞测绘不如学测绘的,搞地理不如学地理的,我现在开始学习编程再怎么学也赶不上学计算机的怎么办? 这个问题简单,如果现在去睡觉能赶上他们的话,完全可以去睡觉嘛. 1.家人问你学什么专业,别说地理信息系统,他们听不懂,说你是学计算机的. 你就是学计算机的,做GIS,就要有良好的计算机基础,面向对象,数据结构.数据库,计算机网络,操作系统,软件工程,设计模式等,这中间涉及的东西太复杂,其实又很简单,做计算机,要的就是基础,很多东西不是看一遍会做就行,C++的书看过不下十遍,T

@最新CAX/EDA/CFD/GIS/光学/化工/液压软件资源网

最新CAX/EDA/CFD/GIS/光学/化工/液压软件资源网 阳光软件园 所有软件资料都随时更新,急需软件可以去看看,基本上能找到你想要的! http://zhangqg.51.net http://cax2one.f3322.net e-mail: [email protected];[email protected];[email protected] 将以上任意链接连接起来输入IE 窗口即可进入网站 下面是一部分软件,更多软件在我们的软件列表,如需要请到列表中去找! ACTRAN v14

Web前台直接加载GIS格式数据分析

本文以Flex直接加载Shp.DWG和MDB为例. 首先看一份现估测数据: 1)  加载Shp文件,目前直接由前台Flex代码完成: 图1 在ArcCatalog里面的Shp文件 图2 直接在前台加载后的Shp文件 结果显示: Shp文件 大小 加载时间 Shp1 50kb 约3s Shp2 750kb 约10s 分析:未用后台开发,直接使用前台Flex对SHP开放数据加载,省去通讯时间,速度快捷,速度与客户端配置成正比. 说明:直接加载使用了LibertyGIS.swc组件. 2)  加载Dw

开源GIS简介

原文 开源GIS C++开源GIS中间件类库: GDAL(栅格)/OGR(矢量)提供了类型丰富的读写支持 GEOS(Geometry Engine Open Source)是基于C++的空间拓扑分析实现类库,遵循LGPL协议发布.GEOS类库提供了丰富的空间拓扑操作函数,用以判断几何对象间的相互关系,以及空间分析操作之后形成新的几何对象.点.线.面要素的两两相互关系,包括相合.分离.相交.重合.包含.相邻等不同位置关系,都可以通过GEOS类库中提供的函数进行分析和判断.并且GEOS类库提供了缓冲

利用ArcGIS Engine、VS .NET和Windows控件开发GIS应用

原文:利用ArcGIS Engine.VS .NET和Windows控件开发GIS应用 此过程说明适合那些使用.NET建立和部署应用的开发者,它描述了使用ArcGIS控件建立和部署应用的方法和步骤. 你可以在下面的目录下找到相应的样例程序: <安装目录>\DeveloperKit\Samples\Developer_Guide_Scenarios\ ArcGIS_Engine\Building_an_ArcGIS_Control_Application\Map_Viewer 注:ArcGIS样

[转]基于C#的开源GIS项目介绍之SharpMap篇

我是一个刚毕业的GIS本科毕业生,目前在杭州从事GIS软件应用开发.在项目开发中总感觉自己的编程水平还不够,于是想找些开源GIS小项目来研究研究,借以提高自己的编程能力和项目开发能力.在网上搜了一下“GIS开源”发现还不少,下面是一个介绍GIS开源项目的链接: http://www.yuanma.org/data/2008/0526/article_3048.htm 里面介绍了基于各种编程语言的GIS开源项目,并列出了各自的特点和官网链接. 由于在学校时候学的一直都是C#和Visual Stud

Web GIS离线解决方案

1.背景 在离线环境下(局域网中)的GIS系统中如何使用地图?这里的地图主要指的是地图底图,有了底图切片数据,我们就可以看到地图,在上面加上自己的业务数据图层,进行相关操作. 要在离线环境下看到GIS地图,就要有底图切片数据,地图的底图切片数据在一定时间内是不会变化的,可以使用一些地图下载器下载地图切片,如这个地图下载器. 在CS系统中可以基于GMap.Net来做,参考<百度谷歌离线地图解决方案>. 下面介绍下Web系统如何使用GIS切片数据,开发web GIS系统. 2.使用GeoWebCa

GIS开发离线地图应用-初识gis

最新公司需要做一个基于gis地图的应用系统,由于之前公司项目中的电子地图模块都是我开发的,所以这个新系统也自然让我先去了解如何开发,可以说做个简单的调研. 和之前的项目中开发的电子地图模块不同,这次是开发gis地图,是要显示真实的地理位置,能有gps定位功能的.而之前开发过的电子地图功能,都只是基于svg的矢量可配置地图(之前采用batik开发过C/S版,用raphael开发过B/S版,都在项目中正常使用). 下面描述下我开始开发前做的准备和了解工作,希望对首次接触并想要开发gis离线地图应用的

使用OpenGL绘制 shapefile文件 完成最基本的gis操作

主要内容概述 1.解析shapefile(.shp)文件.‘ 2.将经纬度数据按照墨卡托投进行投影(调用proj.4库)完成. 3.将数据用OpenGL的方式进行绘制. 上述3方面只是完成初步的绘制,对于要完成一个复杂的地理信息系统还有很大的差距, 下面介绍我设计的简单的地理信息框架(用于交流,进步).先上个图,根据图来说更加易懂. 其中地图框架中包含多个物理地理图层,是真实存在的图层,所有在该图层下的数据都会被绘制到图层上(即一张图片). 物理图层中包含了多个逻辑图层(Layer),是为了方便