西安80北京54,2000和WGS84互转C#程序

24 篇文章 3 订阅

调用方法

static double dmslon, dmslat;

ICoordinate ic = CoordinateFactory.CreateCoordinate(Spheroid.Xian80);

ic.GaussPrjInvCalculate(36543406.9394, 3357113.48175, out dmslon, out dmslat);

输出

108.4513888325,30.3329762486

注意区分3度带和6度带

六度带中央经线经度的计算:当地中央经线经度=6°×当地带号-3°

三度带中央经线经度的计算:中央经线经度=3°×当地带号

这块注意在国内范围的话:

在中华人民共和国陆地范围内,坐标(Y坐标,8位数,前两位是带号)带号小于等于23的肯定是6度带,大于等于24的肯定是3度带。

longitude0 = ProjNo >= 24 ? ProjNo * 3 : ProjNo * 6 - 3;

类文件

using System;
using System.Collections.Generic;
using System.Text;

namespace GIS.Common
{
    public abstract class GaussPrjBase: ICoordinate
    {        
        protected double _a;
        protected double _f;
        private int _zoneWide = 3;
        protected readonly double PI = 3.14159265353846;

        #region ICoordinate Members
        public int ZoneWide
        {
            get { return _zoneWide; }
            set { _zoneWide = value; } 
        }

        public void DD2DMS(double DecimalDegree, out int Degree, out int Minute, out double Second)
        {
            Degree = (int)DecimalDegree;
            Minute = (int)((DecimalDegree - Degree) * 60.0);
            Second = Math.Round((DecimalDegree * 60 - Degree * 60 - Minute) * 60 * 100) / 100.0;
        
        }

        public double DistanceOfTwoPoints(double lng1, double lat1, double lng2, double lat2, Spheroid gs)
        {
            double radLat1 = Rad(lat1);
            double radLat2 = Rad(lat2);
            double a = radLat1 - radLat2;
            double b = Rad(lng1) - Rad(lng2);
            double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) +
            Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2)));
            s = s * (gs == Spheroid.WGS84 ? 6378137.0 : (gs == Spheroid.Xian80 ? 6378140.0 : 6378245.0));
            s = Math.Round(s * 10000) / 10000;
            return s;
        }      

        private  double Rad(double d)
        {
            return d * Math.PI / 180.0;
        }
        public void DMS2DD(int Degree, int Minute, double Second, out double DecimalDegree)
        {
            DecimalDegree = Degree + Minute / 60.0 + Second / 60.0 / 60.0;
        }

        
        public void GaussPrjCalculate(double longitude, double latitude, out double X, out double Y)
        {           
            int ProjNo = 0;
            double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
            double e2, ee, NN, T, C, A, M, iPI;
            iPI = 0.0174532925199433; //3.1415926535898/180.0; 
            //ZoneWide = 6; //6度带宽 
            ProjNo = (int)(longitude / ZoneWide);
            longitude0 = ZoneWide == 3 ? ProjNo * 3 : ProjNo * 6 - 3;
            longitude0 = longitude0 * iPI;
            latitude0 = 0;
            longitude1 = longitude * iPI; //经度转换为弧度
            latitude1 = latitude * iPI; //纬度转换为弧度
            e2 = 2 * _f - _f * _f;
            ee = e2 * (1.0 - e2);
            NN = _a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1));
            T = Math.Tan(latitude1) * Math.Tan(latitude1);
            C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
            A = (longitude1 - longitude0) * Math.Cos(latitude1);
            M = _a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1));
            xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);
            yval = M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);
            X0 = 1000000L * (ProjNo + 1) + 500000L;
            Y0 = 0;
            X = Math.Round((xval + X0) * 10000000000) / 10000000000.0;
            Y = Math.Round((yval + Y0) * 10000000000) / 10000000000.0;
        }

        public void GaussPrjInvCalculate(double X, double Y, out double longitude, out double latitude)
        {
            int ProjNo;
            double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval;
            double e1, e2, ee, NN, T, C, M, D, R, u, fai, iPI;
            iPI = 0.0174532925199433; //3.1415926535898/180.0; 
            ProjNo = (int)(X / 1000000L); //查找带号
            //longitude0 = ProjNo * 3;
            longitude0 = ZoneWide == 3 ? ProjNo * 3 : ProjNo * 6 - 3;
            
            longitude0 = longitude0 * iPI; //中央经线
            X0 = ProjNo * 1000000L + 500000L;
            Y0 = 0;
            xval = X - X0;
            yval = Y - Y0; //带内大地坐标
            e2 = 2 * _f - _f * _f;
            e1 = (1.0 - Math.Sqrt(1 - e2)) / (1.0 + Math.Sqrt(1 - e2));
            ee = e2 / (1 - e2);
            M = yval;
            u = M / (_a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
            fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(4 * u)
            + (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.Sin(8 * u);
            C = ee * Math.Cos(fai) * Math.Cos(fai);
            T = Math.Tan(fai) * Math.Tan(fai);
            NN = _a / Math.Sqrt(1.0 - e2 * Math.Sin(fai) * Math.Sin(fai));
            R = _a * (1 - e2) / Math.Sqrt((1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)));
            D = xval / NN;
            //计算经度(Longitude) 纬度(Latitude)
            longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / Math.Cos(fai);
            latitude1 = fai - (NN * Math.Tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
            //转换为度 DD
            longitude = Math.Round((longitude1 / iPI) * 10000000000) / 10000000000.0;
            latitude = Math.Round((latitude1 / iPI) * 10000000000) / 10000000000.0; 
        }
        #endregion        
    }

    /// <summary>
    /// 高斯投影选用的参考椭球
    /// </summary>
    public enum Spheroid
    {
        Xian80,
        Beijing54,
        WGS84,
    }
}
 public class CoordinateFactory
    {
        public static ICoordinate CreateCoordinate(Spheroid s)
        {
            ICoordinate coordinate;
            switch (s)
            {
                case Spheroid.Beijing54:
                    coordinate = new Beijing54();
                    break;
                case Spheroid.WGS84:
                    coordinate = new WGS84();
                    break;
                case Spheroid.Xian80:
                    coordinate = new Xian80();
                    break;
                default:
                    coordinate = null;
                    break;
            }
            return coordinate;
        }

        public static ICoordinate CreateCoordinate()
        {
            return new WGS84();
        }
    }
