pktools 2.6.7
Processing Kernel for geospatial data
pkfillnodata.cc
1/**********************************************************************
2pkfillnodata.cc: program to fill holes in raster image
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include "cpl_string.h"
21#include "gdal_priv.h"
22#include "gdal.h"
23extern "C" {
24#include "gdal_alg.h"
25}
26#include <string>
27#include "base/Optionpk.h"
28
29/******************************************************************************/
66using namespace std;
67
68int main(int argc,char **argv) {
69 Optionpk<std::string> input_opt("i", "input", "Input raster dataset");
70 Optionpk<int> band_opt("b", "band", "band(s) to process (Default is -1: process all bands)");
71 Optionpk<std::string> mask_opt("m", "mask", "Mask raster dataset indicating pixels to be interpolated (zero valued) ");
72 Optionpk<std::string> output_opt("o", "output", "Output image file");
73 Optionpk<double> distance_opt("d", "distance", "Maximum number of pixels to search in all directions to find values to interpolate from", 0);
74 Optionpk<int> iteration_opt("it", "iteration", "Number of 3x3 smoothing filter passes to run (default 0)", 0);
75 Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
76
77 distance_opt.setHide(1);
78 iteration_opt.setHide(1);
79
80 bool doProcess;//stop process when program was invoked with help option (-h --help)
81 try{
82 doProcess=input_opt.retrieveOption(argc,argv);
83 band_opt.retrieveOption(argc,argv);
84 output_opt.retrieveOption(argc,argv);
85 mask_opt.retrieveOption(argc,argv);
86 distance_opt.retrieveOption(argc,argv);
87 iteration_opt.retrieveOption(argc,argv);
88 verbose_opt.retrieveOption(argc,argv);
89 }
90 catch(std::string predefinedString){
91 std::cout << predefinedString << std::endl;
92 exit(0);
93 }
94 if(!doProcess){
95 cout << endl;
96 cout << "Usage: pkfillnodata -i input.txt -m mask -o output" << endl;
97 cout << endl;
98 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
99 exit(0);//help was invoked, stop processing
100 }
101
102 assert(input_opt.size());
103 assert(mask_opt.size());
104 assert(output_opt.size());
105 GDALAllRegister();
106 GDALDataset *gds_input;
107 if(verbose_opt[0])
108 std::cout << "opening input file " << input_opt[0] << std::endl;
109 gds_input = (GDALDataset *) GDALOpen(input_opt[0].c_str(), GA_ReadOnly);
110 if(gds_input == NULL){
111 std::string errorString="FileOpenError";
112 throw(errorString);
113 }
114
115 GDALDataset *gds_mask;
116 if(verbose_opt[0])
117 std::cout << "opening mask file " << mask_opt[0] << std::endl;
118 gds_mask = (GDALDataset *) GDALOpen(mask_opt[0].c_str(), GA_ReadOnly );
119 if(gds_mask == NULL){
120 std::string errorString="FileOpenError";
121 throw(errorString);
122 }
123 GDALRasterBand *maskBand;
124 if(verbose_opt[0])
125 std::cout << "get mask raster band" << std::endl;
126 maskBand=gds_mask->GetRasterBand(1);
127
128
129 GDALDriver *poDriver;
130 poDriver = GetGDALDriverManager()->GetDriverByName(gds_input->GetDriver()->GetDescription());
131 if( poDriver == NULL ){
132 std::string errorString="FileOpenError";
133 throw(errorString);
134 }
135 if(verbose_opt[0])
136 std::cout << "copying input file to " << output_opt[0] << std::endl;
137 poDriver->CopyFiles(output_opt[0].c_str(),input_opt[0].c_str());
138 GDALDataset *gds_out;
139 gds_out=(GDALDataset *) GDALOpen(output_opt[0].c_str(), GA_Update);
140
141 if(band_opt.empty()){
142 band_opt.clear();
143 for(int iband=0;iband<gds_input->GetRasterCount();++iband)
144 band_opt.push_back(iband);
145 }
146 GDALRasterBand *targetBand;
147 for(unsigned short iband=0;iband<band_opt.size();++iband){
148 targetBand=gds_out->GetRasterBand(band_opt[iband]+1);
149 if(verbose_opt[0])
150 std::cout << "copying input file to " << output_opt[0] << std::endl;
151 double dfComplete=0.0;
152 const char* pszMessage;
153 void* pProgressArg=NULL;
154 GDALProgressFunc pfnProgress=GDALTermProgress;
155 pfnProgress(dfComplete,pszMessage,pProgressArg);
156 if(GDALFillNodata(targetBand,maskBand,distance_opt[0],0,iteration_opt[0],NULL,pfnProgress,pProgressArg)!=CE_None){
157 std::cerr << CPLGetLastErrorMsg() << std::endl;
158 exit(1);
159 }
160 else{
161 dfComplete=1.0;
162 pfnProgress(dfComplete,pszMessage,pProgressArg);
163 }
164
165 // gds_out=poDriver->CreateCopy(output_opt[0].c_str(),gds_input, FALSE,NULL,NULL,NULL);
166 // char **papszParmList=NULL;
167 // gds_out=poDriver->Create(output_opt[0].c_str(),targetBand->GetXSize(),targetBand->GetYSize(),1,targetBand->GetRasterDataType(),papszParmList);
168 // char **papszMetadata;
169 // papszMetadata = poDriver->GetMetadata();
170 // assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
171 // ostringstream compressList;
172 // gds_out->SetMetadataItem("INTERLEAVE",gds_input->GetMetadataItem( "INTERLEAVE", "IMAGE_STRUCTURE"),"IMAGE_STRUCTURE");
173 // gds_out->SetMetadataItem("COMPRESSION",gds_input->GetMetadataItem( "COMPRESSION", "IMAGE_STRUCTURE"),"IMAGE_STRUCTURE");
174 // if(gds_input->GetProjectionRef()!=NULL){
175 // gds_out->SetProjection(gds_input->GetProjectionRef());
176 // double adfGeoTransform[6];
177 // gds_input->GetGeoTransform(adfGeoTransform);
178 // gds_out->SetGeoTransform(adfGeoTransform);
179 // }
180 }
181 GDALClose(gds_input);
182 GDALClose(gds_mask);
183 GDALClose(gds_out);
184 GDALDumpOpenDatasets(stderr);
185 GDALDestroyDriverManager();
186}
187