GeometriesChangeSpatialReference.cxx 7.46 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

/*===========================================================================*/
/*===============================[ Includes ]================================*/
/*===========================================================================*/
#include "otbGeometriesToGeometriesFilter.h"
#include <string>
#include <iostream>
#include "otbGeometriesSet.h"
#include "otbOGRGeometryWrapper.h"

/*===========================================================================*/
/*=========================[ Reprojection functor ]==========================*/
/*===========================================================================*/
namespace internal
{
33
struct Deleters
34
{
35 36 37 38 39 40 41 42 43 44
  void operator()(OGRCoordinateTransformation *p) const { OCTDestroyCoordinateTransformation(p); }
  void operator()(OGRSpatialReference *p) const
    {
#if GDAL_VERSION_NUM >= 1700
    OGRSpatialReference::DestroySpatialReference(p);
#else
#warning the following resource release may crash, please update your version of GDAL
    delete p; // note there is no garanty
#endif
    }
45 46 47 48 49 50 51 52
};
} // internal namespace

struct ReprojectTransformationFunctor
{
  typedef OGRGeometry TransformedElementType;

  otb::ogr::UniqueGeometryPtr operator()(OGRGeometry const* in) const;
53
  void SetTransformation(OGRSpatialReference & in, OGRSpatialReference & out)
54 55 56
    {
    m_reprojector.reset(OGRCreateCoordinateTransformation(&in, &out));
    }
57 58

private:
59
  boost::interprocess::unique_ptr<OGRCoordinateTransformation, internal::Deleters> m_reprojector;
60 61 62 63 64 65 66 67 68
};

otb::ogr::UniqueGeometryPtr
ReprojectTransformationFunctor::operator()(OGRGeometry const* in) const
{
  otb::ogr::UniqueGeometryPtr out(in ? in->clone() : 0);
  if (out)
    {
    const OGRErr err = out->transform(m_reprojector.get());
69
    if (err != OGRERR_NONE)
70 71 72 73
      {
      itkGenericExceptionMacro(<< "Cannot reproject a geometry: " << CPLGetLastErrorMsg());
      }
    }
74
  return otb::move(out);
75 76
}

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
/*===========================================================================*/
/*==========================[ reprojection filter ]==========================*/
/*===========================================================================*/
class MyReprojectionFilter : public otb::DefaultGeometriesToGeometriesFilter<ReprojectTransformationFunctor>
{
public:
  typedef MyReprojectionFilter            Self;
  typedef otb::DefaultGeometriesToGeometriesFilter<ReprojectTransformationFunctor> Superclass;
  typedef itk::SmartPointer<Self>                                Pointer;
  typedef itk::SmartPointer<const Self>                          ConstPointer;

  itkNewMacro(Self);
  itkTypeMacro(MyReprojectionFilter, otb::DefaultGeometriesToGeometriesFilter);

  void SetOutputSpatialReference(std::string const& srs);

  void UpdateTransformation(OGRSpatialReference &in)
    {
    (*this)->SetTransformation(in, *m_osr);
    }

protected:
  /** Default constructor. */
  MyReprojectionFilter(){ }
  /** Destructor. */
102
  ~MyReprojectionFilter() override{ }
103 104 105

private:

106
  OGRSpatialReference*     DoDefineNewLayerSpatialReference(otb::ogr::Layer const& itkNotUsed(source)) const override
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
    {
    return m_osr.get();
    }

  boost::interprocess::unique_ptr<OGRSpatialReference,         internal::Deleters> m_osr;
};


void
MyReprojectionFilter::SetOutputSpatialReference(std::string const& srs)
{
  OGRSpatialReference outSrs;

  char *name = const_cast<char *>(srs.c_str());

  OGRErr ret = outSrs.importFromWkt(&name);
  if (ret)
    {
    ret = outSrs.importFromProj4(name);
    }
  if (ret)
    {
    ret = outSrs.importFromEPSG(atoi(name));
    }

  if (ret)
    {
    std::cout << "Error : spatial reference system not recognized ("<<srs<<")" << std::endl;
    }

  m_osr .reset (outSrs.Clone());
}

