Orfeo Toolbox  4.0
otbMapProjectionAdapter.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 
20 
21 #include <cassert>
22 
23 #include "otbMacro.h"
24 #include "otbUtils.h"
25 
26 #include "ossim/projection/ossimMapProjection.h"
27 #include "ossim/projection/ossimMapProjectionFactory.h"
28 #include "ossim/base/ossimGpt.h"
29 #include "ossim/base/ossimDpt.h"
30 #include "ossim/projection/ossimProjection.h"
31 #include "ossim/base/ossimEllipsoid.h"
32 #include "ossim/base/ossimEllipsoidFactory.h"
33 #include "ossim/base/ossimString.h"
34 #include "gdal/ossimOgcWktTranslator.h"
35 
36 #include "ossim/projection/ossimUtmProjection.h"
37 #include "ossim/projection/ossimLambertConformalConicProjection.h"
38 #include "ossim/projection/ossimTransMercatorProjection.h"
39 #include "ossim/projection/ossimEckert4Projection.h"
40 #include "ossim/projection/ossimMollweidProjection.h"
41 #include "ossim/projection/ossimSinusoidalProjection.h"
42 
43 #include "ossim/support_data/ossimSpaceImagingGeom.h"
44 #include "ossim/base/ossimKeywordNames.h"
45 
46 
47 namespace otb
48 {
49 
51  m_MapProjection(NULL), m_ProjectionRefWkt(""), m_ReinstanciateProjection(true)
52 {
53 }
54 
56 {
57  if (m_MapProjection != NULL)
58  {
59  delete m_MapProjection;
60  }
61 }
62 
64 {
65  itkDebugMacro("returning MapProjection address " << this->m_MapProjection);
67  {
68  this->InstanciateProjection();
69  }
70 
71  return this->m_MapProjection;
72 }
73 
75 {
76  itkDebugMacro("returning MapProjection address " << this->m_MapProjection);
78  {
79  itkExceptionMacro(<< "m_MapProjection not up-to-date, call InstanciateProjection() first");
80  }
81 
82  return this->m_MapProjection;
83 }
84 
85 std::string MapProjectionAdapter::GetWkt() const
86 {
87  ossimKeywordlist kwl;
88  this->GetMapProjection();
89  if (m_MapProjection == NULL)
90  {
91  return "";
92  }
93  m_MapProjection->saveState(kwl);
94  ossimOgcWktTranslator wktTranslator;
95  std::string wkt;
96  wkt = wktTranslator.fromOssimKwl(kwl).chars();
97  return wkt;
98 }
99 
100 void MapProjectionAdapter::SetWkt(const std::string& projectionRefWkt)
101 {
102  this->m_ProjectionRefWkt = projectionRefWkt;
104  this->InstanciateProjection(); //Should not be needed, but it is...
105  this->Modified();
106 }
107 
108 void MapProjectionAdapter::SetParameter(const std::string& key, const std::string& value)
109 {
110  m_ParameterStore[key] = value;
112  this->InstanciateProjection(); //Should not be needed, but it is...
113  this->Modified();
114 }
115 
116 std::string MapProjectionAdapter::GetParameter(const std::string& key) const
117 {
118  // Please refer to the note in the header filer
119  // we do NOT want to read from m_ParameterStore here!
120 
121  std::string projectionName = this->GetMapProjection()->getClassName().string();
122 
123  // Start by matching generic parameters
124  const ossimMapProjection* projection = dynamic_cast<const ossimMapProjection*>(this->GetMapProjection());
125  if (key.compare("Origin") == 0)
126  {
127  return Utils::ConvertToString(projection->origin());
128  }
129  if (key.compare("FalseNorthing") == 0)
130  {
131  return Utils::ConvertToString(projection->getFalseNorthing());
132  }
133  if (key.compare("FalseEasting") == 0)
134  {
135  return Utils::ConvertToString(projection->getFalseEasting());
136  }
137  if (key.compare("StandardParallel1") == 0)
138  {
139  return Utils::ConvertToString(projection->getStandardParallel1());
140  }
141  if (key.compare("StandardParallel2") == 0)
142  {
143  return Utils::ConvertToString(projection->getStandardParallel2());
144  }
145  if (key.compare("A") == 0)
146  {
147  return Utils::ConvertToString(projection->getA());
148  }
149  if (key.compare("B") == 0)
150  {
151  return Utils::ConvertToString(projection->getB());
152  }
153  if (key.compare("F") == 0)
154  {
155  return Utils::ConvertToString(projection->getF());
156  }
157  if (key.compare("MetersPerPixel") == 0)
158  {
159  return Utils::ConvertToString(projection->getMetersPerPixel());
160  }
161  if (key.compare("DecimalDegreesPerPixel") == 0)
162  {
163  return Utils::ConvertToString(projection->getDecimalDegreesPerPixel());
164  }
165 
166  // Apply parameters to transmercator
167  if (projectionName.compare("ossimTransMercatorProjection") == 0)
168  {
169  const ossimTransMercatorProjection* projection = dynamic_cast<const ossimTransMercatorProjection*>(this->GetMapProjection());
170  if (key.compare("ScaleFactor") == 0)
171  {
172  return Utils::ConvertToString(projection->getScaleFactor());
173  }
174  }
175 
176  // Apply parameters to Utm
177  if (projectionName.compare("ossimUtmProjection") == 0)
178  {
179  const ossimUtmProjection* projection = dynamic_cast<const ossimUtmProjection*>(this->GetMapProjection());
180  if (key.compare("Zone") == 0)
181  {
182  return Utils::ConvertToString(projection->getZone());
183  }
184  if (key.compare("Hemisphere") == 0)
185  {
186  return Utils::ConvertToString(projection->getHemisphere());
187  }
188  }
189 
190  return "";
191 }
192 
194 {
195  if ((this->m_ReinstanciateProjection) || (m_MapProjection == NULL))
196  {
197  ossimKeywordlist kwl;
198  ossimOgcWktTranslator wktTranslator;
199 
200  bool projectionInformationAvailable = wktTranslator.toOssimKwl(m_ProjectionRefWkt, kwl);
201 
202  if (projectionInformationAvailable)
203  {
204  //we don't want to have a ossimEquDistCylProjection here:
205  //see discussion in May 2009 on ossim list;
206  //a better solution might be available...
207  std::string projectionString(kwl.find("type"));
208  if (projectionString.find("ossimEquDistCylProjection") != string::npos)
209  {
210  otbMsgDevMacro(<< "WARNING: Not instanciating a ossimEquDistCylProjection: " << projectionString);
211  otbMsgDevMacro(<< "Wkt was: " << kwl);
212  otbMsgDevMacro(<< "From RefWkt: " << m_ProjectionRefWkt);
213  return false;
214  }
215 
216  m_MapProjection = ossimMapProjectionFactory::instance()->createProjection(kwl);
217  }
218  else
219  {
220  otbMsgDevMacro(<< "WARNING: Impossible to create the projection from Wkt: " << m_ProjectionRefWkt);
221  otbMsgDevMacro(<< "Trying with string as a string (ossimUtmProjection or Utm would qualify");
222  // Trying directly with the m_ProjectionRefWkt (is
223  // ossimUtmProjection for example)
224  ossimString name(m_ProjectionRefWkt);
225  m_MapProjection = ossimMapProjectionFactory::instance()->createProjection(name);
226  if (m_MapProjection == NULL)
227  {
228  // Trying directly extending the m_ProjectionRefWkt (convert the
229  // Utm to ossimUtmProjection for example)
230  ossimString extendedName("ossim");
231  extendedName += m_ProjectionRefWkt;
232  extendedName += "Projection";
233  m_MapProjection = ossimMapProjectionFactory::instance()->createProjection(extendedName);
234  }
235 
236  if (m_MapProjection == NULL) return false;
237 
238  }
239 
240  this->m_ReinstanciateProjection = false;
242  return true;
243  }
244  return false;
245 }
246 
247 void MapProjectionAdapter::InverseTransform(double x, double y, double z,
248  double& lon, double& lat, double& h)
249 {
250  if (m_MapProjection == NULL)
251  {
252  otbMsgDevMacro(<< "WARNING: Using identity");
253  lon = x;
254  lat = y;
255  h = z;
256  return;
257  }
258 
259  ossimDpt ossimDPoint(x, y);
260 
261  //map projection
262  ossimGpt ossimGPoint;
263  ossimGPoint = this->GetMapProjection()->inverse(ossimDPoint);
264  ossimGPoint.changeDatum(ossimDatumFactory::instance()->wgs84());
265 
266  lon = ossimGPoint.lon;
267  lat = ossimGPoint.lat;
268  h = z;
269 }
270 
271 void MapProjectionAdapter::ForwardTransform(double lon, double lat, double h,
272  double& x, double& y, double& z)
273 {
274  if (m_MapProjection == NULL)
275  {
276  otbMsgDevMacro(<< "WARNING: Using identity");
277  x = lon;
278  y = lat;
279  z = h;
280  return;
281  }
282 
283  ossimGpt ossimGPoint(lat, lon, h);
284 
285  //map projection
286  ossimDpt ossimDPoint;
287  ossimDPoint = this->GetMapProjection()->forward(ossimGPoint);
288 
289  x = ossimDPoint.x;
290  y = ossimDPoint.y;
291  z = h;
292 }
293 
295 {
296  // Start by identifying the projection, that will be necessary for
297  // the casting.
298  std::string projectionName = this->GetMapProjection()->getClassName().string();
299 
300  StoreType::const_iterator it;
301 
302  // Apply standard map projection parameters
303  ossimMapProjection* projection = dynamic_cast<ossimMapProjection*>(this->GetMapProjection());
304  // Set up origin
305 
306  const ossimDatum* datum = ossimDatumFactory::instance()->wgs84(); //default value
307  it = m_ParameterStore.find("Datum");
308  if (it != m_ParameterStore.end())
309  {
310  datum = ossimDatumFactory::instance()->create((*it).second);
311  projection->setDatum(datum);
312  }
313 
314  StoreType::const_iterator itX = m_ParameterStore.find("OriginX");
315  StoreType::const_iterator itY = m_ParameterStore.find("OriginY");
316  StoreType::const_iterator itZ = m_ParameterStore.find("OriginZ");
317 
318  if (itX != m_ParameterStore.end() && itY != m_ParameterStore.end())
319  {
320  double originX = atof((*itX).second.c_str());
321  double originY = atof((*itY).second.c_str());
322  double originZ = 0;
323  if (itZ != m_ParameterStore.end())
324  {
325  originZ = atof((*itZ).second.c_str());
326  }
327  ossimGpt origin(originY, originX, originZ, datum);
328  projection->setOrigin(origin);
329  }
330 
331  // Set up resolution
332  StoreType::const_iterator itResMeterX = m_ParameterStore.find("MetersPerPixelX");
333  StoreType::const_iterator itResMeterY = m_ParameterStore.find("MetersPerPixelY");
334  if (itResMeterX != m_ParameterStore.end() && itResMeterY != m_ParameterStore.end())
335  {
336  double resMeterX = atof((*itResMeterX).second.c_str());
337  double resMeterY = atof((*itResMeterY).second.c_str());
338  ossimDpt resMeter(resMeterX, resMeterY);
339  projection->setMetersPerPixel(resMeter);
340  }
341 
342  // Apply parameters to LambertConformalConic
343  if (projectionName.compare("ossimLambertConformalConicProjection") == 0)
344  {
345  ossimLambertConformalConicProjection* projection = dynamic_cast<ossimLambertConformalConicProjection*>(this->GetMapProjection());
346 
347  it = m_ParameterStore.find("FalseNorthing");
348  if (it != m_ParameterStore.end())
349  {
350  double value = atof((*it).second.c_str());
351 
352  projection->setFalseNorthing(value);
353  }
354  it = m_ParameterStore.find("FalseEasting");
355  if (it != m_ParameterStore.end())
356  {
357  double value = atof((*it).second.c_str());
358 
359  projection->setFalseEasting(value);
360  }
361  it = m_ParameterStore.find("StandardParallel1");
362  if (it != m_ParameterStore.end())
363  {
364  double value = atof((*it).second.c_str());
365  projection->setStandardParallel1(value);
366  }
367  it = m_ParameterStore.find("StandardParallel2");
368  if (it != m_ParameterStore.end())
369  {
370  double value = atof((*it).second.c_str());
371  projection->setStandardParallel2(value);
372  }
373  }
374 
375  // Apply parameters to Eckert4
376  if (projectionName.compare("ossimEckert4Projection") == 0)
377  {
378  ossimEckert4Projection* projection = dynamic_cast<ossimEckert4Projection*>(this->GetMapProjection());
379 
380  it = m_ParameterStore.find("FalseNorthing");
381  if (it != m_ParameterStore.end())
382  {
383  double value = atof((*it).second.c_str());
384  projection->setFalseNorthing(value);
385  }
386  it = m_ParameterStore.find("FalseEasting");
387  if (it != m_ParameterStore.end())
388  {
389  double value = atof((*it).second.c_str());
390  projection->setFalseEasting(value);
391  }
392  }
393 
394  // Apply parameters to Mollweid
395  if (projectionName.compare("ossimMollweidProjection") == 0)
396  {
397  ossimMollweidProjection* projection = dynamic_cast<ossimMollweidProjection*>(this->GetMapProjection());
398 
399  it = m_ParameterStore.find("FalseNorthing");
400  if (it != m_ParameterStore.end())
401  {
402  double value = atof((*it).second.c_str());
403  projection->setFalseNorthing(value);
404  }
405  it = m_ParameterStore.find("FalseEasting");
406  if (it != m_ParameterStore.end())
407  {
408  double value = atof((*it).second.c_str());
409  projection->setFalseEasting(value);
410  }
411  }
412 
413  // Apply parameters to Sinusoidal
414  if (projectionName.compare("ossimSinusoidalProjection") == 0)
415  {
416  ossimSinusoidalProjection* projection = dynamic_cast<ossimSinusoidalProjection*>(this->GetMapProjection());
417 
418  it = m_ParameterStore.find("FalseNorthing");
419  if (it != m_ParameterStore.end())
420  {
421  double value = atof((*it).second.c_str());
422  projection->setFalseNorthing(value);
423  }
424  it = m_ParameterStore.find("FalseEasting");
425  if (it != m_ParameterStore.end())
426  {
427  double value = atof((*it).second.c_str());
428  projection->setFalseEasting(value);
429  }
430  }
431 
432  // Apply parameters to transmercator
433  if (projectionName.compare("ossimTransMercatorProjection") == 0)
434  {
435  ossimTransMercatorProjection* projection = dynamic_cast<ossimTransMercatorProjection*> (this->GetMapProjection());
436  it = m_ParameterStore.find("ScaleFactor");
437  if (it != m_ParameterStore.end())
438  {
439  double scale = atof((*it).second.c_str());
440  projection->setScaleFactor(scale);
441  }
442  it = m_ParameterStore.find("FalseNorthing");
443  if (it != m_ParameterStore.end())
444  {
445  double value = atof((*it).second.c_str());
446  projection->setFalseNorthing(value);
447  }
448  it = m_ParameterStore.find("FalseEasting");
449  if (it != m_ParameterStore.end())
450  {
451  double value = atof((*it).second.c_str());
452  projection->setFalseEasting(value);
453  }
454  }
455 
456  // Apply parameters to Utm
457  if (projectionName.compare("ossimUtmProjection") == 0)
458  {
459  ossimUtmProjection* projection = dynamic_cast<ossimUtmProjection*>(this->GetMapProjection());
460  it = m_ParameterStore.find("Zone");
461  if (it != m_ParameterStore.end())
462  {
463  int zone = atoi((*it).second.c_str());
464  projection->setZone(zone);
465  }
466  it = m_ParameterStore.find("Hemisphere");
467  if (it != m_ParameterStore.end())
468  {
469  projection->setHemisphere((*it).second[0]);
470  }
471  }
472 }
473 
475 {
476  m_MapProjection->print(std::cout);
477  std::cout << "Parameter store:\n";
478  for (StoreType::const_iterator it = m_ParameterStore.begin();
479  it != m_ParameterStore.end();
480  ++it)
481  {
482  std::cout << " " << (*it).first << ": " << (*it).second << "\n";
483  }
484 }
485 
486 
487 namespace Utils
488 {
489 
490 int GetZoneFromGeoPoint(double lon, double lat)
491 {
492  //use ossim to handle the special case of UTM
493  ossimGpt point(lat, lon);
494  ossimUtmProjection projection;
495  int zone = projection.computeZone(point);
496  return zone;
497 }
498 
499 }
500 
501 } // namespace otb

Generated at Sat Mar 8 2014 16:06:41 for Orfeo Toolbox with doxygen 1.8.3.1