百度坐标 转 WGS84坐标 并保存到SHP图层中

步骤1、坐标转换代码

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Linq;

using System.Text;

using System.Windows.Forms;

namespace BDTOWGS84

{

    public partial class Form1 : Form

    {

        public Form1()

        {

            InitializeComponent();

        }

        /*

         * BD-09：百度坐标系(百度地图)

         * GCJ-02：火星坐标系（谷歌中国地图、高德地图）

         * WGS84：地球坐标系（国际通用坐标系，谷歌地图）

         */

        public double x_PI = 3.14159265358979324 * 3000.0 / 180.0;

        public double PI = 3.1415926535897932384626;

        public double a = 6378245.0;

        public double ee = 0.00669342162296594323;

        private void Form1_Load(object sender, EventArgs e)

        {

            //从shp中读取

            double bd_lon = 115.962804;//经度

            double bd_lat = 37.3786;//纬度

            //将wgs存入shp

            double[] wgs = bd09togcj02(bd_lon, bd_lat);

            MessageBox.Show(wgs[0].ToString()+"    "+wgs[1].ToString());

        }

        //百度坐标系转火星坐标系

        private double[] bd09togcj02(double bd_lon, double bd_lat)

        {

            double x = bd_lon - 0.0065;

            double y = bd_lat - 0.006;

            double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * x_PI);

            double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * x_PI);

            double gcj_lon = z * Math.Cos(theta);

            double gcj_lat = z * Math.Sin(theta);

            double[] gcj = { gcj_lon, gcj_lat };//火星坐标系值

            //火星坐标系转wgs84

            double[] wgs = gcj02towgs84(gcj[0], gcj[1]);

            return wgs;

        }

        //火星坐标系转wgs84

        private double[] gcj02towgs84(double gcj_lon, double gcj_lat)

        {

           if (out_of_china(gcj_lon, gcj_lat)) 

           {

               MessageBox.Show("不在国内范围！");

               double[] back={gcj_lon,gcj_lat};

               return back;

           } 

           else 

           {

                var dlon = transformlon(gcj_lon - 105.0, gcj_lat - 35.0);

                var dlat = transformlat(gcj_lon - 105.0, gcj_lat - 35.0);

                var radlat = gcj_lat / 180.0 * PI;

                var magic = Math.Sin(radlat);

                magic = 1 - ee * magic * magic;

                var sqrtmagic = Math.Sqrt(magic);

                dlon = (dlon * 180.0) / (a / sqrtmagic * Math.Cos(radlat) * PI);

                dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI);

                double mglon = gcj_lon + dlon;

                double mglat = gcj_lat + dlat;

                double wgs_lon = gcj_lon * 2 - mglon;

                double wgs_lat = gcj_lat * 2 - mglat;

                double[] wgs = { wgs_lon,wgs_lat };//wgs84坐标系值

                return wgs;

            }

        }

        //火星坐标系转百度坐标系

        private double[] gcj02tobd09(double gcj_lon, double gcj_lat)

        {

            double z = Math.Sqrt(gcj_lon * gcj_lon + gcj_lat * gcj_lat) + 0.00002 * Math.Sin(gcj_lat * x_PI);

            double theta = Math.Atan2(gcj_lat, gcj_lon) + 0.000003 * Math.Cos(gcj_lon * x_PI);

            double bd_lon = z * Math.Cos(theta) + 0.0065;

            double bd_lat = z * Math.Sin(theta) + 0.006;

            double[] bd={bd_lon,bd_lat};

            return bd;

        }

        //wgs84转火星坐标系

        private double[] wgs84togcj02(double wgs_lon, double wgs_lat)

        {

            if (out_of_china(wgs_lon, wgs_lat)) 

            {

                MessageBox.Show("不在国内范围！");

                double[] back={wgs_lon,wgs_lat};

                return back;

            } 

            else 

            {

                double dwgs_lon = transformlon(wgs_lon - 105.0, wgs_lat - 35.0);

                double dwgs_lat = transformlat(wgs_lon - 105.0, wgs_lat - 35.0);

                double radwgs_lat = wgs_lat / 180.0 * PI;

                double magic = Math.Sin(radwgs_lat);

                magic = 1 - ee * magic * magic;

                double sqrtmagic = Math.Sqrt(magic);

                dwgs_lon = (dwgs_lon * 180.0) / (a / sqrtmagic * Math.Cos(radwgs_lat) * PI);

                dwgs_lat = (dwgs_lat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI);

                double gcj_lon = wgs_lon + dwgs_lon;

                double gcj_lat = wgs_lat + dwgs_lat;

                double[] gcj = {gcj_lon, gcj_lat};

                return gcj;

            }

        }

        private double transformlon(double lon, double lat)

        {

            var ret = 300.0 + lon + 2.0 * lat + 0.1 * lon * lon + 0.1 * lon * lat + 0.1 * Math.Sqrt(Math.Abs(lon));

            ret += (20.0 * Math.Sin(6.0 * lon * PI) + 20.0 * Math.Sin(2.0 * lon * PI)) * 2.0 / 3.0;

            ret += (20.0 * Math.Sin(lon * PI) + 40.0 * Math.Sin(lon / 3.0 * PI)) * 2.0 / 3.0;

            ret += (150.0 * Math.Sin(lon / 12.0 * PI) + 300.0 * Math.Sin(lon / 30.0 * PI)) * 2.0 / 3.0;

            return ret;

        }

        private double transformlat(double lon, double lat)

        {

            var ret = -100.0 + 2.0 * lon + 3.0 * lat + 0.2 * lat * lat + 0.1 * lon * lat + 0.2 * Math.Sqrt(Math.Abs(lon));

            ret += (20.0 * Math.Sin(6.0 * lon * PI) + 20.0 * Math.Sin(2.0 * lon * PI)) * 2.0 / 3.0;

            ret += (20.0 * Math.Sin(lat * PI) + 40.0 * Math.Sin(lat / 3.0 * PI)) * 2.0 / 3.0;

            ret += (160.0 * Math.Sin(lat / 12.0 * PI) + 320 * Math.Sin(lat * PI / 30.0)) * 2.0 / 3.0;

            return ret;

        }

        //判断是否在国内，不在国内则不做偏移

        private Boolean out_of_china(double lon, double lat)

        {

            return (lon < 72.004 || lon > 137.8347) || ((lat < 0.8293 || lat > 55.8271) || false);

        }

    }

}

