22 #include "imageclasses/ImgReaderGdal.h" 23 #include "imageclasses/ImgWriterGdal.h" 24 #include "base/Optionpk.h" 28 int main(
int argc,
char *argv[])
32 Optionpk<string> output_opt(
"o",
"output",
"output image file containing ndvi");
33 Optionpk<short> band_opt(
"b",
"band",
"Bands to be used for vegetation index (see rule option)", 0);
34 Optionpk<string> rule_opt(
"r",
"rule",
"Rule for index. ndvi (b1-b0)/(b1+b0), ndvi2 (b1-b0)/(b2+b3), gvmi (b0+0.1)-(b1+0.02))/((b0+0.1)+(b1+0.02))), vari (b1-b2)/(b1+b2-b0), osavi, mcari, tcari, diff (b1-b0), scale, ratio.",
"ndvi");
35 Optionpk<double> invalid_opt(
"t",
"invalid",
"Mask value where image is invalid.", 0);
36 Optionpk<int> nodata_opt(
"nodata",
"nodata",
"Flag value to put in image if not valid (0)", 0);
37 Optionpk<string> colorTable_opt(
"ct",
"ct",
"color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
38 Optionpk<string> description_opt(
"d",
"description",
"Set image description");
39 Optionpk<double> min_opt(
"min",
"min",
"minimum value for ndvi after scaling (set all values smaller than min to min)", 0);
40 Optionpk<double> max_opt(
"max",
"max",
"maximum value for ndvi after scaling (limit all values to max)");
41 Optionpk<double> eps_opt(
"e",
"eps",
"epsilon, contraint division by zero", 0);
42 Optionpk<double> src_scale_opt(
"src_s",
"src_scale",
"scale used for input, scale[1] is used for output: DN=scale[1]*ndvi+offset[1]", 1);
43 Optionpk<double> dst_scale_opt(
"dst_s",
"src_scale",
"scale used for output: DN=dst_s*ndvi+dst_offset", 1);
44 Optionpk<double> src_offset_opt(
"src_o",
"src_offset",
"offset used for input", 0);
45 Optionpk<double> dst_offset_opt(
"dst_o",
"dst_offset",
"offset is used for output: DN=dst_s*ndvi+dst_offset", 0);
46 Optionpk<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",
"Byte");
47 Optionpk<string> oformat_opt(
"of",
"oformat",
"Output image format (see also gdal_translate). Empty string: inherit from input image",
"GTiff");
48 Optionpk<string> option_opt(
"co",
"co",
"Creation option for output file. Multiple options can be specified.");
53 doProcess=input_opt.retrieveOption(argc,argv);
54 output_opt.retrieveOption(argc,argv);
55 band_opt.retrieveOption(argc,argv);
56 rule_opt.retrieveOption(argc,argv);
57 invalid_opt.retrieveOption(argc,argv);
58 nodata_opt.retrieveOption(argc,argv);
59 colorTable_opt.retrieveOption(argc,argv);
60 description_opt.retrieveOption(argc,argv);
61 min_opt.retrieveOption(argc,argv);
62 max_opt.retrieveOption(argc,argv);
63 eps_opt.retrieveOption(argc,argv);
64 src_scale_opt.retrieveOption(argc,argv);
65 src_offset_opt.retrieveOption(argc,argv);
66 dst_scale_opt.retrieveOption(argc,argv);
67 dst_offset_opt.retrieveOption(argc,argv);
68 otype_opt.retrieveOption(argc,argv);
69 oformat_opt.retrieveOption(argc,argv);
70 option_opt.retrieveOption(argc,argv);
71 verbose_opt.retrieveOption(argc,argv);
73 catch(
string predefinedString){
74 std::cout << predefinedString << std::endl;
78 std::cout <<
"short option -h shows basic options only, use long option --help to show all options" << std::endl;
82 if(input_opt.empty()){
83 std::cerr <<
"No input file provided (use option -i). Use --help for help information" << std::endl;
86 if(output_opt.empty()){
87 std::cerr <<
"No output file provided (use option -o). Use --help for help information" << std::endl;
92 if(rule_opt[0]==
"scale")
94 else if(rule_opt[0]==
"vari"||rule_opt[0]==
"mcari"||rule_opt[0]==
"tcari")
96 else if(rule_opt[0]==
"ndvi2")
100 while(band_opt.size()<reqBand)
101 band_opt.push_back(band_opt[0]);
103 std::cout << band_opt;
106 while(input_opt.size()<reqBand)
107 input_opt.push_back(input_opt[0]);
109 std::cout << input_opt;
111 vector<ImgReaderGdal> inputReader(reqBand);
112 for(
int ifile=0;ifile<reqBand;++ifile){
113 inputReader[ifile].open(input_opt[ifile]);
114 assert(inputReader[ifile].nrOfBand()>band_opt[ifile]);
118 cout <<
"opening output image file " << output_opt[0] << endl;
119 cout <<
"data type: " << otype_opt[0] << endl;
122 GDALDataType theType=GDT_Unknown;
124 cout <<
"possible output data types: ";
125 for(
int iType = 0; iType < GDT_TypeCount; ++iType){
127 cout <<
" " << GDALGetDataTypeName((GDALDataType)iType);
128 if( GDALGetDataTypeName((GDALDataType)iType) != NULL
129 && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
130 otype_opt[0].c_str()))
131 theType=(GDALDataType) iType;
133 if(theType==GDT_Unknown)
134 theType=inputReader[0].getDataType();
136 cout << endl <<
"Output pixel type: " << GDALGetDataTypeName(theType) << endl;
140 cout <<
"opening output image file " << output_opt[0] << endl;
142 if(option_opt.findSubstring(
"INTERLEAVE=")==option_opt.end()){
143 string theInterleave=
"INTERLEAVE=";
145 option_opt.push_back(theInterleave);
147 outputWriter.
open(output_opt[0],inputReader[0].nrOfCol(),inputReader[0].nrOfRow(),1,theType,oformat_opt[0],option_opt);
150 if(description_opt.size())
155 double ulx,uly,lrx,lry;
156 inputReader[0].getBoundingBox(ulx,uly,lrx,lry);
159 if(colorTable_opt.size()){
160 if(colorTable_opt[0]!=
"none")
163 else if (inputReader[0].getColorTable()!=NULL)
167 vector<double> lineOutput(outputWriter.
nrOfCol());
171 const char* pszMessage;
172 void* pProgressArg=NULL;
173 GDALProgressFunc pfnProgress=GDALTermProgress;
175 pfnProgress(progress,pszMessage,pProgressArg);
176 for(irow=0;irow<inputReader[0].nrOfRow();++irow){
179 if(rule_opt[0]==
"scale")
180 inputReader[0].readData(lineInput[0],GDT_Float64,irow,band_opt[0]);
181 else if(rule_opt[0]==
"vari"||rule_opt[0]==
"tcari"){
182 inputReader[0].readData(lineInput[0],GDT_Float64,irow,band_opt[0]);
183 inputReader[1].readData(lineInput[1],GDT_Float64,irow,band_opt[1]);
184 inputReader[2].readData(lineInput[2],GDT_Float64,irow,band_opt[2]);
186 else if(rule_opt[0]==
"ndvi2"){
187 inputReader[0].readData(lineInput[0],GDT_Float64,irow,band_opt[0]);
188 inputReader[1].readData(lineInput[1],GDT_Float64,irow,band_opt[1]);
189 inputReader[2].readData(lineInput[2],GDT_Float64,irow,band_opt[2]);
190 inputReader[3].readData(lineInput[3],GDT_Float64,irow,band_opt[3]);
193 inputReader[0].readData(lineInput[0],GDT_Float64,irow,band_opt[0]);
194 inputReader[1].readData(lineInput[1],GDT_Float64,irow,band_opt[1]);
197 catch(
string errorstring){
198 cerr << errorstring << endl;
201 assert(invalid_opt.size()==nodata_opt.size());
202 for(icol=0;icol<inputReader[0].nrOfCol();++icol){
203 double ndvi=min_opt[0];
204 double flagValue=nodata_opt[0];
206 for(
int iflag=0;valid&&iflag<invalid_opt.size();++iflag){
207 for(
int iband=0;iband<lineInput.size();++iband){
208 if(lineInput[iband][icol]==invalid_opt[iflag]){
209 flagValue=nodata_opt[iflag];
218 if(rule_opt[0]==
"ndvi"){
230 denom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
231 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]+(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
233 else if(rule_opt[0]==
"ndvi2"){
238 denom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
239 nom=(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]+(lineInput[3][icol]-src_offset_opt[0])/src_scale_opt[0];
241 else if(rule_opt[0]==
"gvmi"){
242 denom=((lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0]+0.1)-((lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]+0.02);
243 nom=((lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0]+0.1)+((lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]+0.02);
245 else if(rule_opt[0]==
"vari"){
246 denom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0];
247 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]+(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
249 else if(rule_opt[0]==
"osavi"){
250 denom=(1.0+0.16)*(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
251 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]+(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0]+0.16;
253 else if(rule_opt[0]==
"mcari"){
254 denom=((lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-0.2*((lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0]))*(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0];
255 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0];
257 else if(rule_opt[0]==
"tcari"){
258 denom=3*((lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]*(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-0.2*((lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0])*(lineInput[2][icol]-src_offset_opt[0])/src_scale_opt[0]);
259 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0];
261 else if(rule_opt[0]==
"diff"){
262 denom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0]-(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
265 else if(rule_opt[0]==
"scale"){
266 denom=(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
269 else if(rule_opt[0]==
"ratio"){
283 denom=(lineInput[0][icol]-src_offset_opt[0])/src_scale_opt[0];
284 nom=(lineInput[1][icol]-src_offset_opt[0])/src_scale_opt[0];
287 std::cout <<
"Error: rule " << rule_opt[0] <<
" not supported" << std::endl;
290 if(nom>eps_opt[0]||nom<-eps_opt[0])
298 lineOutput[icol]=
static_cast<int>(0.5+ndvi*dst_scale_opt[0]+dst_offset_opt[0]);
301 lineOutput[icol]=ndvi*dst_scale_opt[0]+dst_offset_opt[0];
304 if(lineOutput[icol]<min_opt[0])
305 lineOutput[icol]=min_opt[0];
306 else if(max_opt.size()){
307 if(lineOutput[icol]>max_opt[0])
308 lineOutput[icol]=max_opt[0];
312 lineOutput[icol]=flagValue;
316 outputWriter.
writeData(lineOutput,GDT_Float64,irow);
318 catch(
string errorstring){
319 cerr << errorstring << endl;
323 progress=
static_cast<float>(irow+1.0)/outputWriter.
nrOfRow();
324 pfnProgress(progress,pszMessage,pProgressArg);
326 for(
int ifile=0;ifile<inputReader.size();++ifile)
327 inputReader[ifile].close();
328 outputWriter.
close();
void setImageDescription(const std::string &imageDescription)
Set the image description (only for GeoTiff format: TIFFTAG_IMAGEDESCRIPTION)
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)
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
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.
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...
int nrOfCol(void) const
Get the number of columns of this dataset.