pktools  2.6.7
Processing Kernel for geospatial data
pkreclass.cc
1 /**********************************************************************
2 pkreclass.cc: program to replace pixel values in raster image
3 Copyright (C) 2008-2014 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 <map>
22 #include "base/Optionpk.h"
23 #include "imageclasses/ImgReaderOgr.h"
24 #include "imageclasses/ImgWriterOgr.h"
25 #include "imageclasses/ImgReaderGdal.h"
26 #include "imageclasses/ImgWriterGdal.h"
27 
28 /******************************************************************************/
67 using namespace std;
68 
69 int main(int argc, char *argv[])
70 {
71  Optionpk<string> input_opt("i", "input", "Input image");
72  Optionpk<string> mask_opt("m", "mask", "Mask image(s)");
73  Optionpk<string> output_opt("o", "output", "Output mask file");
74  Optionpk<unsigned short> masknodata_opt("msknodata", "msknodata", "Mask value(s) where image has nodata. Use one value for each mask, or multiple values for a single mask.", 1);
75  Optionpk<int> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid (0)", 0);
76  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
77  Optionpk<unsigned short> band_opt("b", "band", "band index(es) to replace (other bands are copied to output)", 0);
78  Optionpk<string> type_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", "");
79  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
80  Optionpk<string> code_opt("code", "code", "Recode text file (2 colums: from to)");
81  Optionpk<string> class_opt("c", "class", "list of classes to reclass (in combination with reclass option)");
82  Optionpk<string> reclass_opt("r", "reclass", "list of recoded classes (in combination with class option)");
83  Optionpk<string> fieldname_opt("n", "fname", "field name of the shape file to be replaced", "label");
84  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
85  Optionpk<string> description_opt("d", "description", "Set image description");
86  Optionpk<short> verbose_opt("v", "verbose", "verbose", 0);
87 
88  bool doProcess;//stop process when program was invoked with help option (-h --help)
89  try{
90  doProcess=input_opt.retrieveOption(argc,argv);
91  mask_opt.retrieveOption(argc,argv);
92  masknodata_opt.retrieveOption(argc,argv);
93  nodata_opt.retrieveOption(argc,argv);
94  code_opt.retrieveOption(argc,argv);
95  class_opt.retrieveOption(argc,argv);
96  reclass_opt.retrieveOption(argc,argv);
97  colorTable_opt.retrieveOption(argc,argv);
98  output_opt.retrieveOption(argc,argv);
99  type_opt.retrieveOption(argc,argv);
100  oformat_opt.retrieveOption(argc,argv);
101  band_opt.retrieveOption(argc,argv);
102  fieldname_opt.retrieveOption(argc,argv);
103  option_opt.retrieveOption(argc,argv);
104  description_opt.retrieveOption(argc,argv);
105  verbose_opt.retrieveOption(argc,argv);
106  }
107  catch(string predefinedString){
108  std::cout << predefinedString << std::endl;
109  exit(0);
110  }
111  if(!doProcess){
112  cout << endl;
113  cout << "Usage: pkreclass -i input [-c from -r to]* -o output" << endl;
114  cout << endl;
115  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
116  exit(0);//help was invoked, stop processing
117  }
118 
119  if(input_opt.empty()){
120  std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
121  exit(0);
122  }
123  if(output_opt.empty()){
124  std::cerr << "No output file provided (use option -o). Use --help for help information" << std::endl;
125  exit(0);
126  }
127 
128  // vector<short> bandVector;
129  // for(int iband=0;iband<band_opt.size();++iband)
130  // bandVector.push_back(band_opt[iband]);
131  map<string,string> codemapString;//map with codes: codemapString[theKey(from)]=theValue(to)
132  map<double,double> codemap;//map with codes: codemap[theKey(from)]=theValue(to)
133  if(code_opt.size()){
134  if(verbose_opt[0])
135  cout << "opening code text file " << code_opt[0] << endl;
136  ifstream codefile;
137  codefile.open(code_opt[0].c_str());
138  string theKey;
139  string theValue;
140  while(codefile>>theKey){
141  codefile >> theValue;
142  codemapString[theKey]=theValue;
143  codemap[string2type<double>(theKey)]=string2type<double>(theValue);
144  }
145  codefile.close();
146  }
147  else{//use combination of class_opt and reclass_opt
148  assert(class_opt.size()==reclass_opt.size());
149  for(int iclass=0;iclass<class_opt.size();++iclass){
150  codemapString[class_opt[iclass]]=reclass_opt[iclass];
151  codemap[string2type<double>(class_opt[iclass])]=string2type<double>(reclass_opt[iclass]);
152  }
153  }
154  assert(codemapString.size());
155  assert(codemap.size());
156  //if verbose true, print the codes to screen
157  if(verbose_opt[0]){
158  map<string,string>::iterator mit;
159  cout << codemapString.size() << " codes used: " << endl;
160  for(mit=codemapString.begin();mit!=codemapString.end();++mit)
161  cout << (*mit).first << " " << (*mit).second << endl;
162  }
163  bool refIsRaster=true;
164  ImgReaderOgr ogrReader;
165  if(refIsRaster){//image file
166  ImgReaderGdal inputReader;
167  vector<ImgReaderGdal> maskReader(mask_opt.size());
168  ImgWriterGdal outputWriter;
169  if(verbose_opt[0])
170  cout << "opening input image file " << input_opt[0] << endl;
171  inputReader.open(input_opt[0]);
172  for(int imask=0;imask<mask_opt.size();++imask){
173  if(verbose_opt[0])
174  cout << "opening mask image file " << mask_opt[imask] << endl;
175  maskReader[imask].open(mask_opt[imask]);
176  }
177  if(verbose_opt[0]){
178  cout << "opening output image file " << output_opt[0] << endl;
179  cout << "data type: " << type_opt[0] << endl;
180  }
181  //create output image with user defined data type
182  GDALDataType theType=GDT_Unknown;
183  if(verbose_opt[0])
184  cout << "possible output data types: ";
185  for(int iType = 0; iType < GDT_TypeCount; ++iType){
186  if(verbose_opt[0])
187  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
188  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
189  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
190  type_opt[0].c_str()))
191  theType=(GDALDataType) iType;
192  }
193  if(theType==GDT_Unknown)
194  theType=inputReader.getDataType();
195  if(verbose_opt[0])
196  cout << endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
197  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
198  string theInterleave="INTERLEAVE=";
199  theInterleave+=inputReader.getInterleave();
200  option_opt.push_back(theInterleave);
201  }
202  outputWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),inputReader.nrOfBand(),theType,oformat_opt[0],option_opt);
203  for(int iband=0;iband<inputReader.nrOfBand();++iband)
204  outputWriter.GDALSetNoDataValue(nodata_opt[0],iband);
205  if(description_opt.size())
206  outputWriter.setImageDescription(description_opt[0]);
207 
208  if(colorTable_opt.size()){
209  if(colorTable_opt[0]!="none")
210  outputWriter.setColorTable(colorTable_opt[0]);
211  }
212  else if (inputReader.getColorTable()!=NULL)//copy colorTable from input image
213  outputWriter.setColorTable(inputReader.getColorTable());
214 
215  //if input image is georeferenced, copy projection info to output image
216  if(inputReader.isGeoRef()){
217  for(int imask=0;imask<mask_opt.size();++imask)
218  assert(maskReader[imask].isGeoRef());
219  }
220  outputWriter.copyGeoTransform(inputReader);
221  outputWriter.setProjection(inputReader.getProjection());
222  double ulx,uly,lrx,lry;
223  inputReader.getBoundingBox(ulx,uly,lrx,lry);
224  outputWriter.copyGeoTransform(inputReader);
225  assert(nodata_opt.size()==masknodata_opt.size());
226  if(verbose_opt[0]&&mask_opt.size()){
227  for(int iv=0;iv<masknodata_opt.size();++iv)
228  cout << masknodata_opt[iv] << "->" << nodata_opt[iv] << endl;
229  }
230 
231  assert(outputWriter.nrOfCol()==inputReader.nrOfCol());
232  // Vector2d<int> lineInput(inputReader.nrOfBand(),inputReader.nrOfCol());
233  Vector2d<double> lineInput(inputReader.nrOfBand(),inputReader.nrOfCol());
234  Vector2d<short> lineMask(mask_opt.size());
235  for(int imask=0;imask<mask_opt.size();++imask)
236  lineMask[imask].resize(maskReader[imask].nrOfCol());
237  Vector2d<double> lineOutput(outputWriter.nrOfBand(),outputWriter.nrOfCol());
238  int irow=0;
239  int icol=0;
240  const char* pszMessage;
241  void* pProgressArg=NULL;
242  GDALProgressFunc pfnProgress=GDALTermProgress;
243  double progress=0;
244  pfnProgress(progress,pszMessage,pProgressArg);
245  double oldRowMask=-1;
246  for(irow=0;irow<inputReader.nrOfRow();++irow){
247  //read line in lineInput buffer
248  for(int iband=0;iband<inputReader.nrOfBand();++iband){
249  try{
250  // inputReader.readData(lineInput[iband],GDT_Int32,irow,iband);
251  inputReader.readData(lineInput[iband],irow,iband);
252  }
253  catch(string errorstring){
254  cerr << errorstring << endl;
255  exit(1);
256  }
257  }
258  double x,y;//geo coordinates
259  double colMask,rowMask;//image coordinates in mask image
260  for(icol=0;icol<inputReader.nrOfCol();++icol){
261  bool masked=false;
262  if(mask_opt.size()>1){//multiple masks
263  for(int imask=0;imask<mask_opt.size();++imask){
264  inputReader.image2geo(icol,irow,x,y);
265  maskReader[imask].geo2image(x,y,colMask,rowMask);
266  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
267  assert(rowMask>=0&&rowMask<maskReader[imask].nrOfRow());
268  try{
269  maskReader[imask].readData(lineMask[imask],static_cast<int>(rowMask));
270  }
271  catch(string errorstring){
272  cerr << errorstring << endl;
273  exit(1);
274  }
275  oldRowMask=rowMask;
276  }
277  short ivalue=0;
278  if(mask_opt.size()==masknodata_opt.size())//one invalid value for each mask
279  ivalue=masknodata_opt[imask];
280  else//use same invalid value for each mask
281  ivalue=masknodata_opt[0];
282  if(lineMask[imask][colMask]==ivalue){
283  for(int iband=0;iband<inputReader.nrOfBand();++iband)
284  lineInput[iband][icol]=nodata_opt[imask];
285  masked=true;
286  break;
287  }
288  }
289  }
290  else if(mask_opt.size()){//potentially more invalid values for single mask
291  inputReader.image2geo(icol,irow,x,y);
292  maskReader[0].geo2image(x,y,colMask,rowMask);
293  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
294  assert(rowMask>=0&&rowMask<maskReader[0].nrOfRow());
295  try{
296  maskReader[0].readData(lineMask[0],static_cast<int>(rowMask));
297  }
298  catch(string errorstring){
299  cerr << errorstring << endl;
300  exit(1);
301  }
302  oldRowMask=rowMask;
303  }
304  for(int ivalue=0;ivalue<masknodata_opt.size();++ivalue){
305  assert(masknodata_opt.size()==nodata_opt.size());
306  if(lineMask[0][colMask]==masknodata_opt[ivalue]){
307  for(int iband=0;iband<inputReader.nrOfBand();++iband)
308  lineInput[iband][icol]=nodata_opt[ivalue];
309  masked=true;
310  break;
311  }
312  }
313  }
314  for(int iband=0;iband<lineOutput.size();++iband){
315  lineOutput[iband][icol]=lineInput[iband][icol];
316  if(find(band_opt.begin(),band_opt.end(),iband)!=band_opt.end()){
317  if(!masked && codemap.find(lineInput[iband][icol])!=codemap.end()){
318  double toValue=codemap[lineInput[iband][icol]];
319  lineOutput[iband][icol]=toValue;
320  }
321  }
322  }
323  }
324  //write buffer lineOutput to output file
325  try{
326  for(int iband=0;iband<outputWriter.nrOfBand();++iband)
327  outputWriter.writeData(lineOutput[iband],irow,iband);
328  }
329  catch(string errorstring){
330  cerr << errorstring << endl;
331  exit(1);
332  }
333  //progress bar
334  progress=static_cast<float>((irow+1.0)/outputWriter.nrOfRow());
335  pfnProgress(progress,pszMessage,pProgressArg);
336  }
337  inputReader.close();
338  outputWriter.close();
339  }
340 }
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.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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.
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
std::string getInterleave() const
Get the band coding (interleave)
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98