OTB  6.7.0
Orfeo Toolbox
otbMapFileProductWriter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbMapFileProductWriter_hxx
22 #define otbMapFileProductWriter_hxx
23 
25 #include "itksys/SystemTools.hxx"
26 #include "otbSpatialReference.h"
27 
28 namespace otb
29 {
30 
34 template <class TInputImage>
37 {
38  m_GenericRSResampler = GenericRSResamplerType::New();
39  m_TileSize = 256;
40  m_CurrentDepth = 0;
41  m_SRID = 26918;
43 
44  // Modify superclass default values, can be overridden by subclasses
45  this->SetNumberOfRequiredInputs(1);
46 }
47 
48 
52 template <class TInputImage>
55 {
56 }
57 
61 template <class TInputImage>
62 void
65 {
66  // Process object is not const-correct so the const_cast is required here
67  this->ProcessObject::SetNthInput(0,
68  const_cast< InputImageType * >( input ) );
69 }
70 
71 
72 template <class TInputImage>
73 void
75 ::SetInput( unsigned int index, const TInputImage * image )
76 {
77  // Process object is not const-correct so the const_cast is required here
78  this->ProcessObject::SetNthInput(index,
79  const_cast< TInputImage *>( image ) );
80 }
81 
85 template <class TInputImage>
89 {
90  if (this->GetNumberOfInputs() < 1)
91  {
92  return nullptr;
93  }
94 
95  return static_cast<const TInputImage * >
96  (this->ProcessObject::GetInput(0) );
97 }
98 
102 template <class TInputImage>
105 ::GetInput(unsigned int idx)
106 {
107  return static_cast< const TInputImage * >
108  (this->ProcessObject::GetInput(idx));
109 }
110 
111 
116 template <class TInputImage>
117 void
120 {
121  // Get the input Image
122  m_VectorImage = const_cast<TInputImage *>(this->GetInput());
123  m_VectorImage->UpdateOutputInformation();
125 
126  if(m_VectorImage->GetProjectionRef().empty()
127  && m_VectorImage->GetImageKeywordlist().GetSize() > 0)
128  {
129  otbMsgDevMacro( "Sensor Model detected : Reprojecting in the targer SRID" );
130  m_GenericRSResampler->SetInput(this->GetInput());
131  m_GenericRSResampler->SetOutputParametersFromMap(otb::SpatialReference::FromEPSG(m_SRID).ToWkt());
132  m_VectorImage = m_GenericRSResampler->GetOutput();
133  m_VectorImage->UpdateOutputInformation();
134  }
135 
136  // Initialize the filename, the vectordatas
137  this->Initialize();
138 
139  // Generate the mapFile
140  this->GenerateMapFile();
141 
142  // Do the tiling
143  this->Tiling();
144 }
145 
150 template <class TInputImage>
151 void
154 {
155  // check that the right extension is given : expected .map
156  if (itksys::SystemTools::GetFilenameLastExtension(m_FileName) != ".map")
157  {
158  itkExceptionMacro(<<itksys::SystemTools::GetFilenameLastExtension(m_FileName)
159  <<" is a wrong Extension FileName : Expected .map");
160  }
161 
162  // Initialize the index vectordata
163  this->InitializeVectorData();
164 }
165 
170 template <class TInputImage>
171 void
174 {
175  // Intitialize the vectordata to build the indexTile
176  m_VectorDataIndexTile = VectorDataType::New();
177  m_VectorDataIndexTile->SetProjectionRef(otb::SpatialReference::FromEPSG(m_SRID).ToWkt());
178  DataNodeType::Pointer root = m_VectorDataIndexTile->GetDataTree()->GetRoot()->Get();
179  DataNodeType::Pointer document = DataNodeType::New();
180  m_Folder = DataNodeType::New();
182 
183  document->SetNodeType(otb::DOCUMENT);
184  m_Folder->SetNodeType(otb::FOLDER);
185 
186  document->SetNodeId("DOCUMENT");
187  m_Folder->SetNodeId("FOLDER");
188 
189  m_VectorDataIndexTile->GetDataTree()->Add(document, root);
190  m_VectorDataIndexTile->GetDataTree()->Add(m_Folder, document);
191 }
192 
196 template <class TInputImage>
197 void
200 {
201  unsigned int numberOfChannel = m_VectorImage->GetNumberOfComponentsPerPixel();
202 
204  typename InputImageType::PixelType inMin(numberOfChannel), inMax(numberOfChannel),
205  outMin(numberOfChannel), outMax(numberOfChannel);
206  outMin.Fill(0);
207  outMax.Fill(255);
209 
210  // Update image base information
211  m_VectorImage->UpdateOutputInformation();
212 
213  // Get the image size
214  SizeType size;
215  size = m_VectorImage->GetLargestPossibleRegion().GetSize();
216 
217  unsigned int sizeX = size[0];
218  unsigned int sizeY = size[1];
219 
220  // Compute max depth
221  unsigned int maxDepth =
222  static_cast<unsigned int>(std::max(std::ceil(std::log(static_cast<float>(sizeX) / static_cast<float>(m_TileSize)) / std::log(2.0)),
223  std::ceil(std::log(static_cast<float>(sizeY) / static_cast<float>(m_TileSize)) / std::log(2.0))));
224 
225  // Extract size & index
226  SizeType extractSize;
227  IndexType extractIndex;
228 
229  for (unsigned int depth = 0; depth <= maxDepth; depth++)
230  {
231  // update the attribute value Current Depth
232  m_CurrentDepth = depth;
233 
234  // Compose the name of the vectordata
235  // the index shapefile filename
236  std::ostringstream tempIndexShapeName;
237  tempIndexShapeName << m_ShapeIndexPath<<"/";
238  tempIndexShapeName << itksys::SystemTools::GetFilenameWithoutExtension(m_FileName)<<"_"<<m_CurrentDepth;
239  tempIndexShapeName << "_index.shp";
240  m_IndexShapeFileName = tempIndexShapeName.str();
241 
242  // Resample image to the max Depth
243  int sampleRatioValue = 1 << (maxDepth - depth); // 2^(maxDepth - depth)
244 
245  if (sampleRatioValue > 1)
246  {
247  m_StreamingShrinkImageFilter = StreamingShrinkImageFilterType::New();
248 
249  m_StreamingShrinkImageFilter->SetShrinkFactor(sampleRatioValue);
250  m_StreamingShrinkImageFilter->SetInput(m_VectorImage);
251  m_StreamingShrinkImageFilter->GetStreamer()->SetAutomaticStrippedStreaming(0);
252  m_StreamingShrinkImageFilter->Update();
253 
254  m_VectorRescaleIntensityImageFilter = VectorRescaleIntensityImageFilterType::New();
255  m_VectorRescaleIntensityImageFilter->SetInput(m_StreamingShrinkImageFilter->GetOutput());
256  m_VectorRescaleIntensityImageFilter->SetOutputMinimum(outMin);
257  m_VectorRescaleIntensityImageFilter->SetOutputMaximum(outMax);
258 
259  if (depth == 0)
260  {
261  m_VectorRescaleIntensityImageFilter->Update();
262  inMin = m_VectorRescaleIntensityImageFilter->GetInputMinimum();
263  inMax = m_VectorRescaleIntensityImageFilter->GetInputMaximum();
264  }
265  else
266  {
267  m_VectorRescaleIntensityImageFilter->SetInputMinimum(inMin);
268  m_VectorRescaleIntensityImageFilter->SetInputMaximum(inMax);
269  m_VectorRescaleIntensityImageFilter->SetAutomaticInputMinMaxComputation(false);
270  }
271 
272  // New resample vector image
273  m_ResampleVectorImage = m_VectorRescaleIntensityImageFilter->GetOutput();
274  }
275  else
276  {
277  m_VectorRescaleIntensityImageFilter = VectorRescaleIntensityImageFilterType::New();
278  m_VectorRescaleIntensityImageFilter->SetInput(m_VectorImage);
279  m_VectorRescaleIntensityImageFilter->SetOutputMinimum(outMin);
280  m_VectorRescaleIntensityImageFilter->SetOutputMaximum(outMax);
281 
282  m_VectorRescaleIntensityImageFilter->SetInputMinimum(inMin);
283  m_VectorRescaleIntensityImageFilter->SetInputMaximum(inMax);
284  m_VectorRescaleIntensityImageFilter->SetAutomaticInputMinMaxComputation(false);
285 
286  m_ResampleVectorImage = m_VectorRescaleIntensityImageFilter->GetOutput();
287  }
288 
289  m_ResampleVectorImage->UpdateOutputInformation();
290 
291  // Get the image size
292  size = m_ResampleVectorImage->GetLargestPossibleRegion().GetSize();
293 
294  sizeX = size[0];
295  sizeY = size[1];
296 
297  unsigned int x = 0;
298  unsigned int y = 0;
299 
300  // Create directory where to store generated tiles
301  // Do it once here outside of the loop
302  std::ostringstream path;
303  path << m_ShapeIndexPath<<"/tiles";
304 
305  if (!itksys::SystemTools::MakeDirectory(path.str()))
306  {
307  itkExceptionMacro(<< "Error while creating cache directory" << path.str());
308  }
309 
310  // Tiling resample image
311  for (unsigned int tx = 0; tx < sizeX; tx += m_TileSize)
312  {
313  for (unsigned int ty = 0; ty < sizeY; ty += m_TileSize)
314  {
315  if ((tx + m_TileSize) >= sizeX)
316  {
317  extractIndex[0] = tx;
318  extractSize[0] = sizeX - tx;
319  }
320  else
321  {
322  extractIndex[0] = tx;
323  extractSize[0] = m_TileSize;
324  }
325 
326  if ((ty + m_TileSize) >= sizeY)
327  {
328  extractIndex[1] = ty;
329  extractSize[1] = sizeY - ty;
330  }
331  else
332  {
333  extractIndex[1] = ty;
334  extractSize[1] = m_TileSize;
335  }
336 
337  // Generate Tile filename
338  std::ostringstream ossFileName;
339  ossFileName << path.str()<<"/";
340  ossFileName << "tile_"<<m_CurrentDepth<<"_";
341  ossFileName << x<<"_"<<y<<".tif";
342 
343  // Extract ROI
344  m_VectorImageExtractROIFilter = VectorImageExtractROIFilterType::New();
345 
346  // Set extract roi parameters
347  m_VectorImageExtractROIFilter->SetStartX(extractIndex[0]);
348  m_VectorImageExtractROIFilter->SetStartY(extractIndex[1]);
349  m_VectorImageExtractROIFilter->SetSizeX(extractSize[0]);
350  m_VectorImageExtractROIFilter->SetSizeY(extractSize[1]);
351 
352  // Set Channel to extract
353  m_VectorImageExtractROIFilter->SetChannel(1);
354  m_VectorImageExtractROIFilter->SetChannel(2);
355  m_VectorImageExtractROIFilter->SetChannel(3);
356 
357  // Set extract roi input
358  m_VectorImageExtractROIFilter->SetInput(m_ResampleVectorImage);
359 
360  // Configure writer
361  m_VectorWriter = VectorWriterType::New();
362  m_VectorWriter->SetFileName(ossFileName.str());
363  m_VectorWriter->SetInput(m_VectorImageExtractROIFilter->GetOutput());
364  m_VectorWriter->Update();
365 
367  // Search Lat/Lon box
368 
369  // Initialize the transform to be used
370  //typename TransformType::Pointer transform =
371  //TransformType::New();
372 
373  m_Transform = TransformType::New();
374  m_Transform->SetInputProjectionRef(m_GenericRSResampler->GetOutputProjectionRef());
375  m_Transform->SetOutputProjectionRef(otb::SpatialReference::FromEPSG(m_SRID).ToWkt());
376  m_Transform->InstantiateTransform();
377 
378  InputPointType inputPoint;
379  SizeType sizeTile;
380  sizeTile = extractSize;
381 
383  // Compute upper left corner
384  itk::ContinuousIndex<double, 2> indexTile(extractIndex);
385  indexTile[0] += -0.5;
386  indexTile[1] += -0.5;
387  m_ResampleVectorImage->TransformContinuousIndexToPhysicalPoint(indexTile, inputPoint);
388  OutputPointType upperLeftCorner = m_Transform->TransformPoint(inputPoint);
389  //std::cout <<"indexTile "<< indexTile <<" --> input Point "<< inputPoint << " upperLeftCorner "<< upperLeftCorner << std::endl;
391 
392  // Compute lower left corner
393  indexTile[1] += static_cast<double>(sizeTile[1]);
394  m_ResampleVectorImage->TransformContinuousIndexToPhysicalPoint(indexTile, inputPoint);
395  OutputPointType lowerLeftCorner = m_Transform->TransformPoint(inputPoint);
396  //std::cout <<"indexTile "<< indexTile <<" --> input Point "<< inputPoint << " lowerLeftCorner "<< lowerLeftCorner << std::endl;
397 
398  // Compute lower right corner
399  indexTile[0] += static_cast<double>(sizeTile[0]);
400  m_ResampleVectorImage->TransformContinuousIndexToPhysicalPoint(indexTile, inputPoint);
401  OutputPointType lowerRightCorner = m_Transform->TransformPoint(inputPoint);
402  //std::cout <<"indexTile "<< indexTile <<" --> input Point "<< inputPoint << " lowerRightCorner "<< lowerRightCorner << std::endl;
403 
404  // Compute upper right corner
405  indexTile[1] -= static_cast<double>(sizeTile[1]);
406  m_ResampleVectorImage->TransformContinuousIndexToPhysicalPoint(indexTile, inputPoint);
407  OutputPointType upperRightCorner = m_Transform->TransformPoint(inputPoint);
408  //std::cout <<"indexTile "<< indexTile <<" --> input Point "<< inputPoint << " upperRightCorner "<< upperRightCorner << std::endl;
409 
410  // Build The indexTile
411  this->AddBBoxToIndexTile(lowerLeftCorner,
412  lowerRightCorner,
413  upperRightCorner,
414  upperLeftCorner, x, y);
415 
417  y++;
418  }
419  x++;
420  y = 0;
421  }
422 
423  // Add the layer in the mapfile
424  this->AddLayer();
425 
426  // Write the IndexTile
427  typename VectorDataFileWriterType::Pointer writer = VectorDataFileWriterType::New();
428  writer->SetFileName(m_IndexShapeFileName);
429  writer->SetInput(m_VectorDataIndexTile);
430  writer->Update();
431 
432  this->InitializeVectorData();
433  }
434 
435  // Close the file
436  m_File <<"END" << std::endl;
437  m_File.close();
438 }
439 
444 template <class TInputImage>
445 void
448  OutputPointType lowerRightCorner,
449  OutputPointType upperRightCorner,
450  OutputPointType upperLeftCorner,
451  unsigned int x, unsigned int y)
452 {
453  // From PointType to VertexType
454  VertexType pLL, pLR, pUR, pUL;
455 
456  pLL[0]=lowerLeftCorner[0];
457  pLL[1]=lowerLeftCorner[1];
458 
459  pLR[0]=lowerRightCorner[0];
460  pLR[1]=lowerRightCorner[1];
461 
462  pUR[0]=upperRightCorner[0];
463  pUR[1]=upperRightCorner[1];
464 
465  pUL[0]=upperLeftCorner[0];
466  pUL[1]=upperLeftCorner[1];
467 
468  // Form a polygon with the vertices
469  PolygonType::Pointer poly = PolygonType::New();
470  poly->AddVertex(pLL);
471  poly->AddVertex(pLR);
472  poly->AddVertex(pUR);
473  poly->AddVertex(pUL);
474 
475  // Add the polygon to the datanode
476  m_Polygon = DataNodeType::New();
477  m_Polygon->SetNodeId("FEATURE_POLYGON");
478  m_Polygon->SetNodeType(otb::FEATURE_POLYGON);
479  m_Polygon->SetPolygonExteriorRing(poly);
480 
481  std::ostringstream oss;
482  oss << "tiles/tile_";
483  oss <<m_CurrentDepth <<"_"<< x <<"_"<< y <<".tif";
484 
485  // Add the field "LOCATION" used in mapserver clients
486  // to get the path to tiles
487  m_Polygon->SetFieldAsString("LOCATION", oss.str());
488 
489  // Add the to vectordata
490  m_VectorDataIndexTile->GetDataTree()->Add(m_Polygon, m_Folder);
491 }
492 
496 template <class TInputImage>
497 void
500 {
501  // create the directory where to store the filemap
502  // if it does not exists
503  std::ostringstream path;
504  path << itksys::SystemTools::GetFilenamePath(m_FileName);
506 
507  if (!itksys::SystemTools::MakeDirectory(path.str()))
508  {
509  itkExceptionMacro(<< "Error while creating cache directory" << path.str());
510  }
511 
512  // Create a mapfile
513  m_File.open(m_FileName);
514  m_File << std::fixed << std::setprecision(6);
515 
516  // Get the name of the layer
517  std::ostringstream tempIndexShapeName;
518  tempIndexShapeName << itksys::SystemTools::GetFilenameWithoutExtension(m_FileName);
519 
520  m_File <<"MAP" << std::endl;
521  m_File <<"\tNAME "<<tempIndexShapeName.str()<< std::endl;
522  m_File <<"\t# Map image size" << std::endl;
523  m_File <<"\tSIZE "<< m_TileSize <<" "<< m_TileSize << std::endl;
524  m_File <<"\tUNITS dd" << std::endl;
525  m_File <<"\tSHAPEPATH '"<<m_ShapeIndexPath<<"'"<<std::endl;
526  m_File <<"\tEXTENT -180 -90 180 90" << std::endl; //TODO : get the
527  //extent from the
528  //vector data
529  m_File <<"\tPROJECTION" << std::endl;
530  m_File <<"\t \"init=epsg:"<<m_SRID<<"\"" << std::endl;
531  m_File <<"\tEND" << std::endl;
532 
533  m_File <<"\t# Background color for the map canvas -- change as desired" << std::endl;
534  m_File <<"\tIMAGECOLOR 192 192 192" << std::endl;
535  m_File <<"\tIMAGEQUALITY 95" << std::endl;
536  m_File <<"\tIMAGETYPE PNG" << std::endl;
537  m_File <<"\tOUTPUTFORMAT" << std::endl;
538  m_File <<"\t\tNAME PNG" << std::endl;
539  m_File <<"\t\tDRIVER 'GD/PNG'" << std::endl;
540  m_File <<"\t\tMIMETYPE 'image/png'" << std::endl;
541  m_File <<"\t\tIMAGEMODE RGB" << std::endl;
542  m_File <<"\t\tFORMATOPTION INTERLACE=OFF" << std::endl;
543  m_File <<"\t\tEXTENSION 'png'" << std::endl;
544  m_File <<"\tEND" << std::endl;
545 
546 
547  m_File <<"\t# Web interface definition. Only the template parameter" << std::endl;
548  m_File <<"\t# is required to display a map. See MapServer documentation" << std::endl;
549  m_File <<"\tWEB" << std::endl;
550  m_File <<"\t\t# Set IMAGEPATH to the path where MapServer should" << std::endl;
551  m_File <<"\t\t# write its output." << std::endl;
552  //m_File <<"#IMAGEPATH \'D:\OSGeo4W_bis2/tmp/ms_tmp/\'" << std::endl;
553 
554  m_File <<"\t\t# Set IMAGEURL to the url that points to IMAGEPATH" << std::endl;
555  m_File <<"\t\t# as defined in your web server configuration" << std::endl;
556  m_File <<"\t\t#IMAGEURL '/ms_tmp/'" << std::endl;
557 
558  m_File <<"\t\t# WMS server settings" << std::endl;
559  m_File <<"\t\t# NOTE : the user must change the path to the mapserver executable in the "<<std::endl;
560  m_File <<"\t\t# wms_onlineresource field"<<std::endl;
561  m_File <<"\t\tMETADATA" << std::endl;
562  m_File <<"\t\t 'wms_title' 'Level0'" << std::endl;
563  m_File <<"\t\t \'wms_onlineresource\' \'"<<m_CGIPath<<"?map="<<m_FileName<<"&\'" << std::endl;
564  m_File <<"\t\t \'wms_srs\' \'EPSG:"<<m_SRID<<" EPSG:900913\'" << std::endl;
565  m_File <<"\t\tEND" << std::endl;
566  m_File <<"\tEND" << std::endl;
567 
568 
569 }
570 
571 
572 template <class TInputImage>
573 void
576 {
577  m_File <<"\tLAYER" << std::endl;
578  m_File <<"\t\tNAME '"<<itksys::SystemTools::GetFilenameWithoutExtension(m_IndexShapeFileName)/*tempIndexShapeName.str()*/<<"'"<< std::endl;
579  m_File <<"\t\t\tOFFSITE 0 0 0"<< std::endl;
580  m_File <<"\t\t\tTYPE RASTER" << std::endl;
581  m_File <<"\t\t\tTILEITEM 'LOCATION'" << std::endl;
582  m_File <<"\t\t\tTILEINDEX \'"<< itksys::SystemTools::GetFilenameName(m_IndexShapeFileName)<<"\'" << std::endl;
583  //file <<"\t\t\tMETADATA" << std::endl;
584  //file <<"\t\t\t 'wms_title' 'earthsat'" << std::endl;
585  //file <<"\t\t\t 'wms_name' 'earthsat'" << std::endl;
586  //file <<"\t\t\tEND" << std::endl;
587  m_File <<"\t\t\tPROCESSING \"RESAMPLE=AVERAGE\"" << std::endl;
588  m_File <<"\t\t\tSTATUS ON" << std::endl;
589  m_File <<"\t\t\tTRANSPARENCY 100" << std::endl;
590  m_File <<"\t\t\tPROJECTION" << std::endl;
591  m_File <<"\t\t\t \"init=epsg:"<<m_SRID<<"\"" << std::endl;
592  m_File <<"\t\t\tEND" << std::endl;
593  m_File <<"\t\t#MINSCALE 250000" << std::endl; // TODO : problem with
594  // the MINSCALE AND
595  // MAXSCALE ..
596  m_File <<"\tEND" << std::endl;
597 }
598 
599 template <class TInputImage>
600 void
602 ::PrintSelf(std::ostream& os, itk::Indent indent) const
603 {
604  // Call the superclass implementation
605  Superclass::PrintSelf(os, indent);
606 }
607 
608 }
609 #endif
virtual void AddBBoxToIndexTile(OutputPointType lowerLeftCorner, OutputPointType lowerRightCorner, OutputPointType upperRightCorner, OutputPointType upperLeftCorner, unsigned int x, unsigned int y)
virtual void SetInput(const InputImageType *image)
const InputImageType * GetInput(void)
void PrintSelf(std::ostream &os, itk::Indent indent) const override
InputImageType::SizeType SizeType
PolygonType::VertexType VertexType
static SpatialReference FromEPSG(unsigned int epsg)
InputImageType::IndexType IndexType
#define otbMsgDevMacro(x)
Definition: otbMacro.h:66