| 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 omit 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}
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::getProjection
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
Definition: ImgRasterGdal.cc:40
ImgRasterGdal::getGeoTransform
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
Definition: ImgRasterGdal.cc:200
ImgRasterGdal::getColorTable
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
Definition: ImgRasterGdal.cc:118
ImgRasterGdal::setGeoTransform
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
Definition: ImgRasterGdal.cc:148
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgRasterGdal::getInterleave
std::string getInterleave() const
Get the band coding (interleave)
Definition: ImgRasterGdal.cc:293
ImgRasterGdal::GDALSetNoDataValue
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...
Definition: ImgRasterGdal.h:142
ImgRasterGdal::getDataType
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
Definition: ImgRasterGdal.cc:92
ImgRasterGdal::setProjection
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
Definition: ImgRasterGdal.cc:79
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderGdal::close
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
Definition: ImgReaderGdal.cc:42
ImgReaderGdal::open
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
Definition: ImgReaderGdal.cc:36
ImgWriterGdal
Definition: ImgWriterGdal.h:38
ImgWriterGdal::close
void close(void)
Close the image.
Definition: ImgWriterGdal.cc:74
ImgWriterGdal::open
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...
Definition: ImgWriterGdal.cc:43
ImgWriterGdal::setColorTable
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)
Definition: ImgWriterGdal.cc:254
Optionpk
Definition: Optionpk.h:107
filter::Filter
Definition: Filter.h:41