1 //*******************************************************************
2 // Copyright (C) 2001 ImageLinks Inc.
3 //
4 // License:  See top level LICENSE.txt file.
5 //
6 // Author:  David Burken
7 //
8 // Description:
9 // Class for encapsulate projection info given a projection, datum, and
10 // output rectangle.
11 //
12 // NOTE:
13 // - Output rectangle should be relative to the center of pixels.
14 //   Shifts will be made for "pixel is area" internally.
15 //
16 //*******************************************************************
17 //  $Id: ossimMapProjectionInfo.cpp 13025 2008-06-13 17:06:30Z sbortman $
18 
19 #include <cstdlib>
20 #include <fstream>
21 #include <time.h>
22 
23 using namespace std;
24 
25 #include <ossim/projection/ossimMapProjectionInfo.h>
26 #include <ossim/projection/ossimMapProjection.h>
27 #include <ossim/base/ossimTrace.h>
28 #include <ossim/base/ossimCommon.h>
29 #include <ossim/base/ossimFilename.h>
30 #include <ossim/base/ossimDatum.h>
31 #include <ossim/base/ossimEllipsoid.h>
32 #include <ossim/base/ossimKeywordlist.h>
33 #include <ossim/base/ossimDrect.h>
34 #include <ossim/base/ossimKeywordNames.h>
35 #include <ossim/base/ossimNotifyContext.h>
36 #include <ossim/base/ossimUnitTypeLut.h>
37 
38 //***
39 // Static trace for debugging.
40 //***
41 static ossimTrace traceDebug("ossimMapProjectionInfo::debug");
42 
43 //***
44 // Keywords for getStateFrom/saveStateTo.
45 //***
46 const char* ossimMapProjectionInfo::OUTPUT_US_FT_INFO_KW =
47 "viewinfo.output_readme_in_us_ft_flag";
48 
49 const char* ossimMapProjectionInfo::PIXEL_TYPE_KW = "viewinfo.pixel_type";
50 
51 const char* ossimMapProjectionInfo::README_IMAGE_STRING_KW =
52 "viewinfo.readme_image_string";
53 
ossimMapProjectionInfo(const ossimMapProjection * proj,const ossimDrect & output_rect)54 ossimMapProjectionInfo::ossimMapProjectionInfo(const ossimMapProjection* proj,
55                                                const ossimDrect& output_rect)
56    :
57       theProjection               (proj),
58       theErrorStatus              (false),
59       theLinesPerImage            (0),
60       thePixelsPerLine            (0),
61       theCornerGroundPt           (),
62       theCornerEastingNorthingPt  (),
63       theCenterGroundPt           (),
64       theCenterEastingNorthingPt  (0.0, 0.0),
65       thePixelType                (OSSIM_PIXEL_IS_POINT),
66       theOutputInfoInFeetFlag     (false),
67       theImageInfoString          ()
68 {
69    if (!theProjection)
70    {
71       theErrorStatus = true;
72       ossimNotify(ossimNotifyLevel_FATAL)
73          << "FATAL ossimMapProjectionInfo::ossimMapProjectionInfo: "
74          << "Null projection pointer passed to constructor!"
75          << "\nError status has been set.  Returning..."
76          << std::endl;
77       return;
78    }
79 
80    initializeMembers(output_rect);
81 
82    if (traceDebug())
83    {
84       ossimNotify(ossimNotifyLevel_DEBUG)
85          << "DEBUG ossimMapProjectionInfo::ossimMapProjectionInfo:\n"
86          << "output_rect:  " << output_rect << "\n"
87          << *this << std::endl;
88    }
89 }
90 
~ossimMapProjectionInfo()91 ossimMapProjectionInfo::~ossimMapProjectionInfo()
92 {
93 }
94 
errorStatus() const95 bool ossimMapProjectionInfo::errorStatus() const
96 {
97    return theErrorStatus;
98 }
99 
initializeMembers(const ossimDrect & rect)100 void ossimMapProjectionInfo::initializeMembers(const ossimDrect& rect)
101 {
102    theBoundingRect = rect;
103    theLinesPerImage  = ossim::round<int>(rect.height());
104    thePixelsPerLine  = ossim::round<int>(rect.width());
105 
106    theProjection->lineSampleToWorld(rect.ul(),
107                                     theCornerGroundPt[0]);
108 
109    theProjection->lineSampleToWorld(rect.ur(),
110                                     theCornerGroundPt[1]);
111 
112    theProjection->lineSampleToWorld(rect.lr(),
113                                     theCornerGroundPt[2]);
114 
115 
116    theProjection->lineSampleToWorld(rect.ll(),
117                                     theCornerGroundPt[3]);
118 
119    theProjection->lineSampleToEastingNorthing(rect.ul(),
120                                               theCornerEastingNorthingPt[0]);
121 
122    theProjection->lineSampleToEastingNorthing(rect.ur(),
123                                               theCornerEastingNorthingPt[1]);
124 
125    theProjection->lineSampleToEastingNorthing(rect.lr(),
126                                               theCornerEastingNorthingPt[2]);
127 
128    theProjection->lineSampleToEastingNorthing(rect.ll(),
129                                               theCornerEastingNorthingPt[3]);
130 
131    theCenterEastingNorthingPt.x = (theCornerEastingNorthingPt[0].x +
132                                    theCornerEastingNorthingPt[1].x +
133                                    theCornerEastingNorthingPt[2].x +
134                                    theCornerEastingNorthingPt[3].x) / 4.0;
135 
136    theCenterEastingNorthingPt.y = (theCornerEastingNorthingPt[0].y +
137                                    theCornerEastingNorthingPt[1].y +
138                                    theCornerEastingNorthingPt[2].y +
139                                    theCornerEastingNorthingPt[3].y) / 4.0;
140 
141    theCenterGroundPt = theProjection->inverse(theCenterEastingNorthingPt);
142 }
143 
144 
getGeom(ossimKeywordlist & kwl,const char * prefix) const145 void ossimMapProjectionInfo::getGeom(ossimKeywordlist& kwl,
146                                      const char* prefix)const
147 {
148    theProjection->saveState(kwl, prefix);
149 
150    if(theProjection->isGeographic())
151    {
152       ossimGpt gpt = ulGroundPt();
153       kwl.add(prefix,
154               ossimKeywordNames::TIE_POINT_XY_KW,
155               ossimDpt(gpt).toString().c_str(),
156               true);
157       kwl.add(prefix,
158               ossimKeywordNames::TIE_POINT_UNITS_KW,
159               ossimUnitTypeLut::instance()->getEntryString(OSSIM_DEGREES),
160               true);
161    }
162    else
163    {
164       ossimDpt dpt = ulEastingNorthingPt();
165       kwl.add(prefix,
166               ossimKeywordNames::TIE_POINT_XY_KW,
167               dpt.toString().c_str(),
168               true);
169       kwl.add(prefix,
170               ossimKeywordNames::TIE_POINT_UNITS_KW,
171               ossimUnitTypeLut::instance()->getEntryString(OSSIM_METERS),
172               true);
173    }
174 }
175 
print(std::ostream & os) const176 std::ostream& ossimMapProjectionInfo::print(std::ostream& os) const
177 {
178    if (!os)
179    {
180       return os;
181    }
182 
183    os << setiosflags(ios::left)
184       << setiosflags(ios::fixed)
185       << "ossimMapProjectionInfo Data Members:\n"
186       << "Projection name:  " << theProjection->getProjectionName()
187       << setw(30) << "\nOutput pixel type:"
188       << ((getPixelType() == OSSIM_PIXEL_IS_POINT) ? "pixel is point" :
189           "pixel is area")
190       << setw(30) << "\nMeters per pixel:"
191       << getMetersPerPixel()
192       << setw(30) << "\nUS survey feet per pixel:"
193       << getUsSurveyFeetPerPixel()
194       << setw(30) << "\nDecimal degrees per pixel:"
195       << getDecimalDegreesPerPixel()
196       << setw(30) << "\nNumber of lines:"
197       << theLinesPerImage
198       << setw(30) << "\nNumber of pixels:"
199       << thePixelsPerLine
200       << setw(30) << "\nUpper left ground point:"
201       << theCornerGroundPt[0]
202       << setw(30) << "\nUpper right ground point:"
203       << theCornerGroundPt[1]
204       << setw(30) << "\nLower right ground point:"
205       << theCornerGroundPt[2]
206       << setw(30) << "\nLower left ground point:"
207       << theCornerGroundPt[3]
208       << setw(30) << "\nUpper left easting_northing:"
209       << theCornerEastingNorthingPt[0]
210       << setw(30) << "\nUpper right easting_northing:"
211       << theCornerEastingNorthingPt[1]
212       << setw(30) << "\nLower right easting_northing:"
213       << theCornerEastingNorthingPt[2]
214       << setw(30) << "\nLower left easting_northing:"
215       << theCornerEastingNorthingPt[3]
216       << setw(30) << "\nCenter ground point:"
217       << theCenterGroundPt
218       << setw(30) << "\nCenter easting_northing:"
219       << theCenterEastingNorthingPt
220       << "\nMap Projection dump:";
221    theProjection->print(os);
222 
223    return os;
224 }
225 
setPixelType(ossimPixelType type)226 void ossimMapProjectionInfo::setPixelType (ossimPixelType type)
227 {
228    thePixelType = type;
229 }
230 
getPixelType() const231 ossimPixelType ossimMapProjectionInfo::getPixelType () const
232 {
233    return thePixelType;
234 }
235 
setOutputFeetFlag(bool flag)236 void ossimMapProjectionInfo::setOutputFeetFlag (bool flag)
237 {
238    theOutputInfoInFeetFlag = flag;
239 }
240 
unitsInFeet() const241 bool ossimMapProjectionInfo::unitsInFeet() const
242 {
243    return theOutputInfoInFeetFlag;
244 }
245 
getImageInfoString() const246 ossimString ossimMapProjectionInfo::getImageInfoString () const
247 {
248    return theImageInfoString;
249 }
250 
setImageInfoString(const ossimString & string)251 void ossimMapProjectionInfo::setImageInfoString (const ossimString& string)
252 {
253    theImageInfoString = string;
254 }
255 
getProjection() const256 const ossimMapProjection* ossimMapProjectionInfo::getProjection() const
257 {
258    return theProjection;
259 }
260 
ulEastingNorthingPt() const261 ossimDpt ossimMapProjectionInfo::ulEastingNorthingPt( ) const
262 {
263    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
264    {
265       ossimDpt mpp = getMetersPerPixel();
266       ossimDpt pt;
267       pt.x = theCornerEastingNorthingPt[0].x - (mpp.x / 2.0);
268       pt.y = theCornerEastingNorthingPt[0].y + (mpp.y / 2.0);
269       return pt;
270    }
271    else
272    {
273       return theCornerEastingNorthingPt[0];
274    }
275 }
276 
urEastingNorthingPt() const277 ossimDpt ossimMapProjectionInfo::urEastingNorthingPt( ) const
278 {
279    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
280    {
281       ossimDpt mpp = getMetersPerPixel();
282       ossimDpt pt;
283       pt.x = theCornerEastingNorthingPt[1].x + (mpp.x / 2.0);
284       pt.y = theCornerEastingNorthingPt[1].y + (mpp.y / 2.0);
285       return pt;
286    }
287    else
288    {
289       return theCornerEastingNorthingPt[1];
290    }
291 }
292 
lrEastingNorthingPt() const293 ossimDpt ossimMapProjectionInfo::lrEastingNorthingPt( ) const
294 {
295    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
296    {
297       ossimDpt mpp = getMetersPerPixel();
298       ossimDpt pt;
299       pt.x = theCornerEastingNorthingPt[2].x + (mpp.x / 2.0);
300       pt.y = theCornerEastingNorthingPt[2].y - (mpp.y / 2.0);
301       return pt;
302    }
303    else
304    {
305       return theCornerEastingNorthingPt[2];
306    }
307 }
308 
llEastingNorthingPt() const309 ossimDpt ossimMapProjectionInfo::llEastingNorthingPt( ) const
310 {
311    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
312    {
313       ossimDpt mpp = getMetersPerPixel();
314       ossimDpt pt;
315       pt.x = theCornerEastingNorthingPt[3].x - (mpp.x / 2.0);
316       pt.y = theCornerEastingNorthingPt[3].y - (mpp.y / 2.0);
317       return pt;
318    }
319    else
320    {
321       return theCornerEastingNorthingPt[3];
322    }
323 }
324 
ulGroundPt() const325 ossimGpt ossimMapProjectionInfo::ulGroundPt( ) const
326 {
327    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
328    {
329       ossimDpt ddpp = getDecimalDegreesPerPixel();
330       ossimGpt gpt;
331       gpt.latd(theCornerGroundPt[0].latd() + (ddpp.y / 2.0));
332       gpt.lond(theCornerGroundPt[0].lond() - (ddpp.x / 2.0));
333       return gpt;
334    }
335    else
336    {
337       return theCornerGroundPt[0];
338    }
339 }
340 
urGroundPt() const341 ossimGpt ossimMapProjectionInfo::urGroundPt( ) const
342 {
343    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
344    {
345       ossimDpt ddpp = getDecimalDegreesPerPixel();
346       ossimGpt gpt;
347       gpt.latd(theCornerGroundPt[1].latd() + (ddpp.y / 2.0));
348       gpt.lond(theCornerGroundPt[1].lond() + (ddpp.x / 2.0));
349       return gpt;
350    }
351    else
352    {
353       return theCornerGroundPt[1];
354    }
355 }
356 
lrGroundPt() const357 ossimGpt ossimMapProjectionInfo::lrGroundPt( ) const
358 {
359    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
360    {
361       ossimDpt ddpp = getDecimalDegreesPerPixel();
362       ossimGpt gpt;
363       gpt.latd(theCornerGroundPt[2].latd() - (ddpp.y / 2.0));
364       gpt.lond(theCornerGroundPt[2].lond() + (ddpp.x / 2.0));
365       return gpt;
366    }
367    else
368    {
369       return theCornerGroundPt[2];
370    }
371 }
372 
llGroundPt() const373 ossimGpt ossimMapProjectionInfo::llGroundPt( ) const
374 {
375    if (getPixelType() == OSSIM_PIXEL_IS_AREA)
376    {
377       ossimDpt ddpp = getDecimalDegreesPerPixel();
378       ossimGpt gpt;
379       gpt.latd(theCornerGroundPt[3].latd() - (ddpp.y / 2.0));
380       gpt.lond(theCornerGroundPt[3].lond() - (ddpp.x / 2.0));
381       return gpt;
382    }
383    else
384    {
385       return theCornerGroundPt[3];
386    }
387 }
388 
centerGroundPt() const389 ossimGpt ossimMapProjectionInfo::centerGroundPt( ) const
390 {
391    // Center is simply center, no shift for pixel is area...
392    return theCenterGroundPt;
393 }
394 
centerEastingNorthingPt() const395 ossimDpt ossimMapProjectionInfo::centerEastingNorthingPt() const
396 {
397    return theCenterEastingNorthingPt;
398 }
399 
linesPerImage() const400 ossim_int32 ossimMapProjectionInfo::linesPerImage() const
401 {
402    return theLinesPerImage;
403 }
404 
pixelsPerLine() const405 ossim_int32 ossimMapProjectionInfo::pixelsPerLine() const
406 {
407    return thePixelsPerLine;
408 }
409 
loadState(const ossimKeywordlist & kwl,const char *)410 bool ossimMapProjectionInfo::loadState(const ossimKeywordlist& kwl,
411                                        const char* )
412 {
413    const char* lookupReturn = kwl.find(README_IMAGE_STRING_KW);
414 
415    //***
416    // "theImageInfoString" goes in the README file "Image:" field and can
417    // be used to identify the image.
418    //***
419    if (lookupReturn) theImageInfoString = lookupReturn;
420 
421    lookupReturn = kwl.find(ossimKeywordNames::PIXEL_TYPE_KW);
422 
423    if (lookupReturn)
424    {
425       ossimString tmp = lookupReturn;
426 
427       tmp.downcase();
428 
429       if (tmp.contains("area"))
430       {
431          thePixelType = OSSIM_PIXEL_IS_AREA;
432       }
433       else
434       {
435          thePixelType = OSSIM_PIXEL_IS_POINT;
436       }
437    }
438 
439    lookupReturn = kwl.find(OUTPUT_US_FT_INFO_KW);
440 
441    if (lookupReturn)
442    {
443       int tmp = atoi(lookupReturn);
444 
445       if (tmp)
446       {
447          theOutputInfoInFeetFlag = true;
448       }
449       else
450       {
451          theOutputInfoInFeetFlag = false;
452       }
453    }
454 
455    if (traceDebug())
456    {
457       ossimNotify(ossimNotifyLevel_DEBUG)
458          << "DEBUG ossimMapProjectionInfo::loadState:"
459          << "\ntheImageInfoString:       " << theImageInfoString
460          << "\nthePixelType:             " << int(thePixelType)
461          << "\ntheOutputInfoInFeetFlag:  " << theOutputInfoInFeetFlag
462          << endl;
463    }
464 
465    return true;
466 }
467 
saveState(ossimKeywordlist & kwl,const char *) const468 bool ossimMapProjectionInfo::saveState(ossimKeywordlist& kwl,
469                                        const char* ) const
470 {
471    kwl.add(README_IMAGE_STRING_KW,
472            theImageInfoString.chars());
473 
474    ossimString tmp;
475 
476    if (thePixelType == OSSIM_PIXEL_IS_POINT)
477    {
478       tmp = "point";
479    }
480    else
481    {
482       tmp = "area";
483    }
484 
485    kwl.add(PIXEL_TYPE_KW,
486            tmp.chars());
487 
488    kwl.add(OUTPUT_US_FT_INFO_KW,
489            int(theOutputInfoInFeetFlag));
490 
491    return true;
492 }
493 
ulEastingNorthingPtInFt() const494 ossimDpt ossimMapProjectionInfo::ulEastingNorthingPtInFt() const
495 {
496    ossimDpt pt = ulEastingNorthingPt();
497 
498    pt.x = ossim::mtrs2usft(pt.x);
499 
500    pt.y = ossim::mtrs2usft(pt.y);
501 
502    return pt;
503 }
504 
urEastingNorthingPtInFt() const505 ossimDpt ossimMapProjectionInfo::urEastingNorthingPtInFt() const
506 {
507    ossimDpt pt = urEastingNorthingPt();
508    pt.x = ossim::mtrs2usft(pt.x);
509    pt.y = ossim::mtrs2usft(pt.y);
510    return pt;
511 }
512 
lrEastingNorthingPtInFt() const513 ossimDpt ossimMapProjectionInfo::lrEastingNorthingPtInFt() const
514 {
515    ossimDpt pt = lrEastingNorthingPt();
516    pt.x = ossim::mtrs2usft(pt.x);
517    pt.y = ossim::mtrs2usft(pt.y);
518    return pt;
519 }
520 
llEastingNorthingPtInFt() const521 ossimDpt ossimMapProjectionInfo::llEastingNorthingPtInFt() const
522 {
523    ossimDpt pt = llEastingNorthingPt();
524    pt.x = ossim::mtrs2usft(pt.x);
525    pt.y = ossim::mtrs2usft(pt.y);
526    return pt;
527 }
528 
getMetersPerPixel() const529 ossimDpt ossimMapProjectionInfo::getMetersPerPixel() const
530 {
531    return theProjection->getMetersPerPixel();
532 }
533 
getUsSurveyFeetPerPixel() const534 ossimDpt ossimMapProjectionInfo::getUsSurveyFeetPerPixel() const
535 {
536    ossimDpt pt = getMetersPerPixel();
537 
538    pt.x = ossim::mtrs2usft(pt.x);
539    pt.y = ossim::mtrs2usft(pt.y);
540 
541    return pt;
542 }
543 
getDecimalDegreesPerPixel() const544 ossimDpt ossimMapProjectionInfo::getDecimalDegreesPerPixel() const
545 {
546    return theProjection->getDecimalDegreesPerPixel();
547 }
548