pktools  2.6.7
Processing Kernel for geospatial data
pksetmask.cc
1 /**********************************************************************
2 pksetmask.cc: program to apply mask image (set invalid values) to 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 
22 #include "imageclasses/ImgReaderGdal.h"
23 #include "imageclasses/ImgWriterGdal.h"
24 #include "base/Optionpk.h"
25 /******************************************************************************/
71 using namespace std;
72 
73 int main(int argc, char *argv[])
74 {
75  //command line options
76  Optionpk<string> input_opt("i", "input", "Input image");
77  Optionpk<string> mask_opt("m", "mask", "Mask image(s)");
78  Optionpk<string> output_opt("o", "output", "Output mask file");
79  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", "");
80  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate)","GTiff");
81  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
82  Optionpk<int> msknodata_opt("msknodata", "msknodata", "Mask value(s) where image has nodata. Use one value for each mask, or multiple values for a single mask.", 1);
83  Optionpk<short> mskband_opt("mskband", "mskband", "Mask band to read (0 indexed). Provide band for each mask.", 0);
84  Optionpk<char> operator_opt("p", "operator", "Operator: < = > !. Use operator for each msknodata option", '=');
85  Optionpk<int> nodata_opt("nodata", "nodata", "nodata value to put in image if not valid", 0);
86  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
87  Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
88 
89  otype_opt.setHide(1);
90  oformat_opt.setHide(1);
91  option_opt.setHide(1);
92  colorTable_opt.setHide(1);
93  mskband_opt.setHide(1);
94 
95  bool doProcess;//stop process when program was invoked with help option (-h --help)
96  try{
97  doProcess=input_opt.retrieveOption(argc,argv);
98  mask_opt.retrieveOption(argc,argv);
99  msknodata_opt.retrieveOption(argc,argv);
100  mskband_opt.retrieveOption(argc,argv);
101  output_opt.retrieveOption(argc,argv);
102  nodata_opt.retrieveOption(argc,argv);
103  operator_opt.retrieveOption(argc,argv);
104  otype_opt.retrieveOption(argc,argv);
105  oformat_opt.retrieveOption(argc,argv);
106  option_opt.retrieveOption(argc,argv);
107  colorTable_opt.retrieveOption(argc,argv);
108  verbose_opt.retrieveOption(argc,argv);
109  }
110  catch(string predefinedString){
111  std::cout << predefinedString << std::endl;
112  exit(0);
113  }
114  if(!doProcess){
115  cout << endl;
116  cout << "Usage: pksetmask -i input -m mask [-m mask]* [-msknodata value -nodata value]* -o output" << endl;
117  cout << endl;
118  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
119  exit(0);//help was invoked, stop processing
120  }
121 
122  if(verbose_opt[0])
123  cout << "number of mask images: " << mask_opt.size() << endl;
124 
125  //duplicate band used for mask if not explicitly provided
126  while(mskband_opt.size()<mask_opt.size())
127  mskband_opt.push_back(mskband_opt[0]);
128 
129  vector<ImgReaderGdal> maskReader(mask_opt.size());
130  for(int imask=0;imask<mask_opt.size();++imask){
131  if(verbose_opt[0])
132  cout << "opening mask image file " << mask_opt[imask] << endl;
133  maskReader[imask].open(mask_opt[imask]);
134  }
135  assert(input_opt.size());
136  if(verbose_opt[0])
137  cout << "opening input image file " << input_opt[0] << endl;
138  ImgReaderGdal inputReader;
139  inputReader.open(input_opt[0]);
140  string imageType;//=inputReader.getImageType();
141  if(oformat_opt.size())//default
142  imageType=oformat_opt[0];
143  GDALDataType theType=GDT_Unknown;
144  if(verbose_opt[0]){
145  std::cout << "Image type: " << imageType << std::endl;
146  std::cout << "possible output data types: ";
147  }
148  for(int iType = 0; iType < GDT_TypeCount; ++iType){
149  if(verbose_opt[0])
150  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
151  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
152  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
153  otype_opt[0].c_str()))
154  theType=(GDALDataType) iType;
155  }
156  if(theType==GDT_Unknown)
157  theType=inputReader.getDataType();
158 
159  assert(output_opt.size());
160  if(verbose_opt[0]){
161  std::cout << std::endl << "Output data type: " << GDALGetDataTypeName(theType) << std::endl;
162  std::cout << "opening output image for writing: " << output_opt[0] << std::endl;
163  }
164  ImgWriterGdal outputWriter;
165  try{
166  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
167  string theInterleave="INTERLEAVE=";
168  theInterleave+=inputReader.getInterleave();
169  option_opt.push_back(theInterleave);
170  }
171  outputWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),inputReader.nrOfBand(),theType,imageType,option_opt);
172  for(int iband=0;iband<inputReader.nrOfBand();++iband)
173  outputWriter.GDALSetNoDataValue(nodata_opt[0],iband);
174  outputWriter.setProjection(inputReader.getProjection());
175  outputWriter.copyGeoTransform(inputReader);
176  }
177  catch(string errorstring){
178  cout << errorstring << endl;
179  exit(1);
180  }
181  // if(verbose_opt[0])
182  // cout << "opening output image file " << output_opt[0] << endl;
183  // outputWriter.open(output_opt[0],inputReader);
184  if(colorTable_opt.size()){
185  if(colorTable_opt[0]!="none")
186  outputWriter.setColorTable(colorTable_opt[0]);
187  }
188  else if (inputReader.getColorTable()!=NULL)//copy colorTable from input image
189  outputWriter.setColorTable(inputReader.getColorTable());
190  if(inputReader.isGeoRef()){
191  for(int imask=0;imask<mask_opt.size();++imask)
192  assert(maskReader[imask].isGeoRef());
193  }
194  assert(nodata_opt.size()==msknodata_opt.size());
195  assert(operator_opt.size()==msknodata_opt.size()||operator_opt.size()==1);
196  if(verbose_opt[0]){
197  cout << " mask files selected: " << mask_opt.size() << endl;
198  for(int iv=0;iv<msknodata_opt.size();++iv){
199  char op=(operator_opt.size()==msknodata_opt.size())?operator_opt[iv]:operator_opt[0];
200  cout << op << " " << msknodata_opt[iv] << "->" << nodata_opt[iv] << endl;
201  }
202  }
203 
204  Vector2d<double> lineInput(inputReader.nrOfBand(),inputReader.nrOfCol());
205  Vector2d<double> lineOutput(outputWriter.nrOfBand(),outputWriter.nrOfCol());
206  assert(lineOutput.size()==lineInput.size());
207  assert(inputReader.nrOfCol()==outputWriter.nrOfCol());
208  // Vector2d<int> lineMask(mask_opt.size());
209  Vector2d<double> lineMask(mask_opt.size());
210  for(int imask=0;imask<mask_opt.size();++imask){
211  if(verbose_opt[0])
212  cout << "mask " << imask << " has " << maskReader[imask].nrOfCol() << " columns and " << maskReader[imask].nrOfRow() << " rows" << endl;
213  lineMask[imask].resize(maskReader[imask].nrOfCol());
214  }
215  int irow=0;
216  int icol=0;
217  const char* pszMessage;
218  void* pProgressArg=NULL;
219  GDALProgressFunc pfnProgress=GDALTermProgress;
220  float progress=0;
221  if(!verbose_opt[0])
222  pfnProgress(progress,pszMessage,pProgressArg);
223  // double oldRowMask=-1;
224  vector<double> oldRowMask(mask_opt.size());
225  for(int imask=0;imask<mask_opt.size();++imask)
226  oldRowMask[imask]=-1;
227  for(irow=0;irow<inputReader.nrOfRow();++irow){
228  //read line in lineInput buffer
229  for(int iband=0;iband<inputReader.nrOfBand();++iband){
230  try{
231  inputReader.readData(lineInput[iband],irow,iband);
232  }
233  catch(string errorstring){
234  cerr << errorstring << endl;
235  exit(1);
236  }
237  }
238  double x,y;//geo coordinates
239  double colMask,rowMask;//image coordinates in mask image
240  for(icol=0;icol<inputReader.nrOfCol();++icol){
241  if(mask_opt.size()>1){//multiple masks
242  for(int imask=0;imask<mask_opt.size();++imask){
243  inputReader.image2geo(icol,irow,x,y);
244  maskReader[imask].geo2image(x,y,colMask,rowMask);
245  colMask=static_cast<int>(colMask);
246  rowMask=static_cast<int>(rowMask);
247  bool masked=false;
248  if(rowMask>=0&&rowMask<maskReader[imask].nrOfRow()&&colMask>=0&&colMask<maskReader[imask].nrOfCol()){
249  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask[imask])){
250  assert(rowMask>=0&&rowMask<maskReader[imask].nrOfRow());
251  try{
252  // maskReader[imask].readData(lineMask[imask],static_cast<int>(rowMask));
253  maskReader[imask].readData(lineMask[imask],static_cast<int>(rowMask),mskband_opt[imask]);
254  }
255  catch(string errorstring){
256  cerr << errorstring << endl;
257  exit(1);
258  }
259  oldRowMask[imask]=rowMask;
260  }
261  }
262  else
263  continue;//no coverage in this mask
264  int ivalue=0;
265  if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
266  ivalue=msknodata_opt[imask];
267  else//use same invalid value for each mask
268  ivalue=msknodata_opt[0];
269  char op=(operator_opt.size()==mask_opt.size())?operator_opt[imask]:operator_opt[0];
270  switch(op){
271  case('='):
272  default:
273  if(lineMask[imask][colMask]==ivalue)
274  masked=true;
275  break;
276  case('<'):
277  if(lineMask[imask][colMask]<ivalue)
278  masked=true;
279  break;
280  case('>'):
281  if(lineMask[imask][colMask]>ivalue)
282  masked=true;
283  break;
284  case('!'):
285  if(lineMask[imask][colMask]!=ivalue)
286  masked=true;
287  break;
288  }
289  if(masked){
290  if(verbose_opt[0]>1)
291  cout << "image masked at (col=" << icol << ",row=" << irow <<") with mask " << mask_opt[imask] << " and value " << ivalue << endl;
292  for(int iband=0;iband<inputReader.nrOfBand();++iband){
293  if(mask_opt.size()==nodata_opt.size())//one flag value for each mask
294  lineInput[iband][icol]=nodata_opt[imask];
295  else
296  lineInput[iband][icol]=nodata_opt[0];
297  }
298  masked=false;
299  break;
300  }
301  }
302  }
303  else{//potentially more invalid values for single mask
304  inputReader.image2geo(icol,irow,x,y);
305  maskReader[0].geo2image(x,y,colMask,rowMask);
306  colMask=static_cast<int>(colMask);
307  rowMask=static_cast<int>(rowMask);
308  bool masked=false;
309  if(rowMask>=0&&rowMask<maskReader[0].nrOfRow()&&colMask>=0&&colMask<maskReader[0].nrOfCol()){
310  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask[0])){
311  assert(rowMask>=0&&rowMask<maskReader[0].nrOfRow());
312  try{
313  // maskReader[0].readData(lineMask[0],static_cast<int>(rowMask));
314  maskReader[0].readData(lineMask[0],static_cast<int>(rowMask),mskband_opt[0]);
315  }
316  catch(string errorstring){
317  cerr << errorstring << endl;
318  exit(1);
319  }
320  oldRowMask[0]=rowMask;
321  }
322  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
323  assert(msknodata_opt.size()==nodata_opt.size());
324  char op=(operator_opt.size()==msknodata_opt.size())?operator_opt[ivalue]:operator_opt[0];
325  switch(op){
326  case('='):
327  default:
328  if(lineMask[0][colMask]==msknodata_opt[ivalue])
329  masked=true;
330  break;
331  case('<'):
332  if(lineMask[0][colMask]<msknodata_opt[ivalue])
333  masked=true;
334  break;
335  case('>'):
336  if(lineMask[0][colMask]>msknodata_opt[ivalue])
337  masked=true;
338  break;
339  case('!'):
340  if(lineMask[0][colMask]!=msknodata_opt[ivalue])
341  masked=true;
342  break;
343  }
344  if(masked){
345  for(int iband=0;iband<inputReader.nrOfBand();++iband)
346  lineInput[iband][icol]=nodata_opt[ivalue];
347  masked=false;
348  break;
349  }
350  }
351  }
352  }
353  for(int iband=0;iband<lineOutput.size();++iband)
354  lineOutput[iband][icol]=lineInput[iband][icol];
355  }
356  //write buffer lineOutput to output file
357  for(int iband=0;iband<outputWriter.nrOfBand();++iband){
358  try{
359  outputWriter.writeData(lineOutput[iband],irow,iband);
360  }
361  catch(string errorstring){
362  cerr << errorstring << endl;
363  exit(1);
364  }
365  }
366  //progress bar
367  progress=static_cast<float>(irow+1.0)/outputWriter.nrOfRow();
368  pfnProgress(progress,pszMessage,pProgressArg);
369  }
370  inputReader.close();
371  for(int imask=0;imask<mask_opt.size();++imask)
372  maskReader[imask].close();
373  outputWriter.close();
374 }
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.
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 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
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...
Definition: ImgWriterGdal.h:96
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.
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