1 #include <ossim/base/ossimNadconNasDatum.h>
2
3 RTTI_DEF1(ossimNadconNasDatum, "ossimNadconNasDatum", ossimNadconGridDatum);
ossimNadconNasDatum(const ossimFilename & nadconDirectory)4 ossimNadconNasDatum::ossimNadconNasDatum(const ossimFilename& nadconDirectory)
5 :
6 ossimNadconGridDatum(nadconDirectory,
7 ossimString("NAS"),
8 ossimString("NORTH AMERICAN 1927 NADCON, CONUS GRID"),
9 ossimEllipsoidFactory::instance()->create(ossimString("CC")),
10 5.000000000000000, 5.000000000000000, 6.000000000000000,
11 -135, -60.0, 15.0, 60.0, -8.0, 160.0, 176.0)
12 {
13 }
14
shift(const ossimGpt & aPt) const15 ossimGpt ossimNadconNasDatum::shift(const ossimGpt& aPt)const
16 {
17 const ossimDatum* datum = aPt.datum();
18 ossimString code = datum->code();
19 ossimString subCode(code.begin(),
20 code.begin() + 3);
21 if(subCode == "NAS")
22 {
23 return aPt;
24 }
25 else
26 {
27 if(subCode == "NAR")
28 {
29 checkGrid(aPt);
30 if(!theLatGrid.getFileOkFlag()||
31 !theLonGrid.getFileOkFlag())
32 {
33 return ossimThreeParamDatum::shift(aPt);
34 }
35 ossimDpt pt;
36 ossimDpt tempPt;
37 ossimDpt deltaPt;
38 double shiftLat;
39 double shiftLon;
40 double minLat = theCurrentGridRect.ll().lat;
41 double maxLat = theCurrentGridRect.ul().lat;
42 double minLon = theCurrentGridRect.ul().lon;
43 double maxLon = theCurrentGridRect.ur().lon;
44 int maxIter = 20;
45 double epsilon = 1.0e-9;
46 int c = 0;
47
48 pt.lat = aPt.latd();
49 pt.lon = aPt.lond();
50 ossimDpt pt2;
51 tempPt = pt;
52 ossimDpt delta2;
53 while(c < maxIter)
54 {
55 if(tempPt.lat < minLat) tempPt.lat = minLat;
56 if(tempPt.lat > maxLat) tempPt.lat = maxLat;
57 if(tempPt.lon < minLon) tempPt.lon = minLon;
58 if(tempPt.lon > maxLon) tempPt.lon = maxLon;
59
60 shiftLat = theLatGrid.getShiftAtLatLon(tempPt.lat, tempPt.lon);
61 shiftLon = theLonGrid.getShiftAtLatLon(tempPt.lat, tempPt.lon);
62
63
64 pt2.lat = tempPt.lat + shiftLat/3600.0;
65 pt2.lon = tempPt.lon - shiftLon/3600.0;
66
67 delta2 = pt2-pt;
68 if(c == 0)
69 {
70 if(fabs(deltaPt.lon) > epsilon)
71 {
72 tempPt.lon = pt.lon - shiftLon/3600.0;
73 }
74 if(fabs(deltaPt.lat) > epsilon)
75 {
76 tempPt.lat = pt.lat + shiftLat/3600.0;
77 }
78 }
79 else
80 {
81 if(fabs(delta2.lon) > epsilon)
82 {
83 tempPt.lon = tempPt.lon - delta2.lon;
84 }
85 if(fabs(delta2.lat) > epsilon)
86 {
87 tempPt.lat =tempPt.lat - delta2.lat;
88 }
89 }
90 if((fabs(delta2.lat) < epsilon)&&
91 (fabs(delta2.lon) < epsilon))
92 {
93 break;
94 }
95 else
96 {
97 ++c;
98 }
99 }
100 return ossimGpt(tempPt.lat,
101 tempPt.lon,
102 aPt.height(),
103 this);
104 }
105 }
106 return ossimThreeParamDatum::shift(aPt);
107 }
108