public interface ICoordinate
    {        
        /// <summary>
        /// 带宽
        /// </summary>
        int ZoneWide { get;set;}

        /// <summary>
        /// 十进制双精度角度转换成度分秒角度格式
        /// </summary>
        /// <param name="Decimal Degree">度,十进制型双精度</param>
        /// <param name="Degree">度,整型</param>
        /// <param name="Minute">分,整型</param>
        /// <param name="Second">秒,双精度型</param>
        void DD2DMS(double DecimalDegree, out int Degree, out int Minute, out double Second);

        /// <summary>
        /// 求两点之间的距离(根据经纬度)
        /// </summary>
        /// <param name="lng1">经度1</param>
        /// <param name="lat1">纬度1</param>
        /// <param name="lng2">经度2</param>
        /// <param name="lat2">纬度2</param>
        /// <param name="gs">高斯投影中所选用的参考椭球</param>
        /// <returns>两点间距离(单位:meters)</returns>
        double DistanceOfTwoPoints(double lng1, double lat1, double lng2, double lat2, Spheroid gs);

        
        
        

        /// <summary>
        /// 求两点之间的距离(大地坐标)
        /// </summary>
        /// <param name="x1"></param>
        /// <param name="y1"></param>
        /// <param name="x2"></param>
        /// <param name="y2"></param>
        /// <returns>单位为meters</returns>
        //double DistanceOfTwoPoints(double x1, double y1, double x2, double y2);

        /// <summary>
        /// 度分秒角度格式转换成十进制度双精度角度格式
        /// </summary>
        /// <param name="Degree">度,整型</param>
        /// <param name="Minute">分,整型</param>
        /// <param name="Second">秒,双精度型</param>
        /// <param name="Decimal Degree">度,十进制型双精度</param>   
        void DMS2DD(int Degree, int Minute, double Second, out double DecimalDegree);

        /// <summary>
        /// 高期投影正算
        /// 由经纬度(单位:Decimal Degree)正算到大地坐标(单位:Metre,含带号)
        /// </summary>
        /// <param name="longitude">经度</param>
        /// <param name="latitude">纬度</param>
        void GaussPrjCalculate(double longitude, double latitude,out double X,out double Y);

        /// <summary>
        /// 高斯投影反算
        /// 大地坐标(单位:Metre,含带号)反算到经纬度坐标(单位,Decimal Degree)
        /// </summary>
        /// <param name="X">大地坐标X值</param>
        /// <param name="Y">大地坐标Y值</param>
        void GaussPrjInvCalculate(double X, double Y,out double longitude,out double latitude);

    }
/// <summary>
    /// Krasovsky参考椭球 参数(北京54坐标系采用此参考椭球)
    /// </summary>
    public class Beijing54 : GaussPrjBase
    {

        public Beijing54()
        {
            _a = 6378245.0;
            _f = 1.0 / 298.3;
        }
    }
/// <summary>
    /// WGS84参考椭球 参数(全球GPS监测网采用此参考椭球)
    /// </summary>
    public class WGS84 : GaussPrjBase
    {
        public WGS84()
        {
            _a = 6378137.0;
            _f = 1.0 / 298.257223563;

        }
    }
/// <summary>
    /// IUGG1975参考椭球 参数(北京54坐标系采用此参考椭球)
    /// </summary>
    public class Xian80 : GaussPrjBase
    {
        public Xian80()
        {
            _a = 6378140.0;
            _f = 1.0 / 298.257;
        }
    }

 

  • 2
    点赞
  • 13
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值