「已注销」 2022-08-13 09:44 采纳率: 0%
浏览 272
已结题

求高斯坐标正反算,要求如图

img


批处理就是文件形式读取数据,四个子程序可以不同函数放一起
ssuididjhshjjdidudh

  • 写回答

4条回答 默认 最新

  • 鹅毛在路上了 Matlab领域优质创作者 2022-08-15 09:22
    关注

    以下是Matlab实现高斯正反算的两个方法:

    方法1,教程上提供的的高斯正反算例程,函数1:高斯正算

    function PRO = GKZS(Lat,Lon,L0,a,f)
    % Lat: Latitude(rad) 
    % Lon: longitude(rad)
    % REF//程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M].
    %    //北京:测绘出版社,2008:144-148.
    Lat = Lat*pi/180;
    Lon = Lon*pi/180;
    MedLon = L0*pi/180; %中央子午线经度
    Eth.R0 = a;%长半轴
    Eth.f = f; % 扁率
    Eth.e12 = 2*Eth.f - Eth.f*Eth.f; % 第一偏心率额e^2
    Eth.e22 = Eth.e12/((1 - Eth.f)*(1 - Eth.f));% 第二偏心率额e'^2
    %% 高斯投影正算公式
    RN = Eth.R0/sqrt(1 - Eth.e12*sin(Lat)*sin(Lat));
    Lon = Lon - MedLon;
    Lon2 = Lon*Lon;
    Lon4 = Lon2*Lon2;
    tnLat = tan(Lat);
    tn2Lat = tnLat*tnLat;
    tn4Lat = tn2Lat*tn2Lat;
    csLat = cos(Lat);
    cs2Lat = csLat*csLat;
    cs4Lat = cs2Lat*cs2Lat;
    Eta2 = Eth.e22*cs2Lat;
    NTBLP = RN*tnLat*cs2Lat*Lon2;
    coe1 = (5 - tn2Lat + 9*Eta2 + 4*Eta2*Eta2)*cs2Lat*Lon2/24;
    coe2 = (61 - 58*tn2Lat + tn4Lat)*cs4Lat*Lon4/720;
    x = Merdian(Eth,Lat) + NTBLP*(0.5 + coe1 + coe2);
    NBLP = RN*csLat*Lon;
    coe3 = (1 - tn2Lat + Eta2)*cs2Lat*Lon2/6;
    coe4 = (5 - 18*tn2Lat + tn4Lat + 14*Eta2 - 58*tn2Lat*Eta2)*cs4Lat*Lon4/120;
    y = NBLP*(1 + coe3 + coe4) + 500000;
    PRO=[x y]
    end
    
    function X0 = Merdian(Eth,Lat)
    % REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.
    S0 = Eth.R0*(1 - Eth.e12);
    e2 = Eth.e12;
    e4 = e2*e2;
    e6 = e4*e2;
    e8 = e6*e2;
    e10 = e8*e2;
    e12 = e10*e2;
    A1 = 1 + 3*e2/4 + 45*e4/64 + 175*e6/256 + 11025*e8/16384 + 43659*e10/65536 + 693693*e12/1048576;
    B1 = 3*e2/8 + 15*e4/32 + 525*e6/1024 + 2205*e8/4096 + 72765*e10/131072 + 297297*e12/524288;
    C1 = 15*e4/256 + 105*e6/1024 + 2205*e8/16384 + 10395*e10/65536 + 1486485*e12/8388608;
    D1 = 35*e6/3072 + 105*e8/4096 + 10395*e10/262144 + 55055*e12/1048576;
    E1 = 315*e8/131072 + 3465*e10/524288 + 99099*e12/8388608;
    F1 = 693*e10/1310720 + 9009*e12/5242880;
    G1 = 1001*e12/8388608;
    X0 = S0*(A1*Lat - B1*sin(2*Lat) + C1*sin(4*Lat) - D1*sin(6*Lat) +...
        E1*sin(8*Lat) - F1*sin(10*Lat) + G1*sin(12*Lat));
    end
    

    函数2:高斯反算

    function BL = GKFS(x,y,L0,a,f)
    
    iPI = 0.0174532925199433; %3.1415926535898/180.0; 
    Eth.a = a ;%长半轴
    Eth.f = f;% 扁率
    e2 = 2*Eth.f-Eth.f*Eth.f;
    e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));
    ee = e2/(1-e2);
    X = y
    Y = x
    longitude0 = L0*pi/180; %中央子午线经度
    %
    ZoneWide = 6; %6度带宽 
    ProjNo = fix(X/1000000) ; %查找带号
    X0 = ProjNo*1000000+500000; 
    Y0 = 0; 
    xval = X-X0; yval = Y-Y0;%带内大地坐标
    M = yval;
    u = M/(Eth.a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
    fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);
    C = ee*cos(fai)*cos(fai);
    T = tan(fai)*tan(fai);
    NN = Eth.a/sqrt(1.0-e2*sin(fai)*sin(fai));
    R = Eth.a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*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)/cos(fai);
    latitude1 = fai -(NN*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
    L = longitude1 / iPI
    B = latitude1 / iPI;
    L = fix(L)+fix((L-fix(L))*60)/100+((L-fix(L))*60-fix((L-fix(L))*60))*60/10000
    B = fix(B)+fix((B-fix(B))*60)/100+((B-fix(B))*60-fix((B-fix(B))*60))*60/10000
    BL = [B L]
    
    end
    

    方法2的高斯正算函数:

    function [ x,y ] = BL2XY( B,L,a,e,k )
    
    ee=e/sqrt(1-e^2);
    % 计算中央子午线经度
    if k==6
        L0=6*(fix(L/6)+1)-3;
    elseif k==3
        L0=3*(fix((L-1.5)/3)+1);
    end
    LP=(L-L0)/180*pi;
    B=B/180*pi;
    
    % 计算子午圈弧长
    ap=1+3/4*e^2+45/64*e^4+175/256*e^6+11025/16384*e^8+43659/65536*e^10;
    bp=3/4*e^2+15/16*e^4+525/512*e^6+2205/2048*e^8+72765/65536*e^10;
    cp=15/64*e^4+105/256*e^6+2205/4096*e^8+10395/16384*e^10;
    dp=35/512*e^6+315/2048*e^8+31185/131072*e^10;
    ep=315/16384*e^8+3465/65536*e^10;
    fp=639/131072*e^10;
    X=a*(1-e^2)*(ap*B-bp/2*sin(2*B)+cp/4*sin(4*B)-dp/6*sin(6*B)+ep/8*sin(8*B)-fp/10*sin(10*B));
    
    % 计算高斯平面坐标
    N=a/sqrt(1-e^2*sin(B)*sin(B));
    eta=ee*cos(B);
    t=tan(B);
    m=cos(B)*LP;
    x=X+N*t*(0.5*m^2+1/24*(5-t^2+9*eta^2+4*eta^4)*m^4+1/720*(61-58*t^2+t^4)*m^6);
    y=N*(m+1/6*(1-t^2+eta^2)*m^3+1/120*(5-18*t^2+t^4+14*eta^2-58*eta^2*t^2)*m^5);
    
    end
    

    高斯反算函数:

    function [ B,L ] = XY2BL( x,y,a,e,L0 )
    
    % CGCS2000椭球参数
    % a=6378137;
    % e=0.081819191042811;
    % ee=0.0820944381518626;
    % 1975国际椭球参数
    % a=6378140;
    % e=0.0818192214555235;
    
    ee=e/sqrt(1-e^2);
    % 计算底点的经度
    ap=1+3/4*e^2+45/64*e^4+175/256*e^6+11025/16384*e^8+43659/65536*e^10;
    bp=3/4*e^2+15/16*e^4+525/512*e^6+2205/2048*e^8+72765/65536*e^10;
    cp=15/64*e^4+105/256*e^6+2205/4096*e^8+10395/16384*e^10;
    dp=35/512*e^6+315/2048*e^8+31185/131072*e^10;
    ep=315/16384*e^8+3465/65536*e^10;
    fp=639/131072*e^10;
    
    B=x/(a*(1-e^2)*ap);
    deltaB=1000;
    while deltaB>10^(-20)
        Bi=1/ap*(x/(a*(1-e^2))+bp/2*sin(2*B)-cp/4*sin(4*B)+dp/6*sin(6*B)-ep/8*sin(8*B)+fp/10*sin(10*B));
        deltaB=abs(Bi-B);
        B=Bi;
    end
    Bf=B;
    
    % 计算大地经纬度BL
    tf=tan(Bf);
    Vf=sqrt(1+ee^2*cos(Bf)*cos(Bf));
    Wf=sqrt(1-e^2*sin(Bf)*sin(Bf));
    Nf=a/Wf;
    etaf=ee*cos(Bf);
    B=Bf-1/2*Vf^2*tf*((y/Nf)*(y/Nf)-1/12*(5+3*tf^2+etaf^2-9*etaf^2*tf^2)*(y/Nf)^4+1/360*(61+90*tf^2+45*tf^4)*(y/Nf)^6);
    LP=1/cos(Bf)*((y/Nf)-1/6*(1+2*tf^2+etaf^2)*(y/Nf)^3+1/120*(5+28*tf^2+24*tf^4+6*etaf^2+8*etaf^2*tf^2)*(y/Nf)^5);
    B=B/pi*180;
    LP=LP/pi*180;
    L=L0+LP;
    
    
    end
    

    高斯正反算批量计算及GUI设计资料,需要的话私信可以发你邮箱

    评论 编辑记录

报告相同问题?

问题事件

  • 已结题 (查看结题原因) 8月16日
  • 创建了问题 8月13日

悬赏问题

  • ¥15 机器学习教材中的例题询问
  • ¥15 求.net core 几款免费的pdf编辑器
  • ¥15 C# P/Invoke的效率问题
  • ¥20 thinkphp适配人大金仓问题
  • ¥20 Oracle替换.dbf文件后无法连接,如何解决?(相关搜索:数据库|死循环)
  • ¥15 数据库数据成问号了,前台查询正常,数据库查询是?号
  • ¥15 算法使用了tf-idf,用手肘图确定k值确定不了,第四轮廓系数又太小才有0.006088746097507285,如何解决?(相关搜索:数据处理)
  • ¥15 彩灯控制电路,会的加我QQ1482956179
  • ¥200 相机拍直接转存到电脑上 立拍立穿无线局域网传
  • ¥15 (关键词-电路设计)