1 | function [alat,alon]=mapxy(x,y,hem,varargin);
|
---|
2 | %MAPXY - compute latitude and longitude from x and y
|
---|
3 | %
|
---|
4 | % hemisphere must be 'n' for north and 's' for south
|
---|
5 | % One can specify the standard parallel and center meridien.
|
---|
6 | % Default values are (70,45) for northern hemisphere and (71,0) for southerm hemisphere.
|
---|
7 | %
|
---|
8 | % Usage:
|
---|
9 | % [latitude,longitude]=mapxy(x,y,hemisphere)
|
---|
10 | % [latitude,longitude]=mapxy(x,y,hemisphere,parallel,meridien)
|
---|
11 | %
|
---|
12 | % Example:
|
---|
13 | % [lat,long]=mapxy(x,y,'s')
|
---|
14 | % [lat,long]=mapxy(x,y,'n',71,39)
|
---|
15 | %
|
---|
16 | % See also MAPLL, LL2XY
|
---|
17 |
|
---|
18 | %check number of arguments.
|
---|
19 | if ~((nargin==3) | (nargin==5)),
|
---|
20 | mapxyerrorusage();
|
---|
21 | end
|
---|
22 |
|
---|
23 | %hem: must be either 's' or 'n'
|
---|
24 | if ~ischar(hem),
|
---|
25 | error('mapxy error message: hemisphere argument should be ''n'' or ''s''');
|
---|
26 | end
|
---|
27 |
|
---|
28 | %check hem is either 'n' or 's':
|
---|
29 | if ~(strcmpi(hem,'n') | strcmpi(hem,'s')),
|
---|
30 | error('mapxy error message: hem should be ''s'' or ''n''');
|
---|
31 | end
|
---|
32 |
|
---|
33 | %set sn:
|
---|
34 | if strcmpi(hem,'s'),
|
---|
35 | slat=71;
|
---|
36 | sn=-1.0;
|
---|
37 | xlam=0;
|
---|
38 | else
|
---|
39 | slat=70;
|
---|
40 | sn=1.0;
|
---|
41 | xlam=45;
|
---|
42 | end
|
---|
43 |
|
---|
44 | %set defaults for standard parallels and centre meridians.
|
---|
45 | if nargin==5,
|
---|
46 | slat=varargin{1};
|
---|
47 | xlam=varargin{2};
|
---|
48 | end
|
---|
49 |
|
---|
50 | re = 6378137;
|
---|
51 | e2 = 0.00669437999015;
|
---|
52 |
|
---|
53 | e=sqrt(e2);
|
---|
54 |
|
---|
55 | slat=slat/180*pi;
|
---|
56 |
|
---|
57 | rho=sqrt(x.^2+y.^2);
|
---|
58 | cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
|
---|
59 | t=tan((pi/4)-(slat/2.))./((1.0-e*sin(slat))./(1.0+e*sin(slat))).^(e/2.);
|
---|
60 |
|
---|
61 | t=rho.*t./(re*cm);
|
---|
62 |
|
---|
63 | chi=(pi/2.)-2.*atan(t);
|
---|
64 |
|
---|
65 | alat=chi+((e2/2.)+(5.0*e2^2/24.)+(e2^3/12.))*sin(2*chi)+...
|
---|
66 | ((7.0*e2^2/48.)+(29.*e2^3/240.))*sin(4.0*chi)+...
|
---|
67 | (7.0*e2^3/120.)*sin(6.0*chi);
|
---|
68 |
|
---|
69 | alat=(sn*alat/pi*180);
|
---|
70 | xpr=sn*x;
|
---|
71 | ypr=sn*y;
|
---|
72 | alon=atan2(xpr,-ypr)/pi*180-sn*xlam;
|
---|
73 | alon=sn*alon;
|
---|
74 |
|
---|
75 | indice=find(alon<0);
|
---|
76 | if length(indice)>0,
|
---|
77 | alon(indice)=alon(indice)+360;
|
---|
78 | end
|
---|
79 |
|
---|
80 | indice=find(alon>360);
|
---|
81 | if length(indice)>0,
|
---|
82 | alon(indice)=alon(indice)-360;
|
---|
83 | end
|
---|
84 |
|
---|
85 |
|
---|
86 | end
|
---|
87 |
|
---|
88 |
|
---|
89 | function mapxyerrorusage()
|
---|
90 | help mapxy
|
---|
91 | end
|
---|