精华内容
下载资源
问答
  • 包括了北京54坐标系、北京80坐标系、WGS-84坐标系的椭球参数
  • 求x,y的值(参考值x=3380330.773,y=320089.969)。
  • 常用坐标系椭球参数整理

    千次阅读 2016-01-06 14:30:00
    我们在进行坐标转换的过程中经常会用到几个常用坐标系的椭球参数(例如高斯正反算、相同坐标系下大地坐标与空间直角坐标的转换等,后续会整理这些计算的基本公式和代码),下面列出几个我们常用的椭球体的基本参数:...

        我们在进行坐标转换的过程中经常会用到几个常用坐标系的椭球参数(例如高斯正反算、相同坐标系下大地坐标与空间直角坐标的转换等,后续会整理这些计算的基本公式和代码),下面列出几个我们常用的椭球体的基本参数:

     

    克拉索夫斯基椭球体

    1975国际椭球

    WGS84椭球体

    国家2000坐标系椭球

    长半轴(a)

    6378245

    6378140

    6378137

    6378137

    短半轴(b)

    6356863.0187730473

    6356755.288157528

    6356752.3142451795

    6356752.3141403558

        从左到右分别可以对应我们常用的北京54坐标系、西安80坐标系、WGS84坐标系以及CGCS2000坐标系。

        我们表示椭球习惯上是用长半轴a以及扁率表示,现在有了长半轴和短半轴参数信息,习惯的表示参数我们可以下面根据公式求得。

    image

    (Word2013的公式无法插入到WindowsLiveWriter当中,在此通过图片进行展示)

     

    展开全文
  • 该坐标转换系统根据目前我国常用测量坐标转换的需求以及当前主流的坐标转换方法和数学模型,利用C#语言编写了测量坐标转换系统的程序,较好地实现了同一参考椭球下和不同参考椭球之间不同坐标形式的相互转换,即高斯...
  • 求B,L的值。 参考值B=109800.0034″(30°30′),L= 411600.0004″(114°20′)
  • 测量椭球参数换算及CAD VBA编程.pdf
  • 2000国家大地坐标系椭球参数与GRS80和WGS84的比较,让你很明白
  • icesat/glas与SRTM参考椭球和基准对比: 资料来源:Evaluation of the Vertical Accuracy of Open Global DEMs over Steep Terrain Regions Using ICESat Data: A Case Study over Hunan Province, China. icesat ...

    icesat/glas与SRTM参考椭球和基准对比:
    在这里插入图片描述
    资料来源:Evaluation of the Vertical Accuracy of Open Global DEMs over Steep Terrain Regions Using ICESat Data: A Case Study over Hunan Province, China.

    icesat 参考椭球:TOPEX/Poseidon椭球体
    水准面:TOPEX/Poseidon
    垂直精度:±14cm

    SRTM参考椭球:WGS84椭球体
    水准面:EGM96模型
    SRTM3绝对高程精度:±16m

    两者基准转换:
    1)椭球体转换
    WGS84椭球体和TOPEX/Poseidon 椭球体参数下的高程数据差异在70~71cm 之间。
    H_WGS84 = H_TP - 0.7

    2)icesat高程转换
    获取的icesat的i_elevation参数是相对TP椭球体的高程,需要转换成正高,即相对geoid的高程。
    icesat数据中有一个参数i_gdHt是大地水准面差距,即geoid height.、
    地面高程(正高)H是指从一地面点沿过此点的重力线到大地水准面的距离(即地面点到geoid的距离)。

    计算方程为:H_TP = i_elevation - i_gdHt

    可以参考该资料理解:https://blog.csdn.net/qq_32649321/article/details/115631775?spm=1001.2014.3001.5501

    因此,H_WGS84 = i_elevation - i_gdHt - 0.7。

    参考文献:
    Evaluation of the Vertical Accuracy of Open Global DEMs over Steep Terrain Regions Using ICESat Data:A Case Study over Hunan Province, China.

    基于ICESat/GLAS数据的中国典型区域SRTM 与 ASTER GDEM 高程精度评价.

    展开全文
  • 测量上和制图上采用一个非常接近大地体的旋转椭球体作为地球的参考形状和大小,称为地球椭球体,简称椭球体。 由于地球是一个两极压扁的椭球体,斯托克斯在理论上证明:若地球表面重力已知,可以推导出地球表面...

    参考椭球体亦称“参考扁球体”或“参考椭圆体”。椭圆绕其短轴旋转所成的形体,是形状大小一定,且经过定位定向地球椭球体。是与某个区域如一个国家大地水准面最为密合椭球面

     

    测量上和制图上采用一个非常接近大地体的旋转椭球体作为地球的参考形状和大小,称为地球椭球体,简称椭球体

    由于地球是一个两极压扁的椭球体,斯托克斯在理论上证明:若地球表面重力已知,可以推导出地球表面理论公式,即与地球表面最接近的重力等位面方程——参考椭球面。也即一个国家地区处理测量成果采用一种与地球大小形状最接近具有一定参数地球椭球。通常所说地球的形状和大小,实际上就是以参考椭球体半长径半短径扁率来表示。

     

    大地测量学中, 参考椭球体是一个数学定义地球表面,它近似大地水准面。 由于其相对简单,参考椭球是大地控制网计算和显示点坐标(如纬度,经度和海拔)的首选的地球表面的几何模型。通常所说地球的形状和大小,实际上就是以参考椭球体的长半轴、短半轴和扁率来表示的。参考椭球面是测量计算的基准面,法线是测量计算的基准线。我国的大地原点,即椭球定位做最佳拟合的参考点位于陕西省泾阳县永乐镇北洪流村(也有教材称石际寺村)。

     

    椭球性质

    椭球是一种二次曲面,是椭圆在三维空间的推广。椭球在xyz-笛卡儿坐标系中的方程是:

    其中a和b是赤道半径(沿着x和y轴),c是极半径(沿着z轴)。这三个数都是固定的正实数,决定了椭球的形状。

    如果三个半径都是相等的,那么就是一个球;如果有两个半径是相等的,则是一个类球面。点(a,0,0)、(0,b,0)和(0,0,c)都在曲面上。从原点到这三个点的线段,称为椭球的半主轴。它们与椭圆的半长轴和半短轴相对应。

     

    记长轴半径,短轴半径  。常用的地球参考椭球在直角坐标系Oxyz中可表示为:

    长短轴半径及扁率{\displaystyle f}之间有如下关系:

    有时还会用到偏心率:

    第一偏心率: 

    第二偏心率: 

     

    地球参考椭球

    1975年国际大地测量与地球物理联合会推荐的数据为:半长径6378140米,半短径6356755米,扁率1∶298.257。1980年国际大地测量与地球物理联合会推荐数据为:长半轴a=6378137,短半轴b=6356752,扁率α=(a-b)/a=1:298.257。而中国在1978年推测的数据为:a=6378143,b=6356758,α=1:298.255。

     

    下图列出了一些最常见的参考椭球:

    最常用的参考椭球,是美国国防部制图局(DMA)在1984年构建的WGS84。

    大陆地区在1954年前曾采用International 1924参考椭球,之后较长一段时间内采用基于克拉索夫斯基(Krasovsky)1940的1954年北京坐标系。1980年开始使用1975年国际大地测量与地球物理联合会第16届大会推荐的参考椭球。

     

    参考椭球作用和意义

    参考椭球的主要作用就是作为定义经度纬度高程的基础。

    参考椭球体与大地水准面具有确定位置关系。它是一个在局部范围内(一个国家或地区)与大地水准面在各个方面都最接近理想椭球体。它的表面称为参考椭球面,是测量计算的基础面。根据参考椭球面,可以建立经纬度系统,以致地球上任何一点的位置可以用经纬度来描述。经度线即与地轴重合的平面与参考椭球面之交线;纬度线即垂直地轴的平面与参考椭球面之交线。

     

    ---------------The End---------------

    微信关注  奔跑的GISer  获取更多GIS学习资源

     

    展开全文
  • (1)椭球参数 (2)大地坐标、平面坐标、直角坐标 3.Theory (1)椭球长半轴和短半轴计算其它参数 (2)高斯投影正反算 4.Environment and Configuration 5.Conclusion 1.Intro 搞了一段时间的WebGIS和接口...

    目录

    1.Intro

    2.Details

    (1)椭球参数

    (2)大地坐标、平面坐标

    3.Theory

    (1)利用椭球长半轴和短半轴计算其它参数

    (2)高斯投影正反算

    ①高斯投影正算:大地坐标/经纬度坐标(L, B) ==> 平面坐标(X, Y)。

    ②高斯投影反算:平面坐标(X, Y) ==> 大地坐标/经纬度坐标(L, B)。

    4.Environment and Configuration

    5.Conclusion


    1.Intro

    搞了一段时间的WebGIS和接口,突然被领导喊来研究松散型切片(PNG)的切片方法,看到某人写的源码后傻眼了,东拼西凑导致最后瓦片拼接成图的时候,会造成偏差但又找不到原因,美名其曰研究,其实就是收拾烂摊子,水仙不开花装蒜呢。好在人帅功夫深(事实),发现是平面坐标转大地坐标/经纬度坐标(反算)的时候,椭球参数计算误差造成的精度丢失,最后大地坐标转瓦片坐标的时候,自然存在了偏移。

    本着能Ctrl+CV绝不自己写的原则,提取了那一滩狗屎代码里的精华,又参考了网上的一些内容(就是copy),单独生成了一个能设置椭球参数和进行高斯正反算的模块。而关于椭球参数(这里指七参),最好找官方或者论文资料的数据,比较可信,有些文章代码里的七参精度存在一定误差。

    这篇文章就总结一下椭球七参数计算和高斯投影生反算的一些公式和方法,高斯投影正反算理论比较晦涩,可能存在某些地方说不清楚或比较乱的,做个参考就ok。


    2.Details

    (1)椭球参数

    常用椭球参数表
    (参数/椭球名称)CGCS 2000WGS 84
    GRS 80
    长半轴 a63781376378137
    6378137
    短半轴 b6356752.31414035586356752.3142451795
    6356752.3141403474
    极点子午圈曲率半径 c6399593.62586402326399593.6257584931
    6399593.6258640316
    椭球扁率 f0.00335281068118230.00335281066475
    0.00335281068118
    椭球扁率倒数 1/f
    298.257222101 
    298.257223563
    298.2572221012063‬
    第一偏心率 e10.08181919104281580.0818191908426215
    0.0818191910428318
    第二偏心率 e20.08209443815191720.0820944379496957
    0.0820944381519334
    椭球第一偏心率平方 e12 0.00669438002290080.0066943799901413
    0.0066943800229034
    椭球第二偏心率平方 e22 0.00673949677547900.0067394967422764
    0.0067394967754816

    椭球参数是参考椭球的几何参数,是高斯正反算和定义经度、纬度以及高程基础的重要依据,这里推荐一篇论文:2000国家大地坐标系椭球参数与GRS80和WGS84的比较_程鹏飞,里面只有三种椭球的参数,但是对于一般的生产和业务开发,CGCS2000和WGS84已经足够了。更多的椭球参数资料请看最后的下载链接。

     

    (2)大地坐标、平面坐标

    以下内容来自百度百科。

    大地坐标大地坐标(Geodetic coordinate)是大地测量中以参考椭球面为基准面的坐标,地面点P的位置用大地经度L、大地纬度B和大地高H表示。大地坐标多应用于大地测量学,测绘学等。而大地坐标系是大地测量中以参考椭球面基准面建立起来的坐标系。地面点的位置用大地经度、大地纬度和大地高度表示。大地坐标系的确立包括选择一个椭球、对椭球进行定位和确定大地起算数据。其实就是经纬度坐标。

    平面坐标表示点的平面位置。中国一般采用以高斯-克吕格投影分带的中央子午线为纵轴和赤道的投影为横轴的高斯-克吕格平面直角坐标系,简称高斯平面坐标系。坐标纵轴为x,自原点向北为正;坐标横轴为y,自原点向东为正。点的平面坐标为(x,y)。平面坐标系统是确定地面点平面位置所采用的参考系。大地测量、矿区控制测量、地形测量和工程测量所采用的平面坐标系主要有:高斯-克吕格平面直角坐标系、地方(矿区) 平面坐标系、独立坐标系。


    3.Theory

    (1)利用椭球长半轴和短半轴计算其它参数

    已知椭球长半轴 a,短半轴b,求点子午圈曲率半径 c、椭球扁率 f、椭球扁率倒数 1/f、第一偏心率 e1、第二偏心率 e2、椭球第一偏心率平方 e12、椭球第二偏心率平方 e22,并对应原始椭球参数表列出公式。

    椭球参数计算公式与对照结果表
    (参数/结果)CGCS 2000***计算公式
    输出结果(Python科学计算)
    差值
    长半轴 a6378137————
    ————
    ————
    短半轴 b6356752.3141403558————————————
    极点子午圈曲率半径 c6399593.6258640232c=\frac{a^{2}}{b}6399593.6258640230.0000000002
    椭球扁率 f0.0033528106811823f=\frac{a-b}{a}0.00335281068118227240.00000000000000003
    椭球扁率倒数 1/f
    298.257222101 
    \frac{1}{f}=\frac{a}{a-b}298.25722210100410.000000000004
    第一偏心率 e10.0818191910428158e=\sqrt{\frac{a^{2} - b^{2}}{a^{2}}}0.081819191042815170.0000000000000007
    第二偏心率 e20.0820944381519172{e}'=\sqrt{\frac{a^{2} - b^{2}}{b^{2}}}0.082094438151916580.0000000000000007
    椭球第一偏心率平方 e12 0.0066943800229008e^{2}=\frac{a^{2} - b^{2}}{a^{2}}0.00669438002290068550.00000000000000012
    椭球第二偏心率平方 e22 0.0067394967754790{e}'^{2}=\frac{a^{2} - b^{2}}{b^{2}}0.0067394967754788570.00000000000000015

    可以看到计算的参数结果和原始的参数结果误差小到可以忽略不计,这里就由取舍的小数位数来决定最后的误差大小,所以公式基本正确,以下为C#代码实现。

    decimal a = 0.0M; // 椭球长半轴 a
    decimal b = 0.0M; // 椭球短半轴 b
    decimal c = pow(a, 2) / b; // 极点子午圈曲率半径 c
    decimal e1 = sqrt(pow(a, 2) - pow(b, 2)) / a; // 第一偏心率 e1
    decimal e2 = sqrt(pow(a, 2) - pow(b, 2)) / b; // 第二偏心率 e2
    decimal e12 = pow(a, 2) - pow(b, 2) / pow(a, 2); // 椭球第一偏心率平方 e12
    decimal e22 = pow(a, 2) - pow(b, 2) / pow(b, 2); // 椭球第二偏心率平方 e22 
    decimal f = (a - b) / a; // 椭球扁率 f 
    decimal f_ =  a / (a - b); // 椭球扁率倒数 1/f 

    计算结果说明只需要知道长半轴 a 和短半轴 b 两个椭球参数值,就可以计算出其他的椭球参数。关于其他更多椭球参数表,之后会放到网盘上,提供下载链接(网上或者某个软件里找的)。

     

    (2)高斯投影正反算

    高斯投影正反算是大地坐标和平面坐标的相互转换方法。不过在查阅论文的时候,发现有论文提到高斯投影正反算也可以做大地坐标和高斯平面的直角坐标之间的转换。由于道行微末,就不做深究了,这里只涉及到大地坐标和平面坐标。

    ①高斯投影正算:大地坐标/经纬度坐标(L, B) ==> 平面坐标(X, Y)。

    \large x=X+\frac{N}{2{\rho }''^{2}}sinBcosB{l}''^{2} + \frac{N}{24{\rho }''^{2}}sinBcos^{3}B(5-t^{2}+9\eta ^{2}+4\eta ^{4}){l}''^{4}+\frac{N}{720{\rho }''^{6}}sinBcos^{5}B(61-58t^{2}+t^{4}){l}''^{6}

    \large y=\frac{N}{ {\rho}''}cosB{l}''+\frac{N}{6 {\rho}''^{3}}cos^{3}B(1-t^{2}+\eta ^{2}){L}''^{5}+\frac{N}{120{\rho }''^{5}}cos^{5}B(5-18t^{2}+t^{4}+14\eta ^{2}-58\eta ^{2}t^{2}){l}''^{5}

    《大地测量学基础》高斯投影正反算章节中,发现可以用查表的方式进行简化。描述:在高斯投影坐标计算的实际工作中,往往采用查表电算两种方式。为了便于查表计算,编有专门的计算用表,例如《高斯-克吕格投影计算表》(纬度0°~60°)。该表适用于高斯投影正算x,y精确至0.001m;高斯投影反算L,B精确至0.0001″。上表及其他有关用表都是对于克拉索夫斯基椭球而言的。公式简化如下:

    \large x=X+{l}'(a_{1}+a_{2}{l}')+\delta _{x}

    \large y={l}''(b_{1}+b_{2}{l}')+\delta _{y}

    其中,在编表时引用下列符号:

    {l}'={l}''^{2}\cdot 10^{-8}

    a_{1}=\frac{N}{2{\rho}''^{2}}sinBcosB10^{8}

    a_{2}=\frac{N}{24{\rho }''^{4}}sinBcos^{3}B(5-t^{2}+9\eta ^{2}+4\eta ^{4})10^{16}

    b_{1}=\frac{N}{​{\rho }''}cosB

    b_{2}=\frac{N}{6{\rho }''^{3}}cos^{3}B(1-t^{2}+\eta ^{2})10^{8}

    \delta _{x}=\frac{N}{720{\rho }''^{6}}sinBcos^{5}B(61-58t^{2}+t^{4}){l}''^{6}

    \delta _{y}=\frac{N}{120{\rho }''^{5}}cos^{5}B(5-18t^{2}+t^{4}+14\eta ^{2}-58\eta ^{2}t^{2}){l}''^{5}

    上式即为查表计算时的使用公式。式中\large X\large a_{1}\large a_{2}\large b_{1}\large b_{2} 都可 以纬度为引数,\large \delta _{x}\large \delta _{y} 以\large B\large l 为引数,在该表中查取。其中 \large X 和 \large b_{1} 需要进行二次内插。

    公式转换为代码,网上的大都大同小异,这里扒了一段作为代码实现,已经验证过了,转换结果基本正确。

    实体代码:

        /// <summary>
        /// 高斯点实体
        /// </summary>
        public class GaussPoint
        {
            public GaussPoint()
            {
            }
    
            /// <summary>
            /// 高斯点实体
            /// </summary>
            /// <param name="x">平面坐标X</param>
            /// <param name="y">平面坐标Y</param>
            public GaussPoint(double x, double y)
            {
                this.X = x;
                this.Y = y;
            }
    
            /// <summary>
            /// 平面坐标X
            /// </summary>
            public double X { get; set; }
    
            /// <summary>
            /// 平面坐标Y
            /// </summary>
            public double Y { get; set; }
        }

    高斯正算代码:

            /// <summary>
            /// 高斯正算方法(B, L) => (x, y)
            /// </summary>
            /// <param name="B">经度</param>
            /// <param name="L">纬度</param>
            /// <returns>高斯点实体</returns>
            public GaussPoint GetXYFromBL(decimal B, decimal L)
            {
                decimal x = 0, y = 0;
                // 计算临时值
                decimal e4 = pow(e12, 2); // e2 * e2;
                decimal e6 = pow(e12, 3); // e4 * e2;
                decimal e8 = pow(e12, 4); // e6 * e2;
                decimal e10 = pow(e12, 5); // e8 * e2;
                decimal e_12 = pow(e12, 6); // e10 * e2;
    
                decimal A1 = 1 + 3 * e12 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 / 1048576;
                decimal B1 = 3 * e12 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288;
                decimal C1 = 15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608;
                decimal D1 = 35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e_12 / 1048576;
                decimal E1 = 315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e_12 / 8388608;
                decimal F1 = 693 * e10 / 1310720 + 9009 * e_12 / 5242880;
                decimal G1 = 1001 * e_12 / 8388608;
    
                // 中央子午线 度转为秒值 如=105*3600
                decimal l0 = L0 * 3600; // L0为中央子午线,需要自己设置
                // 转为秒值
                decimal LL = L * 3600;
                // 转为弧度值 b
                decimal t_B = B * this.pi / 180;
                // 转为秒值 l 
                decimal t_L = (LL - l0) / p;
                decimal L2 = pow(t_L, 2); // t_L * t_L;
                decimal L4 = pow(t_L, 4); // L2 * L2;
    
                decimal sinB = sin(t_B);
                decimal sinB2 = sinB * sinB;
                decimal W = sqrt(1 - e12 * sinB2);
                decimal N = a / W;
                decimal t = tan(t_B);
                decimal t2 = t * t;
                decimal t4 = t2 * t2;
                decimal cosB = cos(t_B);
                decimal cosB2 = cosB * cosB;
                decimal cosB4 = cosB2 * cosB2;
                decimal y2 = e22 * cosB2; // η2
                decimal y4 = y2 * y2;
                decimal l_p = t_L; // t_L/p,上面t_L已经除了p值,这里就不再除p值
                decimal l2_p2 = L2; // l2/(p*p);
                decimal l4_p4 = L4; // l4/(p*p*p*p);
    
                decimal a0 = a * (1 - e12);
                // 计算子午弧长公式xx
                decimal xx = a0 * (A1 * t_B - B1 * sin(2 * t_B) + C1 * sin(4 * t_B) - D1 * sin(6 * t_B) + E1 * sin(8 * t_B) - F1 * sin(10 * t_B) + G1 * sin(12 * t_B));
                // 计度平面坐标值x,y
                x = xx + N * t * cosB2 * l2_p2 * (0.5M + (5 - t2 + 9 * y2 + 4 * y4) * cosB2 * l2_p2 / 24 + (61 - 58 * t2 + t4) * cosB4 * l4_p4 / 720);
                y = N * cosB * l_p * (1 + (1 - t2 + y2) * cosB2 * l2_p2 / 6 + (5 - 18 * t2 + t4 + 14 * y2 - 58 * y2 * t2) * cosB4 * l4_p4 / 120);
    
                // 转为高斯投影是否为大数投影? 即Zone 35带数(小数投影为CM_105E)
                if (IsBigNumberProj == true)
                {
                    y = y + (L0 / (int)Strip) * 1000000M; // L0为中央子午线,Strip为分度带(3或者6)
                }
                y = y + 500000.0M;
                return new GaussPoint(toNum(x), toNum(y));
            }

     

    ②高斯投影反算:平面坐标(X, Y) ==> 大地坐标/经纬度坐标(L, B)。

     

    \large B=B_{f}-\frac{t_{f}}{2M_{f}N_{f}cosB_{f}}y^{2}+\frac{t_{f}}{24M_{f}N^{3}_{f}}(5+3t^{2}_{f}+\eta ^{2}_{f}-9\eta ^{2}_{f}t^{2}_{f})y^{4}-\frac{t_{f}}{720M_{f}N^{5}_{f}}(61+90t^{2}_{f}+45t^{4}_{f})y^{6}

    \large L=\frac{1}{N_{f}cosB_{f}}y-\frac{1}{6N^{3}_{f}cosB_{f}}(1+2t^{2}_{f}+\eta ^{2}_{f})y^{3}+\frac{1}{120N^{5}_{f}cosB_{f}}(5+28t^{2}_{f}+24t^{4}_{f}+6\eta ^{2}_{f}+8\eta ^{2}_{f}t^{2}_{f})y^{5}

    B和L的单位为弧度(度分秒)。当L<3.5°时,上式换算经度达0.0001″。要让换算精确度至0.01″,可以简化如下:

          \large B=B_{f}-\frac{t_{f}}{2M_{f}N_{f}}y^{2}+\frac{t^{3}_{f}}{24M_{f}N^{3}_{f}}(5+3t^{2}_{f}+\eta ^{2}_{f}-9\eta ^{2}_{f}t^{2}_{f})y^{4}

    \large L=\frac{1}{N_{f}cosB_{f}}y-\frac{1}{6N^{3}_{f}cosB_{f}}(1+2t^{2}_{f}+\eta ^{2}_{f})y^{3}+\frac{1}{120N^{5}_{f}cosB_{f}}(5+28t^{2}_{f}+24t^{4}_{f})y^{5}

    公式转换为代码,网上的大都大同小异,这里扒了一段作为代码实现,已经验证过了,转换结果基本正确。

    实体代码:

        /// <summary>
        /// 经纬度实体
        /// </summary>
        public class Lnglat
        {
            public Lnglat()
            {
            }
    
            /// <summary>
            /// 经纬度实体
            /// </summary>
            /// <param name="lng">经度</param>
            /// <param name="lat">纬度</param>
            public Lnglat(double lng, double lat)
            {
                this.Lat = lat;
                this.Lng = lng;
            }
    
            /// <summary>
            /// 经度
            /// </summary>
            public double Lng { get; set; }
    
            /// <summary>
            /// 纬度
            /// </summary>
            public double Lat { get; set; }
        }

    高斯反算代码:

            /// <summary>
            /// 高斯反算方法 (x,y)=>(B,L)
            /// </summary>
            /// <param name="x">高斯点X</param>
            /// <param name="y">高斯点Y</param>
            /// <returns>返回经纬度实体</returns>
            public Lnglat GetBLFromXY(decimal x, decimal y)
            {
                decimal B = 0, L = 0;
    
                // 去掉大数和东移500公里
                decimal y1 = y - 500000.0M;
                if (IsBigNumberProj == true)
                {
                    y1 = y1 - (L0 / (int)Strip) * 1000000M; // L0为中央子午线,Strip为分度带(3或者6)
                }
                y = y1;
                // 中央子午线转为秒值 如=105*3600
                decimal l0 = L0 * 3600; // L0为中央子午线,需要自己设置
    
                // 计算临时值
                decimal e4 = pow(e12, 2); // e2 * e2;
                decimal e6 = pow(e12, 3); // e4 * e2;
                decimal e8 = pow(e12, 4); // e6 * e2;
                decimal e10 = pow(e12, 5); // e8 * e2;
                decimal e_12 = pow(e12, 6); // e10 * e2;
    
                decimal A1 = 1 + 3 * e12 / 4 + 45 * e4 / 64 + 175 * e6 / 256 + 11025 * e8 / 16384 + 43659 * e10 / 65536 + 693693 * e_12 / 1048576;
                decimal B1 = 3 * e12 / 8 + 15 * e4 / 32 + 525 * e6 / 1024 + 2205 * e8 / 4096 + 72765 * e10 / 131072 + 297297 * e_12 / 524288;
                decimal C1 = 15 * e4 / 256 + 105 * e6 / 1024 + 2205 * e8 / 16384 + 10395 * e10 / 65536 + 1486485 * e_12 / 8388608;
                decimal D1 = 35 * e6 / 3072 + 105 * e8 / 4096 + 10395 * e10 / 262144 + 55055 * e_12 / 1048576;
                decimal E1 = 315 * e8 / 131072 + 3465 * e10 / 524288 + 99099 * e_12 / 8388608;
                decimal F1 = 693 * e10 / 1310720 + 9009 * e_12 / 5242880;
                decimal G1 = 1001 * e_12 / 8388608;
    
                // 求底点纬度值Bf
                decimal B0 = x / (a * (1 - e12) * A1);
                decimal Bf = 0.0M;
                decimal FB = 0.0M;
                decimal FB1 = 0.0M;
                decimal a0 = a * (1 - e12);
                decimal delta = Math.Abs(Bf - B0);
                while (delta > 4.8E-11M) // 0.000000000048M
                {
                    Bf = B0;
                    FB = a0 * (A1 * Bf - B1 * sin(2 * Bf) + C1 * sin(4 * Bf) - D1 * sin(6 * Bf) + E1 * sin(8 * Bf) - F1 * sin(10 * Bf) + G1 * sin(12 * Bf));
                    FB1 = a0 * (A1 - 2 * B1 * cos(2 * Bf) + 4 * C1 * cos(4 * Bf) - 6 * D1 * cos(6 * Bf) + 8 * E1 * cos(8 * Bf) - 10 * F1 * cos(10 * Bf) + 12 * G1 * cos(12 * Bf));
                    B0 = Bf + (x - FB) / FB1;
                    //
                    delta = Math.Abs(Bf - B0);
                }
    
                decimal sinBf = sin(Bf);
                decimal sinBf2 = sinBf * sinBf;
                decimal W = sqrt(1 - e12 * sinBf2);
                decimal W3 = W * W * W;
                decimal N = a / W;
                decimal t = tan(Bf);
                decimal t2 = t * t;
                decimal t4 = t2 * t2;
                decimal cosBf = cos(Bf);
                decimal cosBf2 = cosBf * cosBf;
                decimal yy = e22 * cosBf2; // η2
                decimal Mf = a0 / W3;
                decimal y_N = y / N;
                decimal y_N2 = y_N * y_N;
                decimal y_N4 = y_N2 * y_N2;
    
                // 计算经伟度值B,L
                decimal t_B = Bf * p - (p * t / (2 * Mf) * y * y_N) * (1 - (5 + 3 * t2 + yy - 9 * yy * t2) * y_N2 + (61 + 90 * t2 + 45 * t4) * y_N4 / 360);
                decimal t_L = (p / cosBf) * y_N * (1 - (1 + 2 * t2 + yy) * y_N2 / 6 + (5 + 28 * t2 + 24 * t4 + 6 * yy + 8 * yy * t2) * y_N4 / 120);
    
                // 计算经纬度
                L = t_L + l0;
                // 转为度
                B = t_B / 3600;
                // 转为度
                L = L / 3600;
                return new Lnglat(toNum(L), toNum(B));
            }

    *正反算推算过程和实例都在参考资料中,具体可以参考《大地测量学基础》,里面有非常详细的推演过程,但是原理讲解我觉得不是特别明白,有点让人一头雾水,可能是我道行微末的原因。

    注意1:正反算方法中的e12,e22等参数,对应的是上面椭球参数。

    注意2:正反算方法里的pow、sqrt、tan、toNum等方法,对应的是Math.Pow、Math.Sqrt等C#内置函数的方法,这里只是做了封装。

     

    *自己实现了下七参数和高斯正反算的功能


    4.Environment and Configuration

    Environment:Windows 7及以上

    Language:C#

    IDE:Visual Studio 2019


    5.Conclusion

    高斯正反算原理弄不清楚的,明白个大概就好(不包括大佬),关键是要知道怎么用、用来干嘛的、怎么调整参数等等。

    附上所用的参考资料:https://pan.baidu.com/s/1Asw1hYVR42GD9oBI3QdJqw 提取码:5qrc

     

     

    展开全文
  • 参考椭球的 a,b可变带来的直接问题便是大地计算模型中X和B f 参数难于求解 。本文提出了参考椭球的 a,b可变的一种大地正算和大地反算的新算法 。推导了算法的基本公式并设计了计算流程,在此基础上,编制了相应的...
  • 将笛卡尔坐标转换为大地坐标在三轴椭球或双轴椭球或球面上(x/a)^2+(y/b)^2+(z/c)^2=1 三轴(一般)椭球方程笛卡尔到大地坐标 xyz ==> BL h 参数: * X, [xyz] - 笛卡尔坐标数据,nx 3 矩阵或三个 nx 1 向量* 轴,[a;...
  • 一、不同椭球基准的坐标转换 ①原坐标需先转换成所在椭球基准的空间直角坐标(X1Y1Z1) ②通过“七参转换”转换成目标所在椭球基准的空间直角坐标(X2Y2Z2) ③相同椭球基准的坐标转换,得到所需坐标 ...
  • 1975 年,Vincenty 发表了一种快速收敛算法,用于计算椭球体地球上点之间的距离。 该算法精确到几毫米以内。 从那以后,他的算法在大地测量学和工程学中得到了重要的应用。... 请参阅代码注释以获取参考
  • 针对由于大地高所引起的地面区域实际球面...得到以地面区域的平均大地高为参数将其在参考椭球面上的面积改化为地面实际球面面积的实用方法.研究结果表明:该方法计算简便、精度高,面积修正的相对误差达百万分之一量级.
  • 参考 [1]和Aude Billard,“协方差矩阵的几何变换不变非参数聚类及其在无监督联合分割和动作发现中的应用”。 机器学习研究。电工技术学报。 依存关系 :汤姆·明卡(Tom Minka)的库,其中包含高度优化的数学函数...
  • 在ArcGIS10.2 的ArcTools/Data Management Tools/Projections and Transformations/CreateCustomGeoTransformation可以自定义椭球转换参数,可以选用CoordFrame(七参数)或者GeoCentric_Translation(三参数) ...
  • //获取椭球名称 string strProjectionName=projectedCoordinateSystem.Projection.Name;//投影名称 double dblScaleFactor=projectedCoordinateSystem.ScaleFactor;//尺度因子 double dblLatitudeOfOrigin=...
  • CGCS2000国家大地坐标系参数

    万次阅读 2018-07-27 11:26:04
    国家大地坐标系的定义包括坐标系的原点、三个坐标轴的指向、尺度以及地球椭球的4个基本参数的定义。 2000国家大地坐标系 原点:包括海洋和大气的整个地球的质量中心 Z轴:由原点指向历元2000.0的地球参考极的方向...
  • 而全球定位系统(GPS)应用的是WGS-84系椭球参数。 1.大地基准面(引用至 http://wenku.baidu.com/view/061c40c7aa00b52acfc7cad1.html )  大地基准面(Geodetic datum),设计用为最密合部份或全部...
  • 本章讲述地球椭球与参考椭球的概念,进而介绍椭球的基本几何参数,基本坐标系及其相互关系。同时,讲述椭球面同地面之间的关系,如何将地面观测元素(水平方向及斜距等)归算至椭球面上。在对本章的学习中,要建立起...
  • ArcGIS空间参考

    千次阅读 2016-07-15 10:36:56
    ArcGIS空间参考   ... 空间参考描述了一个地物在地球上的真实位置。为了正确的对位置进行描述,需要引入一个可供...而地球是一个不规则形状的椭球体,那么使用什么样的方法来模拟地球的形状,又该如何将球面上
  • GIS椭球转换

    千次阅读 2013-04-18 14:37:18
    继上篇文章所讲,GIS数据有是空间参考的,否则该数据也就失去了GIS...地理坐标系有参考基准面描述,基准面由参考椭球体和七参数描述;投影坐标系由地理坐标系和投影参数描述。  为什么要定义七参数才可以定义基准面?
  • 文章目录参考介绍等角投影公式 参考 参考书目: 《地图投影原理与方法》。 英文版的一本书《与大地测量和制图有关的纬度发展》,里面有P85页有下面提到的公式。 介绍 在推求地图投影方程时,通常有两种情况(参考...
  • : 链接 长方形 圆 圆角长方形 菱形 关于 Mermaid 语法,参考 这儿, FLowchart流程图 我们依旧会支持flowchart的流程图: Created with Raphaël 2.3.0 开始 我的操作 确认? 结束 yes no 关于 Flowchart流程图 语法...
  • 比如我们常用的四种坐标系北京54、西安80、WGS84、CGCS2000所对应的椭球体,它们的椭球参数就各不相同。而不同空间直角坐标系之间的转换一般通过七参数或者四参数来实现坐标转换。 二.四参数转换 两个不同的二...
  • GIS的空间参考系统

    千次阅读 2018-11-26 17:03:30
    2.2、参考椭球体 2.3、大地原点 三、地图投影 3.1、地图投影类型 3.2、地图投影参数 3.2.1、标准线 3.2.2、比例尺 3.2.3、中心线 3.2.4、坐标偏移 3.3、常用地图投影 3.3.1、横轴墨卡托投影 3.3.2、...
  • 现阶段GPS-RTK以WGS-84 坐标系统为主流,所发布的星历参数也是基于此坐标系统,但随着北斗导航系统的逐步完善,我国测量仪器正在向国家2000 椭球过渡。现阶段我国常用的大地平面坐标系统有国家2000坐标系、西安80...
  • 电子罗盘的椭球拟合与椭球变换 椭球拟合——最小二乘法 由电子罗盘采集数据如何得到椭球方程?
  • GIS中的空间参考

    2019-03-15 09:26:22
    但如果应用不一样,在实际中选择的空间参考系也是不一样的。例如我们使用GIS系统在做一个房间的布置的时候,就应该不会经纬度和高程数据来标识物体的位置。 但其实我们在GIS行业中主要研究的还是比较大的一片区域,...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,530
精华内容 1,012
关键字:

参考椭球参数