2020-08-11

Albers投影转化为wgs

x=pop[i,0]
y=pop[i,1]
xd=x
yd=y
pi=4atan(1);
print, pi
fi0=0pi/180;
lanbuda0=110pi/180;
fi1=25pi/180;
fi2=47pi/180;
EF=0;
NF=0;
a=6378137;
b=6356752.3142;
e1=sqrt(1-b/ab/a);
e2=1-e1e1;
alfa0=e2(sin(fi0)/(1-e1e1sin(fi0)sin(fi0))-0.5/e1alog((1-e1sin(fi0))/(1+e1sin(fi0))));
alfa1=e2*(sin(fi1)/(1-e1e1sin(fi1)sin(fi1))-0.5/e1alog((1-e1sin(fi1))/(1+e1sin(fi1))));
alfa2=e2*(sin(fi2)/(1-e1e1sin(fi2)sin(fi2))-0.5/e1alog((1-e1sin(fi2))/(1+e1sin(fi2))));
m1=cos(fi1)/sqrt(1-e1e1sin(fi1)sin(fi1));
m2=cos(fi2)/sqrt(1-e1e1sin(fi2)sin(fi2));
n=(m1m1-m2m2)/(alfa2-alfa1);
c=m1m1+nalfa1;
rou0=asqrt(c-nalfa0)/n;
sita=atan((xd-EF)/(rou0-(yd-NF)));
rou=sqrt((xd-EF)(xd-EF)+(rou0-(yd-NF))(rou0-(yd-NF)));
alfa=(c-(roun/a)(roun/a))/n;
beta=asin(alfa/(1-(1-(e1e1))/(2e1)alog((1-e1)/(1+e1))));
lon=lanbuda0+sita/n;
lat=beta+((e1e1)/3+31(e1e1e1e1)/180+517(e1e1e1e1e1e1)/5040)sin(2beta)+(23(e1e1e1e1)/360+251(e1e1e1e1e1e1)/3780)sin(4beta)+(761(e1e1e1e1e1e1)/45360)sin(6beta);
lon=lon180/pi-5;
lat=lat*180/pi;