步骤2、转换后的XY坐标生成新的shp图层（ArcToolbox下的创建XY事件图层）

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using ESRI.ArcGIS.Controls;

using ESRI.ArcGIS.Carto;

using ESRI.ArcGIS.esriSystem;

using ESRI.ArcGIS.Geometry;

using System.Windows.Forms;

using ESRI.ArcGIS.Display;

using ESRI.ArcGIS.Geodatabase;

using ESRI.ArcGIS.DataManagementTools;

using ESRI.ArcGIS.Geoprocessor;

using ESRI.ArcGIS.Geoprocessing;

using ESRI.ArcGIS.DataSourcesGDB;

namespace BDTOWGS84

{

    public class ThematicMap

    {

        Form1 f1 = new Form1();

        private IMapControl3 m_mapControl = null;//为axMapControl1

        ILayer pLayer = null;

        public double[] itemValue_X;//各X值,MainForm会调用，然后传到FormPie中去计算落在各个区间的个数

        public double[] itemValue_Y;

        public ThematicMap(IMapControl3 mapControl, AxTOCControl axTOCControl1, ILayer ppLayer, string strName_X, string strName_Y,int numall)

         {

             m_mapControl = mapControl;

             pLayer = ppLayer;//获取图层 ，并把它设置成IGeoFeatureLayer的实例 

             IFeatureLayer pFeatureLayer = pLayer as IFeatureLayer;

             IGeoFeatureLayer pGeoFeatureLayer = pLayer as IGeoFeatureLayer;

             IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;获取图层上的feature

             // 通过查找features的所有字段的值

             ITable pTable;

             int fieldNumber_X=0;

             int fieldNumber_Y = 0;

             pTable = pGeoFeatureLayer as ITable;

             fieldNumber_X = pTable.FindField(strName_X);

             fieldNumber_Y = pTable.FindField(strName_Y);

             if (fieldNumber_X == -1)

             {

                 MessageBox.Show(ppLayer.Name + "图层中" + strName_X + "字段未被找到，请检查SHP文件！", "Message", MessageBoxButtons.OK, MessageBoxIcon.Information);

                 return;

             }

             if (fieldNumber_Y == -1)

             {

                 MessageBox.Show(ppLayer.Name + "图层中" + strName_Y + "字段未被找到，请检查SHP文件！", "Message", MessageBoxButtons.OK, MessageBoxIcon.Information);

                 return;

             }

             const int numFields = 1;// 设置bars的个数

             int[] fieldIndecies_X = new int[1];

             int[] fieldIndecies_Y = new int[1];

             double fieldValue_X=0;//X字段值

             double fieldValue_Y = 0;//Y字段值

             fieldIndecies_X[0] = pTable.FindField(strName_X);

             fieldIndecies_Y[0] = pTable.FindField(strName_Y);

             int n = pFeatureClass.FeatureCount(null);//列数，关键

             itemValue_X = new double[n];

             itemValue_Y = new double[n];

             for (int i = 0; i < numFields; i++)

             {

                 IFeatureCursor pFeatureCursor = pFeatureClass.Search(null, false);

                 for (int j = 0; j < n; j++)

                 {

                     IFeature pFeature = pFeatureCursor.NextFeature();

                     fieldValue_X = Convert.ToDouble(pFeature.get_Value(fieldIndecies_X[i]));

                     fieldValue_Y = Convert.ToDouble(pFeature.get_Value(fieldIndecies_Y[i]));

                     double[] wgs = f1.bd09togcj02(fieldValue_X, fieldValue_Y);//坐标转换

                     if (wgs[0] == 0)

                         continue;

                     itemValue_X[j] = wgs[0];

                     itemValue_Y[j] = wgs[1];

                     pFeature.set_Value(fieldIndecies_X[i],wgs[0]);//更改shp属性字段的值

                     pFeature.set_Value(fieldIndecies_Y[i], wgs[1]);

                     pFeature.Store();

                     if (j == n - 1)

                     {

                         //创建XY事件图层

                         MakeXYEventLayer MxyLayer = new MakeXYEventLayer();

                         MxyLayer.table = pLayer;

                         MxyLayer.in_x_field = "X";

                         MxyLayer.in_y_field = "Y";

                         MxyLayer.out_layer = pLayer.Name + "_AfterTransform";

                         MxyLayer.spatial_reference = "4326";//不能写GCS_WGS_1984,要写4326

                         Geoprocessor GP = new Geoprocessor();

                         GP.Execute(MxyLayer, null);//执行MxyLayer事件

                         //删除gdb下存在的shp文件,首先用Catalog在指定目录下新建gdb文件

                         IWorkspaceFactory worFact = new FileGDBWorkspaceFactory();

                         IWorkspace workspace = worFact.OpenFromFile(Application.StartupPath + "\\ BDTOWGS84 .gdb", 0);

                         IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;

                         IFeatureWorkspaceManage featureWorkspaceMange = (IFeatureWorkspaceManage)featureWorkspace;

                         IEnumDatasetName enumDatasetName = workspace.get_DatasetNames(esriDatasetType.esriDTFeatureClass);

                         IDatasetName datasetName = null;

                         while ((datasetName = enumDatasetName.Next()) != null)

                         {

                             if (datasetName.Name.Equals(MxyLayer.out_layer))

                             {

                                 featureWorkspaceMange.DeleteByName(datasetName);//删除指定要素类

                                 break;

                             }

                         }

                         //在gdb下新建shp文件，首先必须要删除原有重名文件，否则会报错

                         ESRI.ArcGIS.DataManagementTools.CopyFeatures cf = new ESRI.ArcGIS.DataManagementTools.CopyFeatures();

                         cf.in_features = MxyLayer.out_layer;

                         cf.out_feature_class = Application.StartupPath + "\\ BDTOWGS84 .gdb\\" + MxyLayer.out_layer;

                         //cf.out_feature_class = "D:\\test3.gdb\\" + MxyLayer.out_layer;

                         GP.Execute(cf, null);

                         MessageBox.Show("图层" + pLayer.Name + "转换完毕!");

                     }        

                 }

             }

             if (numall!=0)

                MessageBox.Show("全部图层转换完毕！！！", "Message", MessageBoxButtons.OK, MessageBoxIcon.Information);

         }

    }

}