WGS84转BD09坐标系--Matlab程序

一、本算法结合网上算法做了参数改变，实际误差几十米

function m=transformlng( lat, lon)
        a = 300.0 + lat + 2.0 * lon + 0.1 * lat * lat + 0.1 * lat * lon + 0.1 * sqrt(abs(lat));
        b= a+ ((20.0 * sin(6.0 * lat * pi) + 20.0 * sin(2.0 * lat * pi)) * 2.0 / 3.0);
        c= b+ ((20.0 * sin(lat * pi) + 40.0 * sin(lat / 3.0 * pi)) * 2.0 / 3.0);
        d= c+((150.0 * sin(lat / 12.0 * pi) + 300.0 * sin(lat / 30.0 * pi)) * 2.0 / 3.0);

       m=d;
end

 

function n= transformlat( lat, lon)     
         a = -100.0 + 2.0 * lat + 3.0 * lon + 0.2 * lon * lon + 0.1 * lat * lon + 0.2 * sqrt(abs(lat));
         b = a+((20.0 * sin(6.0 * lat * pi) + 20.0 * sin(2.0 * lat * pi)) * 2.0 / 3.0);
         c = b+((20.0 * sin(lon * pi) + 40.0 * sin(lon / 3.0 * pi)) * 2.0 / 3.0);
         d = c+ ((160.0 * sin(lon / 12.0 * pi) + 320 * sin(lon * pi / 30.0)) * 2.0 / 3.0);

       n= d;
end
 

  function y=WGS824TOBD09( lat_f,  lng_f)
  
           g=floor(lat_f/100);
           b=lat_f-g*100;
           c=b/60;

           lat=g+c;
           

           d=floor(lng_f/100);
           e=lng_f-d*100;
           f=e/60;

           lng=d+f;
        
           x_PI = 3.14159265358979324 * 3000.0 / 180.0;
    
           a = 6378245.0;
           ee = 0.00669342162296594323;

          
           dlat = transformlat(lng - 105.0, lat - 35.0);
           dlng = transformlng(lng - 105.0, lat - 35.0);
           radlat = lat / 180.0 * pi;
           magic = sin(radlat);
           magic = 1 - ee * magic * magic;
           sqrtmagic = sqrt(magic);
           dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi);
           dlng = (dlng * 180.0) / (a / sqrtmagic *cos(radlat) * pi);
           mglat = lat + dlat;
           mglng = lng + dlng;
 
          
           z = sqrt(mglng * mglng + mglat * mglat) + 0.00002 * sin(mglat * x_PI);
           theta = atan2(mglat, mglng) + 0.000003 * cos(mglng * x_PI);
           bd_lng = z * cos(theta) + 0.0065;
           bd_lat = z * sin(theta) + 0.006;
           
            x = bd_lng - 0.002558999; y = bd_lat - 0.0031078;    
            c = sqrt(x * x + y * y) - 0.00002 * sin(y * x_PI);    
            theta = atan2(y, x) - 0.000003 * cos(x * x_PI);    
            gg_lon = c * cos(theta);    
            gg_lat = c * sin(theta);
            
            y=[gg_lat,gg_lon];
    
  end