140 141 142
/*===========================================================================*/
/*=================================[ main ]==================================*/
/*===========================================================================*/
143
typedef MyReprojectionFilter FilterType;
144

145 146 147 148 149 150 151 152 153
struct Options
{
  Options(int argc, char **argv)
    {
    workingInplace = true;
    outputIsStdout = false;

    for (int a=1; a!=argc; ++a)
      {
154 155 156 157 158 159 160 161
      if (!strcmp(argv[a], "-srs"))
        {
        if (a+1 < argc)
          {
          ++a;
          srs = argv[a];
          }
        }
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
      else if (!strcmp(argv[a], "-"))
        outputIsStdout = true;
      else if (inputFile.empty())
        inputFile = argv[a];
      else if (outputFile.empty())
        {
        outputFile = argv[a];
        workingInplace = false;
        }
      else
        throw std::runtime_error(usage(argv[0])+"Too many (unexpected) arguments");
      }

    if (!outputFile.empty() && outputIsStdout)
      throw std::runtime_error(usage(argv[0])+"An output file ("+outputFile+") has already been set, cannot dump to stdout");
177
    if (inputFile.empty() || srs.empty())
178 179 180 181
      throw std::runtime_error(usage(argv[0])+"Not enough parameters");
    }
  static std::string usage(std::string const& progname)
    {
182
    return progname + " <inputGeometriesFile> [<outputGeometriesFile>] -srs <id>\n";
183 184 185
    }
  std::string inputFile;
  std::string outputFile;
186
  std::string srs;
187 188 189 190
  bool        workingInplace;
  bool        outputIsStdout;
};

191 192 193 194 195 196 197 198
int main (int argc, char **argv)
{
  if (argc < 2)
    {
    return EXIT_FAILURE;
    }
  try
    {
199
    const Options options(argc, argv);
200 201

    otb::ogr::DataSource::Pointer input = otb::ogr::DataSource::New(
202
      options.inputFile,
203
      options.workingInplace ? otb::ogr::DataSource::Modes::Update_LayerOverwrite : otb::ogr::DataSource::Modes::Read);
204 205

    otb::ogr::DataSource::Pointer output
206 207
      = options.workingInplace ? input
      : options.outputIsStdout ? 0
208
      : otb::ogr::DataSource::New( options.outputFile, otb::ogr::DataSource::Modes::Update_LayerCreateOnly);
209
    std::cout << "input: " << otb::ogr::version_proxy::GetFileListAsStringVector(&input->ogr())[0] << " should be: " << options.inputFile << "\n";
210 211
    if (output)
      {
212
      std::cout << "output: " << otb::ogr::version_proxy::GetFileListAsStringVector(&output->ogr())[0] << " should be: " << options.outputFile << "\n";
213 214 215 216
      }
    // std::cout << "\n";

    FilterType::Pointer filter = FilterType::New();
217
    // TODO: this make no sense for in-place transformations
218 219
    filter->SetOutputSpatialReference(options.srs);
    filter->UpdateTransformation( *(const_cast<OGRSpatialReference*>(input->GetLayer(0).GetSpatialRef())) );
220 221

    if (!options.workingInplace)
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
      {
      otb::GeometriesSet::Pointer in_set = otb::GeometriesSet::New(input);
      filter->SetInput(in_set);
      }
    if (output)
      {
      otb::GeometriesSet::Pointer out_set = otb::GeometriesSet::New(output);
      filter->SetOutput(out_set);
      out_set->Update();
      }
    else
      {
      filter->UpdateOutputData(0);
      assert(filter->GetOutput() && "output not set");
      filter->GetOutput()->Print(std::cout, 0);
      }

    return EXIT_SUCCESS;
    }
  catch (std::exception const& e)
    {
    std::cerr << e.what() << "\n";
    }
  return EXIT_FAILURE;
}