Changeset c46a7b2 in tol


Ignore:
Timestamp:
Jan 5, 2010, 7:18:01 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
a693055
Parents:
36977fc
Message:

División temática de ficheros

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

Location:
tol_pkg/GisTools
Files:
4 added
1 edited

Legend:

Unmodified
Added
Removed

  • tol_pkg/GisTools/GisTools.tol

    r36977fc rc46a7b2  
    1919   "vdebuen@tol-project.org"   ]];
    2020
    21 
    22 
    23 //////////////////////////////////////////////////////////////////////////////
    24 //Abstract class to handle with Universal Transverse Mercator System
    25 //http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM
    26 Class @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 initialize(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::initialize(?);
    122     Eval(name+"=aux")
    123   };
    124 
    125   ////////////////////////////////////////////////////////////////////////////
    126   Static Set _datum = Copy(Empty);
    127   Static Set Datum(Real void)
    128   ////////////////////////////////////////////////////////////////////////////
    129   {
    130     If(Card(_datum),_datum, {
    131     Set _datum := [[
    132     @UTMFormula::New("NAD83",              6378137,   6356752.3142, "Global");
    133     @UTMFormula::New("WGS84",              6378137,   6356752.3142, "Global");
    134     @UTMFormula::New("GRS80",              6378137,   6356752.3141, "US");
    135     @UTMFormula::New("WGS72",              6378135,   6356750.5,    "NASA, DOD");
    136     @UTMFormula::New("Australian.1965",    6378160,   6356774.7,    "Australia");
    137     @UTMFormula::New("Krasovsky.1940",     6378245,   6356863.0,    "Soviet Union");
    138     @UTMFormula::New("International.1924", 6378388,   6356911.9,    "Global except as listed");
    139     @UTMFormula::New("Hayford.1909",       6378388,   6356911.9,    "Global except as listed");
    140     @UTMFormula::New("Clake.1880",         6378249.1, 6356514.9,    "France, Africa");
    141     @UTMFormula::New("Clarke.1866",        6378206.4, 6356583.8,    "North America");
    142     @UTMFormula::New("Airy.1830",          6377563.4, 6356256.9,    "Great Britain");
    143     @UTMFormula::New("Bessel.1841",        6377397.2, 6356079.0,    "Central Europe, Chile, Indonesia");
    144     @UTMFormula::New("Everest.1830",       6377276.3, 6356075.4,    "South Asia") ]];
    145     Real SetIndexByName(_datum);
    146     _datum})
    147   };
    148 
    149   ////////////////////////////////////////////////////////////////////////////
    150   Static @UTMFormula GetDatum(Text name)
    151   ////////////////////////////////////////////////////////////////////////////
    152   {
    153     Set dtm = Datum(0);
    154     dtm[name]
    155   };
    156 
    157   ////////////////////////////////////////////////////////////////////////////
    158   Static Text _.autodoc.member.get.zone.deg= "Returns UTM zone of longitude "
    159   "in degrees";
    160   Static VMatrix get.zone.deg(Real long)
    161   ////////////////////////////////////////////////////////////////////////////
    162   {
    163     If(long<0,Floor((180+long)/6)+1,Floor(long/6)+31)
    164   };
    165 
    166   ////////////////////////////////////////////////////////////////////////////
    167   Static Text _.autodoc.member.get.zone.rad = "Returns UTM zone of longitude "
    168   "in radians";
    169   Static VMatrix get.zone.rad(Real long)
    170   ////////////////////////////////////////////////////////////////////////////
    171   {
    172     @UTMFormula::get.zone.deg(long/deg2rad)
    173   };
    174 
    175   ////////////////////////////////////////////////////////////////////////////
    176   Static Text _.autodoc.member.get.zoneCentralMeridian.deg=
    177   "Returns central meridian of UTM zone for given longitude in degrees";
    178   Static VMatrix get.zoneCentralMeridian.deg(Real long)
    179   ////////////////////////////////////////////////////////////////////////////
    180   {
    181     6*@UTMFormula::get.zone.deg(long)-183
    182   };
    183 
    184   ////////////////////////////////////////////////////////////////////////////
    185   Static Text _.autodoc.member.get.zoneCentralMeridian.rad =
    186   "Returns entral meridian of UTM zone for given longitude in radians";
    187   Static VMatrix get.zoneCentralMeridian.rad(Real long)
    188   ////////////////////////////////////////////////////////////////////////////
    189   {
    190     @UTMFormula::get.zoneCentralMeridian.deg(long/deg2rad)
    191   };
    192 
    193   ////////////////////////////////////////////////////////////////////////////
    194   Static Text _.autodoc.member.rad2utm = "Formulas for converting "
    195   "coordinates of latitude and longitude in radians to UTM";
    196   VMatrix rad2utm(
    197     VMatrix lonlat, //Matrix with two columns, first one is longitud, second
    198                     //one is latitude, both meassured in in radians
    199     VMatrix lonZone) //UTM longitud zone
    200   ////////////////////////////////////////////////////////////////////////////
    201   {
    202     //Central meridian of zone
    203     VMatrix long0 = (lonZone*6-183)*deg2rad;
    204     //longitude of points
    205     VMatrix long = SubCol(lonlat,[[1]]);
    206     //latitude of points
    207     VMatrix lat = SubCol(lonlat,[[2]]);
    208     Real e2 = _.e^2;
    209     VMatrix e2sinlat2_1 = -(Sin(lat)^2)*e2+1;
    210     //This is the radius of curvature  of the earth in the meridian plane.
    211     VMatrix rho = (e2sinlat2_1^(-3/2))*(_.a*(1-e2));
    212     //This is the radius of curvature of the earth perpendicular to the
    213     //meridian plane. It is also the distance from the point in question
    214     //to the polar axis, measured perpendicular to the earthŽs surface.
    215     VMatrix nu = (e2sinlat2_1^(-1/2))*_.a;
    216     VMatrix p = long-long0;
    217     VMatrix k0nu = nu*_.k0;
    218     VMatrix sinlat  = Sin(lat);
    219     VMatrix sin2lat = Sin(lat*2);
    220     VMatrix sin4lat = Sin(lat*4);
    221     VMatrix sin6lat = Sin(lat*6);
    222     VMatrix sin8lat = Sin(lat*8);
    223     VMatrix coslat  = Cos(lat);
    224     VMatrix coslat2  = coslat^2;
    225     VMatrix coslat3  = coslat^3;
    226     VMatrix coslat4  = coslat^4;
    227     VMatrix tanlat   = Tan(lat);
    228     VMatrix tanlat2  = tanlat^2;
    229     VMatrix tanlat3  = tanlat^3;
    230     VMatrix S =
    231       lat*_.A +
    232       sin2lat*_.B +
    233       sin4lat*_.C +
    234       sin6lat*_.D +
    235       sin8lat*_.E;
    236     VMatrix K1 = S*_.k0;
    237     VMatrix K2 = sin2lat$*k0nu/4;
    238     VMatrix K3 =
    239      (sinlat$*coslat3$*(k0nu/24))$*
    240      (-tanlat2 + coslat2*(9*_.e.2) + coslat4*(4*_.e.2^2) + 5);
    241     VMatrix K4 = coslat$*k0nu;
    242     VMatrix K5 = (coslat3$*k0nu/6)$*
    243                  (-tanlat2 + coslat2*_.e.2 +1);
    244     //northing of points
    245     VMatrix y = (K1 + K2$*p^2 + K3$*p^4);
    246     //easting of points
    247     VMatrix x = (K4$*p + K5$*p^3)+500000;
    248     x|IfVMat(LT(lat,lat*0),y+1.E+7,y)
    249   };
    250 
    251   ////////////////////////////////////////////////////////////////////////////
    252   Static Text _.autodoc.member.deg2utm = "Formulas for converting "
    253   "coordinates of latitude and longitude in degrees to UTM";
    254   VMatrix deg2utm(
    255     VMatrix lonlat, //Matrix with two columns, first one is longitud, second
    256                     //one is latitude, both meassured in degrees
    257     VMatrix lonZone) //UTM longitud zone
    258   ////////////////////////////////////////////////////////////////////////////
    259   {
    260     rad2utm(lonlat*deg2rad,lonZone)
    261   };
    262 
    263   ////////////////////////////////////////////////////////////////////////////
    264   Static Text _.autodoc.member.utm2rad = "Formulas for converting "
    265   "UTM coordinates to latitude and longitude in radians";
    266   VMatrix utm2rad(
    267     VMatrix xy, //Matrix with two columns
    268     VMatrix lonZone, //UTM longitud zone
    269     VMatrix latZone) //UTM latitud zone
    270   ////////////////////////////////////////////////////////////////////////////
    271   {
    272   //WriteLn("TRACE zone="<<zone);
    273     //Central meridian of zone
    274     VMatrix long0 = (lonZone*6-183)*deg2rad;
    275   //WriteLn("TRACE long0="<<long0);
    276     //easting of points
    277     VMatrix x = SubCol(xy,[[1]])-500000;
    278     //northing of points
    279     VMatrix isSouthern = LT(latZone,latZone*0);
    280     VMatrix y =
    281     {
    282       VMatrix aux = SubCol(xy,[[2]]);
    283       IfVMat(isSouthern,-aux+1.E+7,aux)
    284     };
    285     Real e2 = _.e^2;
    286     Real e4 = _.e^4;
    287     Real e6 = _.e^6;
    288   //WriteLn("TRACE x="<<VMat2Mat(x));
    289   //WriteLn("TRACE y="<<VMat2Mat(y));
    290     //Calculate the Meridional Arc
    291     VMatrix M = y/_.k0;
    292   //WriteLn("TRACE M="<<VMat2Mat(M));
    293     //Calculating Footprint Latitude
    294     VMatrix mu = M/(_.a*(1-  e2/4 -3*e4/64  - 5*e6/256));
    295   //WriteLn("TRACE mu="<<VMat2Mat(mu));
    296     VMatrix fp =
    297       mu +
    298       Sin(mu*2)*_.J1 +
    299       Sin(mu*4)*_.J2 +
    300       Sin(mu*6)*_.J3 +
    301       Sin(mu*8)*_.J4;
    302   //WriteLn("TRACE fp="<<VMat2Mat(fp));
    303     VMatrix sinfp = Sin(fp);
    304     VMatrix cosfp = Cos(fp);
    305     VMatrix tanfp = Tan(fp);
    306     VMatrix C1 = (cosfp^2)*_.e.2;
    307   //WriteLn("TRACE C1="<<VMat2Mat(C1));
    308     VMatrix T1 = tanfp^2;
    309   //WriteLn("TRACE T1="<<VMat2Mat(T1));
    310     VMatrix sinfp2e21 = -(sinfp^2)*_.e^2+1;
    311     VMatrix R1 = sinfp2e21^(-3/2)*(_.a*(1-_.e^2));
    312   //WriteLn("TRACE R1="<<VMat2Mat(R1));
    313     VMatrix N1 = sinfp2e21^(-1/2)*_.a;
    314   //WriteLn("TRACE N1="<<VMat2Mat(N1));
    315     VMatrix D = x$/(N1*_.k0);
    316   //WriteLn("TRACE D="<<VMat2Mat(D));
    317     VMatrix Q1 = N1$*tanfp$/R1;
    318   //WriteLn("TRACE Q1="<<VMat2Mat(Q1));
    319     VMatrix Q2 = (D^2)/2;
    320   //WriteLn("TRACE Q2="<<VMat2Mat(Q2));
    321     VMatrix Q3 = (T1* 3 + C1* 10           - C1^2*4 +(5-   9*_.e.2))$*D^4/24;
    322   //WriteLn("TRACE Q3="<<VMat2Mat(Q3));
    323     VMatrix Q4 = (T1*90 + C1*298 + T1^2*45 - C1^2*3 +(61-252*_.e.2))$*D^6/720;
    324   //WriteLn("TRACE Q4="<<VMat2Mat(Q4));
    325     VMatrix Q5 = D;
    326   //WriteLn("TRACE Q5="<<VMat2Mat(Q5));
    327     VMatrix Q6 = (T1*2  + C1 + 1)$*D^3/6;
    328   //WriteLn("TRACE Q6="<<VMat2Mat(Q6));
    329     VMatrix Q7 = (T1*28 - C1*2   + T1^2*24 - C1^2*3 +(5 +  8*_.e.2))$*D^5/120;
    330   //WriteLn("TRACE Q7="<<VMat2Mat(Q7));
    331     VMatrix lat = fp - Q1$*(Q2 - Q3 + Q4);
    332   //WriteLn("TRACE lat="<<VMat2Mat(lat));
    333     VMatrix long = (Q5 - Q6 + Q7)$/cosfp+long0;
    334   //WriteLn("TRACE long="<<VMat2Mat(long));
    335     long | IfVMat(isSouthern,-lat,lat)
    336   };
    337   ////////////////////////////////////////////////////////////////////////////
    338   Static Text _.autodoc.member.utm2deg = "Formulas for converting "
    339   "UTM coordinates to latitude and longitude in degrees";
    340   VMatrix utm2deg(
    341     VMatrix xy, //Matrix with two columns
    342     VMatrix lonZone, //UTM longitud zone
    343     VMatrix latZone) //UTM latitud zone
    344   ////////////////////////////////////////////////////////////////////////////
    345   {
    346     utm2rad(xy,lonZone,latZone)/deg2rad
    347   }
    348 };
    349 
    350 //////////////////////////////////////////////////////////////////////////////
    351 //Abstract class to handle with a bidimensional point.
    352 Class @Point2D
    353 //////////////////////////////////////////////////////////////////////////////
    354 {
    355   ////////////////////////////////////////////////////////////////////////////
    356   Static Text _.autodoc.member.numberOfPoints =
    357   "Returns the number of stored points";
    358   Real numberOfPoints(Real void);
    359   ////////////////////////////////////////////////////////////////////////////
    360 
    361   ////////////////////////////////////////////////////////////////////////////
    362   Static Text _.autodoc.member.x =
    363   "Returns the first coordinate of the k-th point";
    364   Real x(Real k);
    365   ////////////////////////////////////////////////////////////////////////////
    366 
    367   ////////////////////////////////////////////////////////////////////////////
    368   Static Text _.autodoc.member.y =
    369   "Returns the second coordinate of the k-th point";
    370   Real y(Real k)
    371   ////////////////////////////////////////////////////////////////////////////
    372 };
    373 
    374 //////////////////////////////////////////////////////////////////////////////
    375 //Class to handle with just one bidimensional point that
    376 //isnŽt externally stored
    377 Class @Point2D.Simple : @Point2D
    378 //////////////////////////////////////////////////////////////////////////////
    379 {
    380   Real _.x;
    381   Real _.y;
    382 
    383   //Just one point is stored
    384   Real numberOfPoints(Real void) { 1 };
    385   //Argument k is unused
    386   Real x(Real k) { _.x };
    387   //Argument k is unused
    388   Real y(Real k) { _.y };
    389   ////////////////////////////////////////////////////////////////////////////
    390   Static Text _.autodoc.member.New = "Creates an instance of @Point2D.Simple";
    391   Static @Point2D.Simple New(Real x, Real y)
    392   ////////////////////////////////////////////////////////////////////////////
    393   {
    394     @Point2D.Simple new = [[
    395       Real _.x = x;
    396       Real _.y = y ]];
    397     new
    398   };
    399   ////////////////////////////////////////////////////////////////////////////
    400   Static Text _.autodoc.member.FromRow = "Creates an instance of "
    401   "@Point2D.Simple using the 2 first columns of selected row of given matrix";
    402   Static @Point2D.Simple FromRow(VMatrix M, Real numRow)
    403   ////////////////////////////////////////////////////////////////////////////
    404   {
    405     @Point2D.Simple new = [[
    406       Real _.x = VMatDat(M,numRow,1);
    407       Real _.y = VMatDat(M,numRow,2) ]];
    408     new
    409   };
    410   ////////////////////////////////////////////////////////////////////////////
    411   Static Text _.autodoc.member.FromCol = "Creates an instance of "
    412   "@Point2D.Simple using the 2 first rows of selected column of given matrix";
    413   Static @Point2D.Simple FromCol(VMatrix M, Real numCol)
    414   ////////////////////////////////////////////////////////////////////////////
    415   {
    416     @Point2D.Simple new = [[
    417       Real _.x = VMatDat(M,1,numCol);
    418       Real _.y = VMatDat(M,2,numRow) ]];
    419     new
    420   }
    421 };
    422 
    423 //////////////////////////////////////////////////////////////////////////////
    424 //Class to store al large set of points as rows of a matrix with 2 columns.
    425 Class @Point2D.Stored
    426 //////////////////////////////////////////////////////////////////////////////
    427 {
    428   VMatrix _.xy;
    429 
    430   Real numberOfPoints(Real void) { VRows(_.xy) };
    431   Real x(Real k) { VMatDat(_.xy,k,1) };
    432   Real y(Real k) { VMatDat(_.xy,k,1) };
    433 
    434   @Point2D.Simple center(Real void)
    435   {
    436     Real m = VRows(_.xy);
    437     VMatrix c = Rand(m,1,1/m,1/m)*_.xy;
    438     @Point2D.Simple::New(VMatDat(c,1,1), VMatDat(c,1,2))
    439   };
    440 
    441   ////////////////////////////////////////////////////////////////////////////
    442   Static Text _.autodoc.member.New = "Creates an instance of @Point2D.Stored";
    443   Static @Point2D.Stored New(VMatrix xy)
    444   ////////////////////////////////////////////////////////////////////////////
    445   {
    446     If(VColumns(xy)>2,
    447       WriteLn("[@Point2D.Stored::New] Only 2 first of "<<VColumns(xy)+
    448       " columns of argument xy will be used","W"));
    449     If(VColumns(xy)<2,
    450       WriteLn("[@Point2D.Stored::New] Almost 2 columns instead of "<<
    451       VColumns(xy)+" are needed for argument xy","E"));
    452     @Point2D.Stored new = [[
    453       VMatrix _.xy = xy ]];
    454     new
    455   }
    456 };
    457 
    458 
    459 //////////////////////////////////////////////////////////////////////////////
    460 //Abstract class to inherite bidimensional interpolating classes.
    461 Class @Interp2D
    462 //////////////////////////////////////////////////////////////////////////////
    463 {
    464   //Bidimensional points in cartesian coordinates. Dim: _.m x 2
    465   @Point2D.Stored _.P;
    466   //Tarjet data to be interpolated, Dim : _.m x _.n
    467   VMatrix _.F;
    468 
    469   ////////////////////////////////////////////////////////////////////////////
    470   //Auxiliar members
    471   ////////////////////////////////////////////////////////////////////////////
    472   //Number of interpolation points
    473   Real _.m = ?;
    474   //Number of independent interpolation functions.
    475   //When _.n==1 itŽs called the scalar interpolating case
    476   //When _.n>=2 itŽs called the vectorial interpolating case
    477   Real _.n = ?;
    478 
    479   ////////////////////////////////////////////////////////////////////////////
    480   Static Text _.autodoc.member.build =
    481   "Initializes all members which are needed to build the interpolating "
    482   "methods";
    483   Real build(Real void);
    484   ////////////////////////////////////////////////////////////////////////////
    485 
    486   ////////////////////////////////////////////////////////////////////////////
    487   Static Text _.autodoc.member.eval.vector =
    488   "Evaluates all interpolation function at point x,y";
    489   VMatrix eval.vector(@Point2D p);
    490   ////////////////////////////////////////////////////////////////////////////
    491 
    492   ////////////////////////////////////////////////////////////////////////////
    493   Static Text _.autodoc.member.eval =
    494   "Evaluates the first interpolation function at point x,y. ItŽs useful in "
    495   "the scalar case, when _.F has just one column (_.n==1)";
    496   Real eval(@Point2D p)
    497   ////////////////////////////////////////////////////////////////////////////
    498   {
    499     VMatDat(eval.vector(p),1,1)
    500   }
    501 };
    502 
    503 //////////////////////////////////////////////////////////////////////////////
    504 //Abstract class to inherite bidimensional weighted interpolating classes.
    505 Class @Interp2D.Gravit : @Interp2D
    506 //////////////////////////////////////////////////////////////////////////////
    507 {
    508   //In order to be more efficient, the points farther of this distance will
    509   //be skipped due its weight is almost null
    510   Real _.maxDistance;
    511 
    512   ////////////////////////////////////////////////////////////////////////////
    513   //Auxiliar members
    514   ////////////////////////////////////////////////////////////////////////////
    515   Set     _.point.idx = Copy(Empty);
    516   Set     _.x.range   = Copy(Empty);
    517   VMatrix _.x.class   = Rand(0,0,0,0);
    518   Set     _.y.range   = Copy(Empty);
    519   VMatrix _.y.class   = Rand(0,0,0,0);
    520   VMatrix _.grid.map  = Rand(0,0,0,0);
    521   Set     _.grid      = Copy(Empty);
    522  
    523   ////////////////////////////////////////////////////////////////////////////
    524   Real build(Real void)
    525   ////////////////////////////////////////////////////////////////////////////
    526   {
    527     Real _.m := VRows(_.P::_.xy);
    528     Real _.n := VRows(_.F);
    529     If(VColumns(_.P::_.xy)!=2,
    530       WriteLn("Interp2D.Gravit::_.P::_.xy must have 2 columns instead "
    531               "of "<<VColumns(_.P::_.xy)));
    532     If(VRows(_.F)!=_.m,
    533       WriteLn("Interp2D.Gravit::_.xy must have 2 columns instead of "
    534               <<VColumns(_.P::_.xy)));
    535 
    536     Set _.point.idx := Range(1,_.m,1);
    537     VMatrix x = SubCol(_.P::_.xy,[[1]]);
    538     VMatrix y = SubCol(_.P::_.xy,[[2]]);
    539     Set _.x.range := Range(VMatMin(x),VMatMax(x),_.maxDistance);
    540     Set _.y.range := Range(VMatMin(y),VMatMax(y),_.maxDistance);
    541     VMatrix _.x.class := MatQuery::ClassifyRowsByRank(x,_.x.range);
    542     VMatrix _.y.class := MatQuery::ClassifyRowsByRank(y,_.y.range);
    543     Set xy.class = Classify(_.point.idx, Real(Real a, Real b)
    544     {
    545       Real x.cmp = Compare(VMatDat(_.x.class,a,1),VMatDat(_.x.class,b,1));
    546       If(x.cmp, x.cmp,
    547       Real y.cmp = Compare(VMatDat(_.y.class,a,1),VMatDat(_.y.class,b,1)))
    548     });
    549     Set _.grid := For(1,Card(_.x.range)+1,Set(Real i)
    550     {
    551       For(1,Card(_.y.range)+1,Set(Real j)
    552       {
    553         Copy(Empty)
    554       })
    555     });
    556     VMatrix _.grid.map := Group("ConcatRows",EvalSet(xy.class, VMatrix(Set cl)
    557     {
    558       Real k = cl[1];
    559       Real x.class = VMatDat(_.x.class,k,1);
    560       Real y.class = VMatDat(_.y.class,k,1);
    561       Set _.grid[x.class][y.class] := cl;
    562       SetRow([[x.class,y.class]])
    563     }));
    564     True
    565   };
    566 
    567   ////////////////////////////////////////////////////////////////////////////
    568   Static Text _.autodoc.member.get.distance =
    569   "Evaluates the weight between a given point p and each point of a 2 "
    570   "columns matrix of coordinates";
    571   VMatrix get.distance(@Point2D p, VMatrix xy)
    572   ////////////////////////////////////////////////////////////////////////////
    573   {
    574     Real m = VRows(xy);
    575     VMatrix p.x = Rand(m,1,p::x(1),p::x(1));
    576     VMatrix p.y = Rand(m,1,p::y(1),p::y(1));
    577     VMatrix x = SubCol(xy,[[1]]);
    578     VMatrix y = SubCol(xy,[[2]]);
    579     RPow(RPow(x-p.x,2)+RPow(y-p.y,2),0.5)
    580   };
    581  
    582   ////////////////////////////////////////////////////////////////////////////
    583   NameBlock evaluate(@Point2D p)
    584   ////////////////////////////////////////////////////////////////////////////
    585   {[[
    586     VMatrix p.x = Rand(1,1,p::x(1),p::x(1));
    587     VMatrix p.y = Rand(1,1,p::y(1),p::y(1));
    588     VMatrix p.xy = p.x | p.y;
    589     Real p.x.class = VMatDat(MatQuery::ClassifyRowsByRank(p.x,_.x.range),1,1);
    590     Real p.y.class = VMatDat(MatQuery::ClassifyRowsByRank(p.y,_.y.range),1,1);
    591     Real i.min = Max(1,p.x.class-1);
    592     Real i.max = Min(Card(_.x.range)+1,p.x.class+1);
    593     Real j.min = Max(1,p.y.class-1);
    594     Real j.max = Min(Card(_.y.range)+1,p.y.class+1);
    595     Set p.candidates = BinGroup("<<",For(i.min, i.max, Set(Real i)
    596     {
    597       BinGroup("<<",For(j.min, j.max, Set(Real j)
    598       {
    599         _.grid[i][j] 
    600       }))
    601     }));
    602     VMatrix xy.candidates = SubRow(_.P::_.xy, p.candidates);
    603     VMatrix d.candidates = get.distance(p,xy.candidates);
    604     Set neighborhood = MatQuery::SelectRowsInInterval(d.candidates,
    605       False, 0, _.maxDistance, True);
    606     Matrix d = VMat2Mat(SubRow(d.candidates, neighborhood));
    607     VMatrix v = Mat2VMat(RPow(d,2));
    608     VMatrix w = v/VMatSum(v);
    609     VMatrix f = Tra(w) * SubRow(SubRow(_.F,p.candidates),neighborhood)
    610   ]]};
    611   ////////////////////////////////////////////////////////////////////////////
    612   VMatrix eval.vector(@Point2D p)
    613   ////////////////////////////////////////////////////////////////////////////
    614   {
    615     VMatrix p.x = Rand(1,1,p::x(1),p::x(1));
    616     VMatrix p.y = Rand(1,1,p::y(1),p::y(1));
    617     VMatrix p.xy = p.x | p.y;
    618     Real p.x.class = VMatDat(MatQuery::ClassifyRowsByRank(p.x,_.x.range),1,1);
    619     Real p.y.class = VMatDat(MatQuery::ClassifyRowsByRank(p.y,_.y.range),1,1);
    620     Real i.min = Max(1,p.x.class-1);
    621     Real i.max = Min(Card(_.x.range)+1,p.x.class+1);
    622     Real j.min = Max(1,p.y.class-1);
    623     Real j.max = Min(Card(_.y.range)+1,p.y.class+1);
    624     Set p.candidates = BinGroup("<<",For(i.min, i.max, Set(Real i)
    625     {
    626       BinGroup("<<",For(j.min, j.max, Set(Real j)
    627       {
    628         _.grid[i][j] 
    629       }))
    630     }));
    631     VMatrix xy.candidates = SubRow(_.P::_.xy, p.candidates);
    632     VMatrix d.candidates = get.distance(p,xy.candidates);
    633     Set neighborhood = MatQuery::SelectRowsInInterval(d.candidates,
    634       False, 0, _.maxDistance, True);
    635     Matrix d = VMat2Mat(SubRow(d.candidates, neighborhood),1);
    636     VMatrix v = Mat2VMat(RPow(d,2));
    637     VMatrix w = v/VMatSum(v);
    638     w * SubRow(SubRow(_.F,p.candidates),neighborhood)
    639   };
    640 
    641   Static Text _.autodoc.member.New = "Creates an instance of "
    642   "@Interp2D.Gravit";
    643 
    644   Static @Interp2D.Gravit New(VMatrix XY, VMatrix F, Real maxDistance)
    645   {
    646     @Interp2D.Gravit new =
    647     [[
    648       @Point2D.Stored _.P = @Point2D.Stored::New(XY);
    649       VMatrix _.F = F;
    650       Real _.maxDistance = maxDistance
    651     ]]
    652   }
    653 }
     21#Embed "UTMFormula.tol";
     22#Embed "KML.tol";
     23#Embed "Point2D.tol";
     24#Embed "Interp2D.tol"
    65425
    65526]];
    656 
    657 
    658 
    659 
    660 
Note: See TracChangeset for help on using the changeset viewer.