| 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}
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::isGeoRef
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
Definition: ImgRasterGdal.h:104
ImgRasterGdal::getProjection
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
Definition: ImgRasterGdal.cc:40
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::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::copyGeoTransform
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
Definition: ImgRasterGdal.cc:166
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::image2geo
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
Definition: ImgRasterGdal.cc:416
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::readData
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
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
ImgWriterGdal::writeData
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
Optionpk
Definition: Optionpk.h:107
Vector2d< double >