pktools  2.6.7
Processing Kernel for geospatial data
pkstatprofile.cc
1 /**********************************************************************
2 pkstatprofile.cc: program to calculate statistics in temporal or spectral profile
3 Copyright (C) 2008-2015 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <assert.h>
21 #include <iostream>
22 #include <string>
23 #include <fstream>
24 #include <math.h>
25 #include <sys/types.h>
26 #include <stdio.h>
27 #include "base/Optionpk.h"
28 #include "base/Vector2d.h"
29 #include "algorithms/Filter2d.h"
30 #include "algorithms/Filter.h"
31 #include "imageclasses/ImgReaderGdal.h"
32 #include "imageclasses/ImgWriterGdal.h"
33 #include "algorithms/StatFactory.h"
34 
35 /******************************************************************************/
94 using namespace std;
95 /*------------------
96  Main procedure
97  ----------------*/
98 int main(int argc,char **argv) {
99  Optionpk<std::string> input_opt("i","input","input image file");
100  Optionpk<std::string> output_opt("o", "output", "Output image file");
101  Optionpk<std::string> function_opt("f", "function", "Statistics function (mean, median, var, stdev, min, max, sum, mode (provide classes), ismin, ismax, proportion (provide classes), percentile, nvalid");
102  Optionpk<double> percentile_opt("perc","perc","Percentile value(s) used for rule percentile");
103  Optionpk<short> class_opt("class", "class", "class value(s) to use for mode, proportion");
104  Optionpk<double> nodata_opt("nodata", "nodata", "nodata value(s)");
105  Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
106  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate)","GTiff");
107  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
108  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
109  // Optionpk<short> down_opt("d", "down", "down sampling factor. Use value 1 for no downsampling). Use value n>1 for downsampling (aggregation)", 1);
110  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
111 
112  percentile_opt.setHide(1);
113  class_opt.setHide(1);
114  otype_opt.setHide(1);
115  oformat_opt.setHide(1);
116  colorTable_opt.setHide(1);
117 
118  bool doProcess;//stop process when program was invoked with help option (-h --help)
119  try{
120  doProcess=input_opt.retrieveOption(argc,argv);
121  output_opt.retrieveOption(argc,argv);
122  function_opt.retrieveOption(argc,argv);
123  percentile_opt.retrieveOption(argc,argv);
124  nodata_opt.retrieveOption(argc,argv);
125  otype_opt.retrieveOption(argc,argv);
126  oformat_opt.retrieveOption(argc,argv);
127  colorTable_opt.retrieveOption(argc,argv);
128  verbose_opt.retrieveOption(argc,argv);
129  }
130  catch(string predefinedString){
131  std::cout << predefinedString << std::endl;
132  exit(0);
133  }
134  if(!doProcess){
135  cout << endl;
136  cout << "Usage: pkstatprofile -i input -o output [-function]*" << endl;
137  cout << endl;
138  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
139  exit(0);//help was invoked, stop processing
140  }
141 
142  if(function_opt.empty()){
143  cerr << "Error: no function selected, use option -f" << endl;
144  exit(1);
145  }
146 
147  filter::Filter filter1d;
148  filter1d.setNoDataValues(nodata_opt);
149  for(int iclass=0;iclass<class_opt.size();++iclass)
150  filter1d.pushClass(class_opt[iclass]);
151 
152  filter1d.setThresholds(percentile_opt);
153 
154  ImgReaderGdal input;
155  ImgWriterGdal output;
156  if(input_opt.empty()){
157  cerr << "Error: no input file selected, use option -i" << endl;
158  exit(1);
159  }
160  if(output_opt.empty()){
161  cerr << "Error: no output file selected, use option -o" << endl;
162  exit(1);
163  }
164  try{
165  input.open(input_opt[0]);
166  }
167  catch(string errorstring){
168  cout << errorstring << endl;
169  exit(1);
170  }
171 
172  GDALDataType theType=GDT_Unknown;
173  if(verbose_opt[0])
174  cout << "possible output data types: ";
175  for(int iType = 0; iType < GDT_TypeCount; ++iType){
176  if(verbose_opt[0])
177  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
178  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
179  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
180  otype_opt[0].c_str()))
181  theType=(GDALDataType) iType;
182  }
183  if(theType==GDT_Unknown)
184  theType=input.getDataType();
185 
186  if(verbose_opt[0])
187  std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
188 
189  string imageType;//=input.getImageType();
190  if(oformat_opt.size())
191  imageType=oformat_opt[0];
192 
193  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
194  string theInterleave="INTERLEAVE=";
195  theInterleave+=input.getInterleave();
196  option_opt.push_back(theInterleave);
197  }
198  if(verbose_opt[0])
199  cout << "Calculating statistic metrics: " << function_opt.size() << endl;
200  try{
201  output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),function_opt.size(),theType,imageType,option_opt);
202  }
203  catch(string errorstring){
204  cout << errorstring << endl;
205  exit(4);
206  }
207 
208  output.setProjection(input.getProjection());
209  double gt[6];
210  input.getGeoTransform(gt);
211  output.setGeoTransform(gt);
212 
213  if(colorTable_opt.size()){
214  if(colorTable_opt[0]!="none"){
215  if(verbose_opt[0])
216  cout << "set colortable " << colorTable_opt[0] << endl;
217  assert(output.getDataType()==GDT_Byte);
218  output.setColorTable(colorTable_opt[0]);
219  }
220  }
221  else if(input.getColorTable()!=NULL)
222  output.setColorTable(input.getColorTable());
223 
224  if(nodata_opt.size()){
225  for(int iband=0;iband<output.nrOfBand();++iband)
226  output.GDALSetNoDataValue(nodata_opt[0],iband);
227  }
228 
229  try{
230  filter1d.stats(input,output,function_opt);
231  }
232  catch(string errorstring){
233  cerr << errorstring << endl;
234  }
235  input.close();
236  output.close();
237  return 0;
238 }
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98