| 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 columns: 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}
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
ImgRasterGdal::getBoundingBox
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
Definition: ImgRasterGdal.cc:333
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
ImgReaderOgr
Definition: ImgReaderOgr.h: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::setImageDescription
void setImageDescription(const std::string &imageDescription)
Set the image description (only for GeoTiff format: TIFFTAG_IMAGEDESCRIPTION)
Definition: ImgWriterGdal.h:55
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 >