Index: /issm/trunk-jpl/src/m/coordsystems/lambert2xy.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/lambert2xy.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/lambert2xy.m (revision 15139) @@ -0,0 +1,69 @@ +function [x,y] = lambert2xy(lat,lon,sgn,projection_center_lat,projection_center_lon) +%LAMBERT2XY - converts lat long from Lambert Azimuthal to Polar Stereographic +% +% Converts from geodetic latitude and longitude that are +% in Lambert Azimuthal (equal area) projections to Polar +% Stereographic (X,Y) coordinates for the polar regions. +% +% Usage: +% [x,y] = lambert2xy(lat,lon,sgn) +% [x,y] = lambert2xy(lat,lon,sgn,projection_center_lat,projection_center_lon) +% +% - provide lat in [-90,90] and lon in [-180,180]. + +% - sgn = +1 N hemisphere [default projection center lat = 90 lon=0] +% -1 S hemisphere [default projection center lat = -90 lon=0] + +%Get projection_center_lat and projection_center_lon +if nargin==5, + latitude0 = projection_center_lat; + longitude0 = projection_center_lon; +elseif nargin==3, + if sgn==1, + latitude0 = 90; longitude0 = 0; + disp('Info: creating coordinates in polar stereographic (Projection center lat: 90N lon: 0)'); + elseif sgn==-1, + latitude0 = -90; longitude0 = 0; + disp('Info: creating coordinates in polar stereographic (Projection center lat: 90S lon: 0)'); + else + error('Sign should be either +1 or -1'); + end +else + help lambert2xy + error('bad usage'); +end + +% Radius of the earth in meters +a = 6378137.0; +% Eccentricity of the Hughes ellipsoid squared +e = 0.081819191; + +% Projection center latitude and longitude in radians +phi0 = latitude0 * pi/180; +lam0 = longitude0 * pi/180; + +% Some constant based on phi0 and lam0 +qp= (1-e^2)*((1/(1-e^2))-((1/(2*e))*log((1-e)/(1+e)))); +q0=(1-e^2)*((sin(phi0)/(1-e^2*sin(phi0)*sin(phi0)))-((1/(2*e))*log((1-e*sin(phi0))/(1+e*sin(phi0))))); +Rq=a*sqrt(qp/2); +b0=asin(q0/qp); +D =a*(cos(phi0)/sqrt(1-e^2*sin(phi0)*sin(phi0)))/(Rq*cos(b0)); + +% Latitude and longitude in radians +phi = lat*pi/180; +lam = lon*pi/180; + +% Some other phi,lam dependent parameters +q=(1-e^2)*((sin(phi)/(1-e^2*sin(phi)*sin(phi)))-((1/(2*e))*log((1-e*sin(phi))/(1+e*sin(phi))))); +b =asin(q/qp); +B =Rq*sqrt(2/(1+sin(b0)*sin(b)+(cos(b0)*cos(b)*cos(lam-lam0)))); + +% Calculation of x and y +if(abs(lat)==90) + x=0.0; + y=0.0; +else + x=(B*D)*(cos(b)*sin(lam-lam0)); + y=(B/D)*((cos(b0)*sin(b))-(sin(b0)*cos(b)*cos(lam-lam0))); +end + Index: /issm/trunk-jpl/src/m/coordsystems/ll2mercator.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/ll2mercator.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/ll2mercator.m (revision 15139) @@ -0,0 +1,14 @@ +function [mx my]=ll2mercator(lat, lon), +%LL2MERCATOR - transform lat long to mercator projection +% +% Usage: +% [mx my]=ll2mercator(lat, lon) + +EARTH_RADIUS = 6378137; +EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS; +INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0; +ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0; + +mx = (lon * ORIGIN_SHIFT) / 180.0; +my = log(tan((90 + lat) * pi/360.0))/(pi/180.0); +my = (my * ORIGIN_SHIFT) /180.0; Index: /issm/trunk-jpl/src/m/coordsystems/ll2xy.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/ll2xy.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/ll2xy.m (revision 15139) @@ -0,0 +1,66 @@ +function [x,y] = ll2xy(lat,lon,sgn,central_meridian,standard_parallel) +%LL2XY - converts lat long to polar stereographic +% +% Converts from geodetic latitude and longitude to Polar +% Stereographic (X,Y) coordinates for the polar regions. +% Author: Michael P. Schodlok, December 2003 (map2ll) +% +% Usage: +% [x,y] = ll2xy(lat,lon,sgn) +% [x,y] = ll2xy(lat,lon,sgn,central_meridian,standard_parallel) +% +% - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70) +% -1 : south latitude (default is mer=0 lat=71) + +%Get central_meridian and standard_parallel depending on hemisphere +if nargin==5, + delta = central_meridian; + slat = standard_parallel; +elseif nargin==3 + if sgn == 1, + delta = 45; slat = 70; + disp('Info: creating coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); + elseif sgn==-1, + delta = 0; slat = 71; + disp('Info: creating coordinates in polar stereographic (Std Latitude: 71ºS Meridian: 0º)'); + else + error('Sign should be either +1 or -1'); + end +else + help ll2xy + error('bad usage'); +end + +% Conversion constant from degrees to radians +cde = 57.29577951; +% Radius of the earth in meters +re = 6378.273*10^3; +% Eccentricity of the Hughes ellipsoid squared +ex2 = .006693883; +% Eccentricity of the Hughes ellipsoid +ex = sqrt(ex2); + +latitude = abs(lat) * pi/180.; +longitude = (lon + delta) * pi/180.; + +% compute X and Y in grid coordinates. +T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2); + +if (90 - slat) < 1.e-5 + rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex)); +else + sl = slat*pi/180.; + tc = tan(pi/4.-sl/2.)/((1.-ex*sin(sl))/(1.+ex*sin(sl)))^(ex/2.); + mc = cos(sl)/sqrt(1.0-ex2*(sin(sl)^2)); + rho = re*mc*T/tc; +end + +y = -rho .* sgn .* cos(sgn.*longitude); +x = rho .* sgn .* sin(sgn.*longitude); + +[cnt1,cnt2] = find(latitude >= pi / 2.); + +if cnt1 + x(cnt1,1) = 0.0; + y(cnt1,1) = 0.0; +end Index: /issm/trunk-jpl/src/m/coordsystems/utm2ll.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/utm2ll.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/utm2ll.m (revision 15139) @@ -0,0 +1,117 @@ +function [Lat,Lon] = utm2ll(xx,yy,utmzone) +% ------------------------------------------------------------------------- +% [Lat,Lon] = utm2ll(x,y,utmzone) +% +% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84). +% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez. +% +% Inputs:on) +% -3.485713 7.801235 -119.955246 -17.759537 -94.799019 121.640266 +% +% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds +% [Lat, Lon]=utm2ll(x,y,utmzone); +% LatDMS=dms2mat(deg2dms(Lat)) +%LatDMS = +% 40.00 18.00 55.55 +% 46.00 17.00 2.01 +% 37.00 34.00 40.17 +% 28.00 38.00 44.33 +% 38.00 51.00 19.96 +% 25.00 3.00 42.41 +% LonDMS=dms2mat(deg2dms(Lon)) +%LonDMS = +% -3.00 29.00 8.61 +% 7.00 48.00 4.40 +% -119.00 57.00 18.93 +% -17.00 45.00 34.33 +% -94.00 47.00 56.47 +% 121.00 38.00 24.96 +% +% Author: +% Rafael Palacios +% Universidad Pontificia Comillas +% Madrid, Spain +% Version: Apr/06, Jun/06, Aug/06 +% Aug/06: corrected m-Lint warnings +%------------------------------------------------------------------------- + +% Argument checking +% +error(nargchk(3, 3, nargin)); %3 arguments required +n1=length(xx); +n2=length(yy); +n3=size(utmzone,1); +if (n1~=n2 || n1~=n3) + error('x,y and utmzone vectors should have the same number or rows'); +end +c=size(utmzone,2); +if (c~=4) + error('utmzone should be a vector of strings like "30 T"'); +end + +% Memory pre-allocation +% +Lat=zeros(n1,1); +Lon=zeros(n1,1); + +% Main Loop +% +for i=1:n1 + if (utmzone(i,4)>'X' || utmzone(i,4)<'C') + fprintf('utm2ll: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n'); + end + if (utmzone(i,4)>'M') + hemis='N'; % Northern hemisphere + else + hemis='S'; + end + + x=xx(i); + y=yy(i); + zone=str2double(utmzone(i,1:2)); + + sa = 6378137.000000 ; sb = 6356752.314245; + + % e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa; + e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb; + e2cuadrada = e2 ^ 2; + c = ( sa ^ 2 ) / sb; + % alpha = ( sa - sb ) / sa; %f + % ablandamiento = 1 / alpha; % 1/f + + X = x - 500000; + + if hemis == 'S' || hemis == 's' + Y = y - 10000000; + else + Y = y; + end + + S = ( ( zone * 6 ) - 183 ); + lat = Y / ( 6366197.724 * 0.9996 ); + v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996; + a = X / v; + a1 = sin( 2 * lat ); + a2 = a1 * ( cos(lat) ) ^ 2; + j2 = lat + ( a1 / 2 ); + j4 = ( ( 3 * j2 ) + a2 ) / 4; + j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3; + alfa = ( 3 / 4 ) * e2cuadrada; + beta = ( 5 / 3 ) * alfa ^ 2; + gama = ( 35 / 27 ) * alfa ^ 3; + Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 ); + b = ( Y - Bm ) / v; + Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2; + Eps = a * ( 1 - ( Epsi / 3 ) ); + nab = ( b * ( 1 - Epsi ) ) + lat; + senoheps = ( exp(Eps) - exp(-Eps) ) / 2; + Delt = atan(senoheps / (cos(nab) ) ); + TaO = atan(cos(Delt) * tan(nab)); + longitude = (Delt *(180 / pi ) ) + S; + latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ... + (180 / pi); + + Lat(i)=latitude; + Lon(i)=longitude; + +end Index: /issm/trunk-jpl/src/m/coordsystems/xy2lambert.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/xy2lambert.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/xy2lambert.m (revision 15139) @@ -0,0 +1,69 @@ +function [lat,lon] = xy2lambert(x,y,sgn,projection_center_lat,projection_center_lon) +%XY2LAMBERT - converts xy to lat lon in Lambert Azimuthal +% +% Converts from Ploar Stereographic (X,Y) coordinates to geodetic +% lat lon that are in Lambert Azimuthal (equal area) projection. +% +% Usage: +% [lat,lon] = xy2lambert(x,y,sgn) +% [lat,lon] = xy2lambert(x,y,sgn,projection_center_lat,projection_center_lon) +% +% - provide lat in [-90,90] and lon in [-180,180]. +% +% - sgn = +1 N hemisphere [default projection center lat = 90 lon=0] +% -1 S hemisphere [default projection center lat = -90 lon=0] + +%Get projection_center_lat and projection_center_lon +if nargin==5, + latitude0 = projection_center_lat; + longitude0 = projection_center_lon; +elseif nargin==3, + if sgn==1, + latitude0 = 90; longitude0 = 0; + disp('Info: creating coordinates in Lambert Azimuthal equal-area (Projection center lat: 90N lon: 0)'); + elseif sgn==-1, + latitude0 = -90; longitude0 = 0; + disp('Info: creating coordinates in Lambert Azimuthal equal-area (Projection center lat: 90S lon: 0)'); + else + error('Sign should be either +1 or -1'); + end +else + help xy2lambert + error('bad usage'); +end + +% Radius of the earth in meters +a = 6378137.0; +% Eccentricity of the Hughes ellipsoid squared +e = 0.081819191; + +% Projection center latitude and longitude in radians +phi0 = latitude0 * pi/180; +lam0 = longitude0 * pi/180; + +% Some constants based on phi0 and lam0 +% (as in forward calculation) +qp= (1-e^2)*((1/(1-e^2))-((1/(2*e))*log((1-e)/(1+e)))); +q0=(1-e^2)*((sin(phi0)/(1-e^2*sin(phi0)*sin(phi0)))-((1/(2*e))*log((1-e*sin(phi0))/(1+e*sin(phi0))))); +Rq=a*sqrt(qp/2); +b0=asin(q0/qp); +D =a*(cos(phi0)/sqrt(1-e^2*sin(phi0)*sin(phi0)))/(Rq*cos(b0)); + +% Some other (x,y) dependent parameters +rho=sqrt((x/D)^2+(D*y)^2); +C=2*asin(rho/(2*Rq)); +b_prime=asin((cos(C)*sin(b0))+((D*y*sin(C)*cos(b0))/rho)); + +% Calculation of lat and lon +dist=sqrt(x^2+y^2); +if(dist<=0.1) + lat=sgn*90.0; + lon=0.0; +else + lat_rad=b_prime+((e^2/3+31*e^4/180+517*e^6/5040)*sin(2*b_prime))+((23*e^4/360+251*e^6/3780)*sin(4*b_prime))+((761*e^6/45360)*sin(6*b_prime)); + lon_rad=lam0+atan(x*sin(C)/(D*rho*cos(b0)*cos(C)-D^2*y*sin(b0)*sin(C))); + % in degrees + lat=lat_rad*180/pi; + lon=lon_rad*180/pi; +end + Index: /issm/trunk-jpl/src/m/coordsystems/xy2ll.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/xy2ll.m (revision 15139) +++ /issm/trunk-jpl/src/m/coordsystems/xy2ll.m (revision 15139) @@ -0,0 +1,72 @@ +function [lat,lon] = xy2ll(x,y,sgn,central_meridian,standard_parallel) +%XY2LL - converts xy to lat long +% +% Converts Polar Stereographic (X,Y) coordinates for the polar regions to +% latitude and longitude Stereographic (X,Y) coordinates for the polar +% regions. +% Author: Michael P. Schodlok, December 2003 (map2xy.m) +% +% Usage: +% [lat,lon] = xy2ll(x,y,sgn); +% [lat,lon] = xy2ll(x,y,sgn,central_meridian,standard_parallel); +% +% - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70) +% -1 : south latitude (default is mer=0 lat=71) + +%Get central_meridian and standard_parallel depending on hemisphere +if nargin==5, + delta = central_meridian; + slat = standard_parallel; +elseif nargin==3 + if sgn == 1, + delta = 45; slat = 70; + disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); + elseif sgn==-1, + delta = 0; slat = 71; + disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 71ºS Meridian: 0º)'); + else + error('Sign should be either +1 or -1'); + end +else + help xy2ll + error('bad usage'); +end + +% Conversion constant from degrees to radians +cde = 57.29577951; +% Radius of the earth in meters +re = 6378.273*10^3; +% Eccentricity of the Hughes ellipsoid squared +ex2 = .006693883; +% Eccentricity of the Hughes ellipsoid +ex = sqrt(ex2); + +sl = slat*pi/180.; +rho = sqrt(x.^2 + y.^2); +cm = cos(sl) / sqrt(1.0 - ex2 * (sin(sl)^2)); +T = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl)))^(ex / 2.0); + +if abs(slat-90.) < 1.e-5 + T = rho * sqrt((1. + ex)^(1. + ex) * (1. - ex)^(1. - ex)) / 2. / re; +else + T = rho * T / (re * cm); +end + +chi = (pi / 2.0) - 2.0 * atan(T); +lat = chi + ((ex2 / 2.0) + (5.0 * ex2^2.0 / 24.0) + (ex2^3.0 / 12.0)) * ... + sin(2 * chi) + ((7.0 * ex2^2.0 / 48.0) + (29.0 * ex2^3 / 240.0)) * ... + sin(4.0 * chi) + (7.0 * ex2^3.0 / 120.0) * sin(6.0 * chi) ; + +lat = sgn * lat; +lon = atan2(sgn * x,-sgn * y); +lon = sgn * lon; + +[res1,res2] = find(rho <= 0.1); +if res1 + lat(res1,1) = 90. * sgn; + lon(res1,1) = 0.0; +end + +lon = lon * 180. / pi; +lat = lat * 180. / pi; +lon = lon - delta;