pktools 2.6.7
Processing Kernel for geospatial data
ImgWriterGdal.cc
1/**********************************************************************
2ImgWriterGdal.cc: class to write raster files using GDAL API library
3Copyright (C) 2008-2016 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 <iostream>
21#include <iomanip>
22#include <time.h>
23#include <algorithm>
24extern "C" {
25#include "gdal_alg.h"
26}
27#include "ImgWriterGdal.h"
28
29#ifdef HAVE_CONFIG_H
30#include <config.h>
31#endif
32
34
36
37
43void ImgWriterGdal::open(const std::string& filename, const ImgReaderGdal& imgSrc, const std::vector<std::string>& options)
44{
45 m_ncol=imgSrc.nrOfCol();
46 m_nrow=imgSrc.nrOfRow();
47 m_nband=imgSrc.nrOfBand();
48 m_dataType=imgSrc.getDataType();
49 m_filename=filename;
50 m_options=options;
51 setCodec(imgSrc);
52}
53
63void ImgWriterGdal::open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::vector<std::string>& options)
64{
65 m_ncol = ncol;
66 m_nrow = nrow;
67 m_nband = nband;
68 m_dataType = dataType;
69 m_filename = filename;
70 m_options=options;
71 setCodec(imageType);
72}
73
75{
77 char **papszOptions=NULL;
78 for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
79 papszOptions=CSLAddString(papszOptions,optionIt->c_str());
80 if(papszOptions)
81 CSLDestroy(papszOptions);
82}
83
84
89 GDALAllRegister();
90 GDALDriver *poDriver;
91 poDriver = GetGDALDriverManager()->GetDriverByName(imgSrc.getDriverDescription().c_str());
92 if( poDriver == NULL ){
93 std::string errorString="FileOpenError";
94 throw(errorString);
95 }
96 char **papszMetadata;
97 papszMetadata = poDriver->GetMetadata();
98 //todo: try and catch if CREATE is not supported (as in PNG)
99 assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
100 char **papszOptions=NULL;
101 for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
102 papszOptions=CSLAddString(papszOptions,optionIt->c_str());
103
104 m_gds=poDriver->Create(m_filename.c_str(),m_ncol,m_nrow,m_nband,imgSrc.getDataType(),papszOptions);
105 double gt[6];
106 imgSrc.getGeoTransform(gt);
107 if(setGeoTransform(gt)!=CE_None)
108 std::cerr << "Warning: could not write geotransform information in " << m_filename << std::endl;
110 if(setProjection(imgSrc.getProjection())!=CE_None)
111 std::cerr << "Warning: could not write projection information in " << m_filename << std::endl;
112
113 if(m_noDataValues.size()){
114 for(int iband=0;iband<nrOfBand();++iband)
116 }
117
118 m_gds->SetMetadata(imgSrc.getMetadata() );
119 m_gds->SetMetadataItem( "TIFFTAG_DOCUMENTNAME", m_filename.c_str());
120 std::string versionString="pktools ";
121 versionString+=VERSION;
122 versionString+=" by Pieter Kempeneers";
123 m_gds->SetMetadataItem( "TIFFTAG_SOFTWARE", versionString.c_str());
124 time_t rawtime;
125 time ( &rawtime );
126
127 time_t tim=time(NULL);
128 tm *now=localtime(&tim);
129 std::ostringstream datestream;
130 //date std::string must be 20 characters long...
131 datestream << now->tm_year+1900;
132 if(now->tm_mon+1<10)
133 datestream << ":0" << now->tm_mon+1;
134 else
135 datestream << ":" << now->tm_mon+1;
136 if(now->tm_mday<10)
137 datestream << ":0" << now->tm_mday;
138 else
139 datestream << ":" << now->tm_mday;
140 if(now->tm_hour<10)
141 datestream << " 0" << now->tm_hour;
142 else
143 datestream << " " << now->tm_hour;
144 if(now->tm_min<10)
145 datestream << ":0" << now->tm_min;
146 else
147 datestream << ":" << now->tm_min;
148 if(now->tm_sec<10)
149 datestream << ":0" << now->tm_sec;
150 else
151 datestream << ":" << now->tm_sec;
152 m_gds->SetMetadataItem( "TIFFTAG_DATETIME", datestream.str().c_str());
153 if(imgSrc.getColorTable()!=NULL)
155}
156
161void ImgWriterGdal::setCodec(const std::string& imageType)
162{
163 GDALAllRegister();
164 GDALDriver *poDriver;
165 poDriver = GetGDALDriverManager()->GetDriverByName(imageType.c_str());
166 if( poDriver == NULL ){
167 std::ostringstream s;
168 s << "FileOpenError (" << imageType << ")";
169 throw(s.str());
170 }
171 char **papszMetadata;
172 papszMetadata = poDriver->GetMetadata();
173 //todo: try and catch if CREATE is not supported (as in PNG)
174 assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
175 char **papszOptions=NULL;
176 for(std::vector<std::string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
177 papszOptions=CSLAddString(papszOptions,optionIt->c_str());
178 m_gds=poDriver->Create(m_filename.c_str(),m_ncol,m_nrow,m_nband,m_dataType,papszOptions);
179 double gt[6];
180 getGeoTransform(gt);
181 if(setGeoTransform(gt)!=CE_None)
182 std::cerr << "Warning: could not write geotransform information in " << m_filename << std::endl;
183 if(setProjection(m_projection)!=CE_None)
184 std::cerr << "Warning: could not write projection information in " << m_filename << std::endl;
185 m_gds->SetMetadataItem( "TIFFTAG_DOCUMENTNAME", m_filename.c_str());
186 std::string versionString="pktools ";
187 versionString+=VERSION;
188 versionString+=" by Pieter Kempeneers";
189 m_gds->SetMetadataItem( "TIFFTAG_SOFTWARE", versionString.c_str());
190 time_t rawtime;
191 time ( &rawtime );
192
193 time_t tim=time(NULL);
194 tm *now=localtime(&tim);
195 std::ostringstream datestream;
196 //date std::string must be 20 characters long...
197 datestream << now->tm_year+1900;
198 if(now->tm_mon+1<10)
199 datestream << ":0" << now->tm_mon+1;
200 else
201 datestream << ":" << now->tm_mon+1;
202 if(now->tm_mday<10)
203 datestream << ":0" << now->tm_mday;
204 else
205 datestream << ":" << now->tm_mday;
206 if(now->tm_hour<10)
207 datestream << " 0" << now->tm_hour;
208 else
209 datestream << " " << now->tm_hour;
210 if(now->tm_min<10)
211 datestream << ":0" << now->tm_min;
212 else
213 datestream << ":" << now->tm_min;
214 if(now->tm_sec<10)
215 datestream << ":0" << now->tm_sec;
216 else
217 datestream << ":" << now->tm_sec;
218 m_gds->SetMetadataItem( "TIFFTAG_DATETIME", datestream.str().c_str());
219}
220
224void ImgWriterGdal::setMetadata(char** metadata)
225{
226 if(m_gds)
227 m_gds->SetMetadata(metadata);
228}
229
230//default projection: ETSR-LAEA
231// std::string ImgWriterGdal::setProjection(void)
232// {
233// std::string theProjection;
234// OGRSpatialReference oSRS;
235// char *pszSRS_WKT = NULL;
236// //// oSRS.importFromEPSG(3035);
237// oSRS.SetGeogCS("ETRS89","European_Terrestrial_Reference_System_1989","GRS 1980",6378137,298.2572221010042,"Greenwich",0,"degree",0.0174532925199433);
238// // cout << setprecision(16) << "major axis: " << oSRS.GetSemiMajor(NULL) << endl;//notice that major axis can be set to a different value than the default to the well known standard corresponding to the name (European_Terrestrial_Reference_System_1989), but that new value, while recognized by GetSemiMajor, will not be written in the geotiff tag!
239// oSRS.SetProjCS( "ETRS89 / ETRS-LAEA" );
240// oSRS.SetLAEA(52,10,4321000,3210000);
241// oSRS.exportToWkt( &pszSRS_WKT );
242// theProjection=pszSRS_WKT;
243// CPLFree( pszSRS_WKT );
244// assert(m_gds);
245// m_gds->SetProjection(theProjection.c_str());
246// return(theProjection);
247// }
248
249
254void ImgWriterGdal::setColorTable(const std::string& filename, int band)
255{
256 //todo: fool proof table in file (no checking currently done...)
257 std::ifstream ftable(filename.c_str(),std::ios::in);
258 std::string line;
259// poCT=new GDALColorTable();
260 GDALColorTable colorTable;
261 short nline=0;
262 while(getline(ftable,line)){
263 ++nline;
264 std::istringstream ist(line);
265 GDALColorEntry sEntry;
266 short id;
267 ist >> id >> sEntry.c1 >> sEntry.c2 >> sEntry.c3 >> sEntry.c4;
268 colorTable.SetColorEntry(id,&sEntry);
269 }
270 if(m_gds)
271 (m_gds->GetRasterBand(band+1))->SetColorTable(&colorTable);
272}
273
278void ImgWriterGdal::setColorTable(GDALColorTable* colorTable, int band)
279{
280 if(m_gds)
281 (m_gds->GetRasterBand(band+1))->SetColorTable(colorTable);
282}
283
284// //write an entire image from memory to file
285// bool ImgWriterGdal::writeData(void* pdata, const GDALDataType& dataType, int band){
286// //fetch raster band
287// GDALRasterBand *poBand;
288// if(band>=nrOfBand()+1){
289// std::ostringstream s;
290// s << "band (" << band << ") exceeds nrOfBand (" << nrOfBand() << ")";
291// throw(s.str());
292// }
293// poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
294// poBand->RasterIO(GF_Write,0,0,nrOfCol(),nrOfRow(),pdata,nrOfCol(),nrOfRow(),dataType,0,0);
295// return true;
296// }
297
314void ImgWriterGdal::rasterizeOgr(ImgReaderOgr& ogrReader, const std::vector<double>& burnValues, const std::vector<std::string>& controlOptions, const std::vector<std::string>& layernames ){
315 std::vector<int> bands;
316 if(burnValues.empty()&&controlOptions.empty()){
317 std::string errorString="Error: either burn values or control options must be provided";
318 throw(errorString);
319 }
320 for(int iband=0;iband<nrOfBand();++iband)
321 bands.push_back(iband+1);
322 std::vector<OGRLayerH> layers;
323 int nlayer=0;
324
325 std::vector<double> burnBands;//burn values for all bands in a single layer
326 std::vector<double> burnLayers;//burn values for all bands and all layers
327 if(burnValues.size()){
328 burnBands=burnValues;
329 while(burnBands.size()<nrOfBand())
330 burnBands.push_back(burnValues[0]);
331 }
332 for(int ilayer=0;ilayer<ogrReader.getLayerCount();++ilayer){
333 std::string currentLayername=ogrReader.getLayer(ilayer)->GetName();
334 if(layernames.size())
335 if(find(layernames.begin(),layernames.end(),currentLayername)==layernames.end())
336 continue;
337 std::cout << "processing layer " << currentLayername << std::endl;
338 layers.push_back((OGRLayerH)ogrReader.getLayer(ilayer));
339 ++nlayer;
340 if(burnValues.size()){
341 for(int iband=0;iband<nrOfBand();++iband)
342 burnLayers.insert(burnLayers.end(),burnBands.begin(),burnBands.end());
343 }
344 }
345 void* pTransformArg;
346 GDALProgressFunc pfnProgress=NULL;
347 void* pProgressArg=NULL;
348
349 char **coptions=NULL;
350 for(std::vector<std::string>::const_iterator optionIt=controlOptions.begin();optionIt!=controlOptions.end();++optionIt)
351 coptions=CSLAddString(coptions,optionIt->c_str());
352
353 if(controlOptions.size()){
354 if(GDALRasterizeLayers( (GDALDatasetH)m_gds,nrOfBand(),&(bands[0]),layers.size(),&(layers[0]),NULL,pTransformArg,NULL,coptions,pfnProgress,pProgressArg)!=CE_None){
355 std::string errorString(CPLGetLastErrorMsg());
356 throw(errorString);
357 }
358 }
359 else if(burnValues.size()){
360 if(GDALRasterizeLayers( (GDALDatasetH)m_gds,nrOfBand(),&(bands[0]),layers.size(),&(layers[0]),NULL,pTransformArg,&(burnLayers[0]),NULL,pfnProgress,pProgressArg)!=CE_None){
361 std::string errorString(CPLGetLastErrorMsg());
362 throw(errorString);
363 }
364 }
365 else{
366 std::string errorString="Error: either attribute fieldname or burn values must be set to rasterize vector dataset";
367 throw(errorString);
368 }
369}
virtual void close(void)
Close the image.
int nrOfRow(void) const
Get the number of rows of this dataset.
std::string getDriverDescription() const
Get the GDAL driver description of this dataset.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
int m_nband
number of bands in this dataset
int m_ncol
number of columns in this dataset
std::string m_filename
filename of this dataset
char ** getMetadata()
Get the metadata of this dataset.
GDALDataset * m_gds
instance of the GDAL dataset of this dataset
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
int m_nrow
number of rows in this dataset
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
std::vector< double > m_noDataValues
no data values for this dataset
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...
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
GDALDataType m_dataType
GDAL data type for this dataset.
void close(void)
Close the image.
void setMetadata(char **metadata)
Set specific metadata (driver specific)
virtual void setCodec(const std::string &imageType)
Register GDAL driver, setting the datatype, imagetype and some metadata.
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...
ImgWriterGdal(void)
default constructor. Image needs to be opened later with one of the open methods.
void rasterizeOgr(ImgReaderOgr &ogrReader, const std::vector< double > &burnValues, const std::vector< std::string > &controlOptions=std::vector< std::string >(), const std::vector< std::string > &layernames=std::vector< std::string >())
Rasterize an OGR vector dataset using the gdal algorithm "GDALRasterizeLayers".
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)
~ImgWriterGdal(void)
destructor