Changeset d0fc203 in tol


Ignore:
Timestamp:
Jan 2, 2010, 12:31:12 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
a8dd8b9
Parents:
59e155b
Message:

New Class @UTMFormula to handle with Universal Transverse Mercator System

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1898 be4846ae-f051-483a-831f-b64e7d08527f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol_pkg/GisTools/GisTools.tol

    r59e155b rd0fc203  
    1818Set _.autodoc.authors = [[
    1919   "vdebuen@tol-project.org"   ]];
     20
     21
     22
     23//////////////////////////////////////////////////////////////////////////////
     24//Abstract class to handle with Universal Transverse Mercator System
     25//http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM
     26Class @UTMFormula
     27//////////////////////////////////////////////////////////////////////////////
     28{
     29  //Convertingfactor from decimal degrees to radians
     30  Static Real deg2rad = Pi/180;
     31
     32  ////////////////////////////////////////////////////////////////////////////
     33  //UTM Datum defining
     34  ////////////////////////////////////////////////////////////////////////////
     35  //Datum name
     36  Text _.name;
     37  //Equatorial Radius, meters
     38  Real _.a;
     39  //Polar Radius, meters
     40  Real _.b;
     41  //Geographic zone of application
     42  Text _.use;
     43
     44  //scale along central meridian of zone. Even though itŽs a constant,
     45  //we retain it as a separate symbol to keep the numerical
     46  //coefficients simpler, also to allow for systems that might use
     47  //a different Mercator projection.
     48  Real _.k0 = 0.9996;
     49  //eccentricity of the earthŽs elliptical cross-section.
     50  Real _.e = ?;
     51  //auxiliar transformed eccentricity.
     52  Real _.e.2 = ?;
     53  Real _.e_1 = ?;
     54  //auxiliar ratio of difference and sum of elliptic radius
     55  Real _.n = ?;
     56  //auxiliar variables used to calculate Meridional Arc
     57  Real _.A = ?;
     58  Real _.B = ?;
     59  Real _.C = ?;
     60  Real _.D = ?;
     61  Real _.E = ?;
     62  //auxiliar variables used to calculate Footprint Latitude
     63  Real _.J1 = ?;
     64  Real _.J2 = ?;
     65  Real _.J3 = ?;
     66  Real _.J4 = ?;
     67
     68  ////////////////////////////////////////////////////////////////////////////
     69  Real initialice(Real void)
     70  ////////////////////////////////////////////////////////////////////////////
     71  {
     72    Real e2 = 1-(_.b/_.a)^2;
     73    Real _.e := Sqrt(e2);
     74    Real e4 = _.e^4;
     75    Real e6 = _.e^6;
     76    Real _.e.2 := e2/(1-e2);
     77    Real _.n := (_.a-_.b)/(_.a+_.b);
     78    Real _.e_1 := (1-Sqrt((1-e2)))/(1+Sqrt((1-e2)));
     79    //Calculating the arc length of an ellipse involves functions called
     80    //elliptic integrals, which donŽt reduce to neat closed formulas.
     81    //So they have to be represented as series.
     82/* */
     83    //Army form
     84    Real n23 = _.n^2-_.n^3;
     85    Real n45 = _.n^4-_.n^5;
     86    Real _.A := +_.a*                  (1-_.n + (5/4)*n23 + (81/64)*n45);
     87    Real _.B := -_.a*((3/2)    *_.n*   (1-_.n + (7/8)*n23 + (55/64)*n45));
     88    Real _.C := +_.a*((15/16)  *_.n^2* (1-_.n + (3/4)*n23));
     89    Real _.D := -_.a*((35/48)  *_.n^3* (1-_.n + (11/16)*n23));
     90    Real _.E := +_.a*((315/512)*_.n^4* (1-_.n));
     91/* * /
     92    //The USGS gives this form, which may be more appealing to some.
     93    Real _.A := +_.a*(1-  e2/4 -3*e4/64  - 5*e6/256 );
     94    Real _.B := -_.a*(  3*e2/8 +3*e4/32  +45*e6/1024);
     95    Real _.C := +_.a*(         15*e4/256 +45*e6/1024);
     96    Real _.D := -_.a*(                    35*e6/3072);
     97    Real _.E := 0;
     98/* */
     99    //Auxiliar variables to calculate Footprint Latitude
     100    Real E1 = _.e_1;
     101    Real E2 = _.e_1^2;
     102    Real E3 = _.e_1^3;
     103    Real E4 = _.e_1^4;
     104    Real _.J1 := 3*E1/2 - 27*E3/32;
     105    Real _.J2 := 21*E2/16 - 55*E4/32;
     106    Real _.J3 := 151*E3/96;
     107    Real _.J4 := 1097*E4/512;
     108
     109    True
     110  };
     111  ////////////////////////////////////////////////////////////////////////////
     112  Static @UTMFormula New(Text name, Real a, Real b, Text use)
     113  ////////////////////////////////////////////////////////////////////////////
     114  {
     115    @UTMFormula aux = [[
     116      Text _.name = name;
     117      Text _.use = use;
     118      Real _.a = a;
     119      Real _.b = b
     120    ]];
     121    Real aux::initialice(?);
     122    Eval(name+"=aux")
     123  };
     124
     125  ////////////////////////////////////////////////////////////////////////////
     126  Static Set datum = Copy(Empty);
     127  Static Real BuildDatum(Real void)
     128  ////////////////////////////////////////////////////////////////////////////
     129  {
     130    Set datum := [[
     131    @UTMFormula::New("NAD83",              6378137,   6356752.3142, "Global");
     132    @UTMFormula::New("WGS84",              6378137,   6356752.3142, "Global");
     133    @UTMFormula::New("GRS80",              6378137,   6356752.3141, "US");
     134    @UTMFormula::New("WGS72",              6378135,   6356750.5,    "NASA, DOD");
     135    @UTMFormula::New("Australian.1965",    6378160,   6356774.7,    "Australia");
     136    @UTMFormula::New("Krasovsky.1940",     6378245,   6356863.0,    "Soviet Union");
     137    @UTMFormula::New("International.1924", 6378388,   6356911.9,    "Global except as listed");
     138    @UTMFormula::New("Hayford.1909",       6378388,   6356911.9,    "Global except as listed");
     139    @UTMFormula::New("Clake.1880",         6378249.1, 6356514.9,    "France, Africa");
     140    @UTMFormula::New("Clarke.1866",        6378206.4, 6356583.8,    "North America");
     141    @UTMFormula::New("Airy.1830",          6377563.4, 6356256.9,    "Great Britain");
     142    @UTMFormula::New("Bessel.1841",        6377397.2, 6356079.0,    "Central Europe, Chile, Indonesia");
     143    @UTMFormula::New("Everest.1830",       6377276.3, 6356075.4,    "South Asia") ]];
     144    Real SetIndexByName(datum);
     145    True
     146  };
     147
     148
     149  ////////////////////////////////////////////////////////////////////////////
     150  Static Text _.autodoc.member.get.zone.deg= "Returns UTM zone of longitude "
     151  "in degrees";
     152  Static VMatrix get.zone.deg(Real long)
     153  ////////////////////////////////////////////////////////////////////////////
     154  {
     155    If(long<0,Floor((180+long)/6)+1,Floor(long/6)+31)
     156  };
     157
     158  ////////////////////////////////////////////////////////////////////////////
     159  Static Text _.autodoc.member.get.zone.rad = "Returns UTM zone of longitude "
     160  "in radians";
     161  Static VMatrix get.zone.rad(Real long)
     162  ////////////////////////////////////////////////////////////////////////////
     163  {
     164    @UTMFormula::get.zone.deg(long/deg2rad)
     165  };
     166
     167  ////////////////////////////////////////////////////////////////////////////
     168  Static Text _.autodoc.member.get.zoneCentralMeridian.deg=
     169  "Returns central meridian of UTM zone for given longitude in degrees";
     170  Static VMatrix get.zoneCentralMeridian.deg(Real long)
     171  ////////////////////////////////////////////////////////////////////////////
     172  {
     173    6*@UTMFormula::get.zone.deg(long)-183
     174  };
     175
     176  ////////////////////////////////////////////////////////////////////////////
     177  Static Text _.autodoc.member.get.zoneCentralMeridian.rad =
     178  "Returns entral meridian of UTM zone for given longitude in radians";
     179  Static VMatrix get.zoneCentralMeridian.rad(Real long)
     180  ////////////////////////////////////////////////////////////////////////////
     181  {
     182    @UTMFormula::get.zoneCentralMeridian.deg(long/deg2rad)
     183  };
     184
     185  ////////////////////////////////////////////////////////////////////////////
     186  Static Text _.autodoc.member.rad2utm = "Formulas for converting "
     187  "coordinates of latitude and longitude in radians to UTM";
     188  VMatrix rad2utm(
     189    VMatrix lonlat, //Matrix with two columns, first one is longitud, second
     190                    //one is latitude, both meassured in in radians
     191    Real zone) //UTM zone
     192  ////////////////////////////////////////////////////////////////////////////
     193  {
     194    //Central meridian of zone
     195    Real long0 = (6*zone-183)*deg2rad;
     196    //longitude of points
     197    VMatrix long = ASin(Sin(SubCol(lonlat,[[1]])));
     198    //latitude of points
     199    VMatrix lat = ASin(Sin(SubCol(lonlat,[[2]])));
     200    Real e2 = _.e^2;
     201    VMatrix e2sinlat2_1 = -(Sin(lat)^2)*e2+1;
     202    //This is the radius of curvature  of the earth in the meridian plane.
     203    VMatrix rho = (e2sinlat2_1^(-3/2))*(_.a*(1-e2));
     204    //This is the radius of curvature of the earth perpendicular to the
     205    //meridian plane. It is also the distance from the point in question
     206    //to the polar axis, measured perpendicular to the earthŽs surface.
     207    VMatrix nu = (e2sinlat2_1^(-1/2))*_.a;
     208    VMatrix p = long-long0;
     209    VMatrix k0nu = nu*_.k0;
     210    VMatrix sinlat  = Sin(lat);
     211    VMatrix sin2lat = Sin(lat*2);
     212    VMatrix sin4lat = Sin(lat*4);
     213    VMatrix sin6lat = Sin(lat*6);
     214    VMatrix sin8lat = Sin(lat*8);
     215    VMatrix coslat  = Cos(lat);
     216    VMatrix coslat2  = coslat^2;
     217    VMatrix coslat3  = coslat^3;
     218    VMatrix coslat4  = coslat^4;
     219    VMatrix tanlat   = Tan(lat);
     220    VMatrix tanlat2  = tanlat^2;
     221    VMatrix tanlat3  = tanlat^3;
     222    VMatrix S =
     223      lat*_.A +
     224      sin2lat*_.B +
     225      sin4lat*_.C +
     226      sin6lat*_.D +
     227      sin8lat*_.E;
     228    VMatrix K1 = S*_.k0;
     229    VMatrix K2 = sin2lat$*k0nu/4;
     230    VMatrix K3 =
     231     (sinlat$*coslat3$*(k0nu/24))$*
     232     (-tanlat2 + coslat2*(9*_.e.2) + coslat4*(4*_.e.2^2) + 5);
     233    VMatrix K4 = coslat$*k0nu;
     234    VMatrix K5 = (coslat3$*k0nu/6)$*
     235                 (-tanlat2 + coslat2*_.e.2 +1);
     236    //northing of points
     237    VMatrix y = (K1 + K2$*p^2 + K3$*p^4);
     238    //easting of points
     239    VMatrix x = (K4$*p + K5$*p^3)+500000;
     240    x|IfVMat(LT(lat,lat*0),y+1.E+7,y)
     241  };
     242
     243  ////////////////////////////////////////////////////////////////////////////
     244  Static Text _.autodoc.member.deg2utm = "Formulas for converting "
     245  "coordinates of latitude and longitude in degrees to UTM";
     246  VMatrix deg2utm(
     247    VMatrix lonlat, //Matrix with two columns, first one is longitud, second
     248                    //one is latitude, both meassured in degrees
     249    Real zone) //UTM zone
     250  ////////////////////////////////////////////////////////////////////////////
     251  {
     252    rad2utm(lonlat*deg2rad,zone)
     253  };
     254
     255  ////////////////////////////////////////////////////////////////////////////
     256  Static Text _.autodoc.member.utm2rad = "Formulas for converting "
     257  "UTM coordinates to latitude and longitude in radians";
     258  VMatrix utm2rad(
     259    VMatrix xy, //Matrix with two columns
     260    Real zone) //UTM zone
     261  ////////////////////////////////////////////////////////////////////////////
     262  {
     263    //Central meridian of zone
     264    Real long0 = (6*zone-183)*deg2rad;
     265    //easting of points
     266    VMatrix x = SubCol(xy,[[1]])-500000;
     267    //northing of points
     268    VMatrix y = SubCol(xy,[[2]]);
     269    Real e2 = _.e^2;
     270    Real e4 = _.e^4;
     271    Real e6 = _.e^6;
     272  //WriteLn("TRACE x="<<VMat2Mat(x));
     273  //WriteLn("TRACE y="<<VMat2Mat(y));
     274    //Calculate the Meridional Arc
     275    VMatrix M = y/_.k0;
     276  //WriteLn("TRACE M="<<VMat2Mat(M));
     277    //Calculating Footprint Latitude
     278    VMatrix mu = M/(_.a*(1-  e2/4 -3*e4/64  - 5*e6/256));
     279  //WriteLn("TRACE mu="<<VMat2Mat(mu));
     280    VMatrix fp =
     281      mu +
     282      Sin(mu*2)*_.J1 +
     283      Sin(mu*4)*_.J2 +
     284      Sin(mu*6)*_.J3 +
     285      Sin(mu*8)*_.J4;
     286  //WriteLn("TRACE fp="<<VMat2Mat(fp));
     287    VMatrix sinfp = Sin(fp);
     288    VMatrix cosfp = Cos(fp);
     289    VMatrix tanfp = sinfp$/cosfp;
     290    VMatrix C1 = (cosfp^2)*_.e.2;
     291  //WriteLn("TRACE C1="<<VMat2Mat(C1));
     292    VMatrix T1 = tanfp^2;
     293  //WriteLn("TRACE T1="<<VMat2Mat(T1));
     294    VMatrix sinfp2e21 = -(sinfp^2)*_.e^2+1;
     295    VMatrix R1 = sinfp2e21^(-3/2)*(_.a*(1-_.e^2));
     296  //WriteLn("TRACE R1="<<VMat2Mat(R1));
     297    VMatrix N1 = sinfp2e21^(-1/2)*_.a;
     298  //WriteLn("TRACE N1="<<VMat2Mat(N1));
     299    VMatrix D = x$/(N1*_.k0);
     300  //WriteLn("TRACE D="<<VMat2Mat(D));
     301    VMatrix Q1 = N1$*tanfp$/R1;
     302  //WriteLn("TRACE Q1="<<VMat2Mat(Q1));
     303    VMatrix Q2 = (D^2)/2;
     304  //WriteLn("TRACE Q2="<<VMat2Mat(Q2));
     305    VMatrix Q3 = (T1* 3 + C1* 10           - C1^2*4 +(5-   9*_.e.2))$*D^4/24;
     306  //WriteLn("TRACE Q3="<<VMat2Mat(Q3));
     307    VMatrix Q4 = (T1*90 + C1*298 + T1^2*45 - C1^2*3 +(61-252*_.e.2))$*D^6/720;
     308  //WriteLn("TRACE Q4="<<VMat2Mat(Q4));
     309    VMatrix Q5 = D;
     310  //WriteLn("TRACE Q5="<<VMat2Mat(Q5));
     311    VMatrix Q6 = (T1*2  + C1 + 1)$*D^3/6;
     312  //WriteLn("TRACE Q6="<<VMat2Mat(Q6));
     313    VMatrix Q7 = (T1*28 - C1*2   + T1^2*24 - C1^2*3 +(5 +  8*_.e.2))$*D^5/120;
     314  //WriteLn("TRACE Q7="<<VMat2Mat(Q7));
     315    VMatrix lat = fp - Q1$*(Q2 - Q3 + Q4);
     316  //WriteLn("TRACE lat="<<VMat2Mat(lat));
     317    VMatrix long = (Q5 - Q6 + Q7)$/cosfp+long0;
     318  //WriteLn("TRACE long="<<VMat2Mat(long));
     319    long | lat
     320  };
     321  ////////////////////////////////////////////////////////////////////////////
     322  Static Text _.autodoc.member.utm2deg = "Formulas for converting "
     323  "UTM coordinates to latitude and longitude in degrees";
     324  VMatrix utm2deg(
     325    VMatrix xy, //Matrix with two columns
     326    Real zone) //UTM zone
     327  ////////////////////////////////////////////////////////////////////////////
     328  {
     329    utm2rad(xy,zone)/deg2rad
     330  }
     331};
     332
     333
     334Real _UTMFormula.BuildDatum = @UTMFormula::BuildDatum(?);
    20335
    21336//////////////////////////////////////////////////////////////////////////////
    … …  
    127442};
    128443
     444
    129445//////////////////////////////////////////////////////////////////////////////
    130446//Abstract class to inherite bidimensional interpolating classes.
    … …  
    162478  ////////////////////////////////////////////////////////////////////////////
    163479  Static Text _.autodoc.member.eval =
    164   "Evaluates the first interpolation function at point x,y. It's useful in "
     480  "Evaluates the first interpolation function at point x,y. ItŽs useful in "
    165481  "the scalar case, when _.F has just one column (_.n==1)";
    166482  Real eval(@Point2D p)
    … …  
    326642
    327643
     644
Note: See TracChangeset for help on using the changeset viewer.