pktools 2.6.7
Processing Kernel for geospatial data
ImgWriterOgr.cc
1/**********************************************************************
2ImgWriterOgr.cc: class to write vector files using OGR API library
3Copyright (C) 2008-2012 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 <fstream>
22#include <sstream>
23#include "ImgReaderOgr.h"
24#include "ImgWriterOgr.h"
25#include "ImgReaderGdal.h"
26#include "cpl_string.h"
27//---------------------------------------------------------------------------
28ImgWriterOgr::ImgWriterOgr(void)
29{}
30
31ImgWriterOgr::~ImgWriterOgr(void)
32{
33}
34
35ImgWriterOgr::ImgWriterOgr(const std::string& filename, const std::string& imageType)
36{
37 open(filename,imageType);
38}
39
40ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderOgr)
41{
42 m_filename=filename;
43 setCodec(imgReaderOgr.getDriver());
44 int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
45 for(int ilayer=0;ilayer<nlayer;++ilayer){
46 std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
47 createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
48 copyFields(imgReaderOgr,ilayer,ilayer);
49 }
50}
51
52ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderOgr, bool copyData)
53{
54 CPLErrorReset();
55 m_filename=filename;
56 setCodec(imgReaderOgr.getDriver());
57 int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
58 for(int ilayer=0;ilayer<nlayer;++ilayer){
59 std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
60 createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
61 copyFields(imgReaderOgr,ilayer,ilayer);
62 if(copyData){
63 OGRFeature *poFeature;
64 while(true){// (poFeature = imgReaderOgr.getLayer()->GetNextFeature()) != NULL ){
65 poFeature = imgReaderOgr.getLayer(ilayer)->GetNextFeature();
66 if( poFeature == NULL )
67 break;
68 OGRFeature *poDstFeature = NULL;
69
70 poDstFeature=createFeature(ilayer);
71 //todo: check here if SetFrom works (experienced segmentation fault)
72 if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
73 const char* fmt;
74 std::string errorString="Unable to translate feature %d from layer %s.\n";
75 fmt=errorString.c_str();
76 CPLError( CE_Failure, CPLE_AppDefined,
77 fmt,
78 poFeature->GetFID(), getLayerName().c_str() );
79 // CPLError( CE_Failure, CPLE_AppDefined,
80 // "Unable to translate feature %d from layer %s.\n",
81 // poFeature->GetFID(), getLayerName().c_str() );
82
83 OGRFeature::DestroyFeature( poFeature );
84 OGRFeature::DestroyFeature( poDstFeature );
85 }
86 poDstFeature->SetFID( poFeature->GetFID() );
87 OGRFeature::DestroyFeature( poFeature );
88
89 CPLErrorReset();
90 if(createFeature( poDstFeature,ilayer ) != OGRERR_NONE){
91 const char* fmt;
92 std::string errorString="Unable to translate feature %d from layer %s.\n";
93 fmt=errorString.c_str();
94 CPLError( CE_Failure, CPLE_AppDefined,
95 fmt,
96 poFeature->GetFID(), getLayerName().c_str() );
97 OGRFeature::DestroyFeature( poDstFeature );
98 }
99 OGRFeature::DestroyFeature( poDstFeature );
100 }
101 }
102 }
103}
104
105//---------------------------------------------------------------------------
106
107void ImgWriterOgr::open(const std::string& filename, ImgReaderOgr& imgReaderOgr)
108{
109 m_filename=filename;
110 setCodec(imgReaderOgr.getDriver());
111 int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
112 for(int ilayer=0;ilayer<nlayer;++ilayer){
113 std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
114 createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
115 copyFields(imgReaderOgr,ilayer,ilayer);
116 }
117}
118
119void ImgWriterOgr::open(const std::string& filename, const std::string& imageType)
120{
121 m_filename = filename;
122 setCodec(imageType);
123}
124
125//---------------------------------------------------------------------------
126void ImgWriterOgr::close(void)
127{
128#if GDAL_VERSION_MAJOR < 2
129 OGRDataSource::DestroyDataSource(m_datasource);
130#else
131 GDALClose(m_datasource);
132#endif
133}
134
135//---------------------------------------------------------------------------
136void ImgWriterOgr::setCodec(const std::string& imageType){
137#if GDAL_VERSION_MAJOR < 2
138 //register the drivers
139 OGRRegisterAll();
140 //fetch the OGR file driver
141 OGRSFDriver *poDriver;
142 poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(imageType.c_str());
143#else
144 //register the drivers
145 GDALAllRegister();
146 GDALDriver *poDriver;
147 poDriver = GetGDALDriverManager()->GetDriverByName(imageType.c_str());
148#endif
149 if( poDriver == NULL ){
150 std::string errorString="FileOpenError";
151 throw(errorString);
152 }
153#if GDAL_VERSION_MAJOR < 2
154 m_datasource = OGRSFDriverRegistrar::Open( m_filename.c_str(), TRUE );
155#else
156 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_UPDATE||GDAL_OF_VECTOR, NULL, NULL, NULL);
157#endif
158 if( m_datasource == NULL ){
159#if GDAL_VERSION_MAJOR < 2
160 m_datasource = OGRSFDriverRegistrar::Open( m_filename.c_str(), FALSE );
161#else
162 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_READONLY||GDAL_OF_VECTOR, NULL, NULL, NULL);
163#endif
164 if ( m_datasource != NULL){// we can only open in not update mode
165 std::string errorString="Update mode not supported, delete output dataset first";
166 throw(errorString);
167#if GDAL_VERSION_MAJOR < 2
168 OGRDataSource::DestroyDataSource(m_datasource);
169#else
170 GDALClose(m_datasource);
171#endif
172 m_datasource = NULL;
173 }
174 else //create the data source
175#if GDAL_VERSION_MAJOR < 2
176 m_datasource=poDriver->CreateDataSource(m_filename.c_str(),NULL);
177#else
178 m_datasource=poDriver->Create(m_filename.c_str(),0,0,0,GDT_Unknown,NULL);
179#endif
180 }
181 else{//datasets exists, always overwrite all layers (no update append for now)
182 int nLayerCount = m_datasource->GetLayerCount();
183 for(int iLayer = 0; iLayer < nLayerCount; ++iLayer){
184 if(m_datasource->DeleteLayer(iLayer)!=OGRERR_NONE){
185 std::string errorstring="DeleteLayer() failed when overwrite requested";
186 throw(errorstring);
187 }
188 }
189 }
190 if(m_datasource==NULL){
191 std::string errorString="Creation of output file failed";
192 throw(errorString);
193 }
194}
195
196#if GDAL_VERSION_MAJOR < 2
197void ImgWriterOgr::setCodec(OGRSFDriver *poDriver){
198 OGRRegisterAll();
199#else
200void ImgWriterOgr::setCodec(GDALDriver *poDriver){
201 GDALAllRegister();
202#endif
203 if( poDriver == NULL ){
204 std::string errorString="FileOpenError";
205 throw(errorString);
206 }
207#if GDAL_VERSION_MAJOR < 2
208 m_datasource = OGRSFDriverRegistrar::Open( m_filename.c_str(), TRUE );
209#else
210 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_UPDATE||GDAL_OF_VECTOR, NULL, NULL, NULL);
211#endif
212 if( m_datasource == NULL ){
213#if GDAL_VERSION_MAJOR < 2
214 m_datasource = OGRSFDriverRegistrar::Open( m_filename.c_str(), FALSE );
215#else
216 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_READONLY||GDAL_OF_VECTOR, NULL, NULL, NULL);
217#endif
218 if ( m_datasource != NULL){// we can only open in not update mode
219 std::string errorString="Update mode not supported, delete output dataset first";
220 throw(errorString);
221#if GDAL_VERSION_MAJOR < 2
222 OGRDataSource::DestroyDataSource(m_datasource);
223#else
224 GDALClose(m_datasource);
225#endif
226 m_datasource = NULL;
227 }
228 else //create the data source
229#if GDAL_VERSION_MAJOR < 2
230 m_datasource=poDriver->CreateDataSource(m_filename.c_str(),NULL);
231#else
232 m_datasource=poDriver->Create(m_filename.c_str(),0,0,0,GDT_Unknown,NULL);
233#endif
234 }
235 else{//datasets exists, always overwrite all layers (no update append for now)
236 int nLayerCount = m_datasource->GetLayerCount();
237 for(int iLayer = 0; iLayer < nLayerCount; ++iLayer){
238 if(m_datasource->DeleteLayer(iLayer)!=OGRERR_NONE){
239 std::string errorstring="DeleteLayer() failed when overwrite requested";
240 throw(errorstring);
241 }
242 }
243 }
244 if(m_datasource==NULL){
245 std::string errorString="Creation of output file failed";
246 throw(errorString);
247 }
248}
249
250// OGRLayer* ImgWriterOgr::copyLayer(OGRLayer* poSrcLayer, const std::string& layername, char** papszOptions)
251// {
252// return(m_datasource->CopyLayer(poSrcLayer, layername.c_str(),papszOptions));
253// }
254
255OGRLayer* ImgWriterOgr::createLayer(const std::string& layername, const std::string& theProjection, const OGRwkbGeometryType& eGType, char** papszOptions)
256{
257 if( !m_datasource->TestCapability( ODsCCreateLayer ) ){
258 std::string errorString="Test capability to create layer failed";
259 throw(errorString);
260 }
261 //papszOptions = CSLSetNameValue( papszOptions, "DIM", "1" );
262 //if points: use wkbPoint
263 //if no constraints on the types geometry to be written: use wkbUnknown
264 OGRLayer* poLayer;
265
266 OGRSpatialReference oSRS;
267
268 if(theProjection!=""){
269 oSRS.SetFromUserInput(theProjection.c_str());
270 poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
271 // if(theProjection.find("EPSPG:")!=std::string::npos){
272 // int epsg_code=atoi(theProjection.substr(theProjection.find_first_not_of("EPSG:")).c_str());
273 // OGRSpatialReference oSRS;
274 // oSRS.importFromEPSG(epsg_code);
275 // poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
276 // }
277 // else{
278 // OGRSpatialReference oSRS(theProjection.c_str());
279 // poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
280 // }
281 // }
282 // oSRS.importFromProj4(theProjection);
283 }
284 else
285 poLayer=m_datasource->CreateLayer( layername.c_str(), NULL, eGType,papszOptions );
286 //check if destroy is needed?!
287 CSLDestroy( papszOptions );
288 if( poLayer == NULL ){
289 std::string errorstring="Layer creation failed";
290 throw(errorstring);
291 }
292 // OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
293 return poLayer;
294}
295
296void ImgWriterOgr::createField(const std::string& fieldname, const OGRFieldType& fieldType, int theLayer)
297{
298 OGRFieldDefn oField( fieldname.c_str(), fieldType );
299 if(fieldType==OFTString)
300 oField.SetWidth(32);
301 if(theLayer<0)
302 theLayer=m_datasource->GetLayerCount()-1;//get back layer
303 if(m_datasource->GetLayer(theLayer)->CreateField( &oField ) != OGRERR_NONE ){
304 std::ostringstream es;
305 es << "Creating field " << fieldname << " failed";
306 std::string errorString=es.str();
307 throw(errorString);
308 }
309}
310
311int ImgWriterOgr::getFields(std::vector<std::string>& fields, int layer) const
312{
313 if(layer<0)
314 layer=m_datasource->GetLayerCount()-1;
315 assert(m_datasource->GetLayerCount()>layer);
316 OGRLayer *poLayer;
317 if((poLayer = m_datasource->GetLayer(layer))==NULL){
318 std::string errorstring="Could not get layer";
319 throw(errorstring);
320 }
321 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
322 fields.clear();
323 fields.resize(poFDefn->GetFieldCount());
324 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
325 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
326 fields[iField]=poFieldDefn->GetNameRef();
327 }
328 return(fields.size());
329}
330
331int ImgWriterOgr::getFields(std::vector<OGRFieldDefn*>& fields, int layer) const
332{
333 if(layer<0)
334 layer=m_datasource->GetLayerCount()-1;
335 assert(m_datasource->GetLayerCount()>layer);
336 OGRLayer *poLayer;
337 if((poLayer = m_datasource->GetLayer(layer))==NULL){
338 std::string errorstring="Could not get layer";
339 throw(errorstring);
340 }
341 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
342 fields.clear();
343 fields.resize(poFDefn->GetFieldCount());
344 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
345 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
346 fields[iField]=poFDefn->GetFieldDefn(iField);
347 }
348 assert(fields.size()==getFieldCount(layer));
349 return(fields.size());
350}
351
352void ImgWriterOgr::copyFields(const ImgReaderOgr& imgReaderOgr, int srcLayer, int targetLayer){
353 if(targetLayer<0)
354 targetLayer=m_datasource->GetLayerCount()-1;//get back layer
355 //get fields from imgReaderOgr
356 std::vector<OGRFieldDefn*> fields;
357
358 imgReaderOgr.getFields(fields,srcLayer);
359 for(int iField=0;iField<fields.size();++iField){
360 if(m_datasource->GetLayer(targetLayer)->CreateField(fields[iField]) != OGRERR_NONE ){
361 std::ostringstream es;
362 es << "Creating field " << fields[iField]->GetNameRef() << " failed";
363 std::string errorString=es.str();
364 throw(errorString);
365 }
366 }
367}
368
369void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, const std::string& theId, int layer){
370 OGRFeature *poFeature;
371 poFeature=createFeature(layer);
372 OGRPoint pt;
373 poFeature->SetField( fieldName.c_str(), theId.c_str());
374 for(std::map<std::string,double>::const_iterator mit=pointAttributes.begin();mit!=pointAttributes.end();++mit){
375 poFeature->SetField((mit->first).c_str(), mit->second);
376 }
377 pt.setX(x);
378 pt.setY(y);
379 poFeature->SetGeometry( &pt );
380 if(createFeature(poFeature,layer)!=OGRERR_NONE){
381 std::string errorString="Failed to create feature";
382 throw(errorString);
383 }
384 OGRFeature::DestroyFeature( poFeature );
385}
386
387void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, int theId, int layer){
388 OGRFeature *poFeature;
389 poFeature = createFeature(layer);
390 OGRPoint pt;
391 if(pointAttributes.size()+1!=poFeature->GetFieldCount()){
392 std::ostringstream ess;
393 ess << "Failed to add feature: " << pointAttributes.size() << " !=" << poFeature->GetFieldCount() << std::endl;
394 throw(ess.str());
395 }
396 assert(pointAttributes.size()+1==poFeature->GetFieldCount());
397 poFeature->SetField( fieldName.c_str(), theId);
398 int fid=0;
399 for(std::map<std::string,double>::const_iterator mit=pointAttributes.begin();mit!=pointAttributes.end();++mit){
400 poFeature->SetField((mit->first).c_str(),mit->second);
401 }
402 pt.setX(x);
403 pt.setY(y);
404 poFeature->SetGeometry( &pt );
405 if(createFeature(poFeature,layer)!=OGRERR_NONE){
406 std::string errorString="Failed to create feature";
407 throw(errorString);
408 }
409 OGRFeature::DestroyFeature( poFeature );
410}
411
412//add a line std::string (polygon), caller is responsible to close the line (end point=start point)
413void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer){
414 OGRFeature *poFeature;
415 poFeature = createFeature(layer);
416 poFeature->SetStyleString("PEN(c:#FF0000,w:5px)");//see also http://www.gdal.org/ogr/ogr_feature_style.html
417 poFeature->SetField( fieldName.c_str(), theId);
418 OGRLineString theLineString;
419 theLineString.setNumPoints(points.size());
420 for(int ip=0;ip<points.size();++ip)
421 theLineString.setPoint(ip,points[ip]);
422 if(poFeature->SetGeometry( &theLineString )!=OGRERR_NONE){
423 std::string errorString="Failed to set line OGRLineString as feature geometry";
424 throw(errorString);
425 }
426 if(createFeature(poFeature,layer)!=OGRERR_NONE){
427 std::string errorString="Failed to create feature";
428 throw(errorString);
429 }
430 OGRFeature::DestroyFeature( poFeature );
431}
432
433//add a ring (polygon), caller is responsible to close the line (end point=start point)?
434void ImgWriterOgr::addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer){
435 OGRFeature *poFeature;
436 poFeature = createFeature(layer);
437 poFeature->SetStyleString("PEN(c:#FF0000,w:5px)");//see also http://www.gdal.org/ogr/ogr_feature_style.html
438 poFeature->SetField( fieldName.c_str(), theId);
439 // OGRLineString theLineString;
440 // theLineString.setNumPoints(points.size());
441 OGRPolygon thePolygon;
442 OGRLinearRing theRing;
443 for(int ip=0;ip<points.size();++ip)
444 theRing.addPoint(points[ip]);
445 // theRing.addPoint(points[0]);//close the ring
446 theRing.closeRings();//redundent with previous line?
447 thePolygon.addRing(&theRing);
448 // SetSpatialFilter(&thePolygon)
449 poFeature->SetGeometry( &thePolygon );
450 if(createFeature(poFeature,layer)!=OGRERR_NONE){
451 std::string errorString="Failed to create feature";
452 throw(errorString);
453 OGRFeature::DestroyFeature( poFeature );
454 }
455 if(poFeature->SetGeometry( &thePolygon )!=OGRERR_NONE){
456 std::string errorString="Failed to set polygon as feature geometry";
457 throw(errorString);
458 }
459 OGRFeature::DestroyFeature( poFeature );
460}
461
462//add a line string (polygon), caller is responsible to close the line (end point=start point)
463void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId, int layer){
464 OGRFeature *poFeature;
465 poFeature = createFeature(layer);
466 poFeature->SetField( fieldName.c_str(), theId.c_str());
467 OGRLineString theLineString;
468 theLineString.setNumPoints(points.size());
469 for(int ip=0;ip<points.size();++ip)
470 theLineString.setPoint(ip,points[ip]);
471 if(poFeature->SetGeometry( &theLineString )!=OGRERR_NONE){
472 std::string errorString="Failed to set line OGRLineString as feature geometry";
473 throw(errorString);
474 }
475 if(createFeature(poFeature,layer)!=OGRERR_NONE){
476 std::string errorString="Failed to create feature";
477 throw(errorString);
478 }
479 OGRFeature::DestroyFeature( poFeature );
480}
481
482OGRFeature* ImgWriterOgr::createFeature(int layer){
483 return(OGRFeature::CreateFeature(m_datasource->GetLayer(layer)->GetLayerDefn()));
484}
485
486OGRErr ImgWriterOgr::createFeature(OGRFeature *theFeature,int layer){
487 return m_datasource->GetLayer(layer)->CreateFeature(theFeature);
488}
489
490int ImgWriterOgr::getFieldCount(int layer) const
491{
492 assert(m_datasource->GetLayerCount()>layer);
493 OGRLayer *poLayer;
494 if((poLayer = m_datasource->GetLayer(layer))==NULL){
495 std::string errorstring="Could not get layer";
496 throw(errorstring);
497 }
498 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
499 return(poFDefn->GetFieldCount());
500}
501
502int ImgWriterOgr::getFeatureCount(int layer) const
503{
504 return(getLayer(layer)->GetFeatureCount());
505}
506
507int ImgWriterOgr::ascii2ogr(const std::string& filename, const std::string &layername, const std::vector<std::string>& fieldName, const std::vector<OGRFieldType>& fieldType, short colX, short colY, const std::string& theProjection, const OGRwkbGeometryType& eGType, const char fs)
508{
509 char **papszOptions=NULL;
510 createLayer(layername, theProjection, eGType, papszOptions);
511 //create attribute fields that should appear on the layer. Fields must be added to the layer before any features are written. To create a field we initialize an OGRField object with the information about the field. In the case of Shapefiles, the field width and precision is significant in the creation of the output .dbf file, so we set it specifically, though generally the defaults are OK
512 int ncol=fieldName.size();
513 assert(colX>=0);
514 assert(colY>=0);
515 assert(colX<ncol+2);
516 assert(colY<ncol+2);
517 for(int ifield=0;ifield<ncol;++ifield)
518 createField(fieldName[ifield],fieldType[ifield]);
519 //create a local OGRFeature, set attributes and attach geometry before trying to write it to the layer. It is imperative that this feature be instantiated from the OGRFeatureDefn associated with the layer it will be written to.
520 //todo: try to open and catch if failure...
521 std::ifstream fpoints(filename.c_str(),std::ios::in);
522 std::string line;
523 OGRPolygon thePolygon;
524 OGRLinearRing theRing;
525 OGRPoint firstPoint;
526 OGRFeature *polyFeature;
527 if(eGType!=wkbPoint)
528 polyFeature=createFeature();
529
530
531 if(fs>' '&&fs<='~'){//field separator is a regular character (minimum ASCII code is space, maximum ASCII code is tilde)
532 std::string csvRecord;
533 while(getline(fpoints,csvRecord)){//read a line
534 OGRFeature *pointFeature;
535 if(eGType==wkbPoint)
536 pointFeature=createFeature();
537 OGRPoint thePoint;
538 bool skip=false;
539 std::istringstream csvstream(csvRecord);
540 std::string value;
541 int colId=0;
542 int fieldId=0;
543 while(getline(csvstream,value,fs)){//read a column
544 if(colId==colX)
545 thePoint.setX(atof(value.c_str()));
546 else if(colId==colY)
547 thePoint.setY(atof(value.c_str()));
548 else{
549 switch(fieldType[fieldId]){
550 case(OFTReal):
551 if(eGType==wkbPoint)
552 pointFeature->SetField(fieldId,atof(value.c_str()));
553 else if(firstPoint.IsEmpty())
554 polyFeature->SetField(fieldId,atof(value.c_str()));
555 break;
556 case(OFTInteger):
557 if(eGType==wkbPoint)
558 pointFeature->SetField(fieldId,atoi(value.c_str()));
559 else if(firstPoint.IsEmpty())
560 polyFeature->SetField(fieldId,atoi(value.c_str()));
561 break;
562 case(OFTString):
563 if(eGType==wkbPoint)
564 pointFeature->SetField(fieldId,value.c_str());
565 else if(firstPoint.IsEmpty())
566 polyFeature->SetField(fieldId,value.c_str());
567 break;
568 default:
569 break;
570 }
571 ++fieldId;
572 }
573 ++colId;
574 }
575 if(colId!=fieldId+2){
576 std::ostringstream ess;
577 ess << "Error: colId = " << colId << " is different from fieldId+2 = " << fieldId;
578 throw(ess.str());
579 }
580 if(eGType==wkbPoint){
581 pointFeature->SetGeometry( &thePoint );
582 if(createFeature(pointFeature)!=OGRERR_NONE){
583 std::string errorString="Failed to create feature";
584 throw(errorString);
585 OGRFeature::DestroyFeature( pointFeature );
586 }
587 }
588 else{
589 if(firstPoint.IsEmpty()){
590 firstPoint=thePoint;
591 }
592 theRing.addPoint(&thePoint);
593 }
594 }
595 }
596 else{//space or tab delimited fields
597 while(getline(fpoints,line)){
598 OGRFeature *pointFeature;
599 if(eGType==wkbPoint)
600 pointFeature=createFeature();
601 OGRPoint thePoint;
602 bool skip=false;
603 std::istringstream ist(line);
604 std::string value;
605 int colId=0;
606 int fieldId=0;
607 while(ist >> value){
608 if(colId==colX)
609 thePoint.setX(atof(value.c_str()));
610 else if(colId==colY)
611 thePoint.setY(atof(value.c_str()));
612 else{
613 switch(fieldType[fieldId]){
614 case(OFTReal):
615 if(eGType==wkbPoint)
616 pointFeature->SetField(fieldId,atof(value.c_str()));
617 else if(firstPoint.IsEmpty())
618 polyFeature->SetField(fieldId,atof(value.c_str()));
619 break;
620 case(OFTInteger):
621 if(eGType==wkbPoint)
622 pointFeature->SetField(fieldId,atoi(value.c_str()));
623 else if(firstPoint.IsEmpty())
624 polyFeature->SetField(fieldId,atoi(value.c_str()));
625 break;
626 case(OFTString):
627 if(eGType==wkbPoint)
628 pointFeature->SetField(fieldId,value.c_str());
629 else if(firstPoint.IsEmpty())
630 polyFeature->SetField(fieldId,value.c_str());
631 break;
632 default:
633 break;
634 }
635 ++fieldId;
636 }
637 ++colId;
638 }
639 if(colId!=fieldId+2){
640 std::ostringstream ess;
641 ess << "Error: colId = " << colId << " is different from fieldId+2 = " << fieldId;
642 throw(ess.str());
643 }
644 if(eGType==wkbPoint){
645 pointFeature->SetGeometry( &thePoint );
646 if(createFeature(pointFeature)!=OGRERR_NONE){
647 std::string errorString="Failed to create feature";
648 throw(errorString);
649 OGRFeature::DestroyFeature( pointFeature );
650 }
651 }
652 else{
653 if(firstPoint.IsEmpty()){
654 firstPoint=thePoint;
655 }
656 theRing.addPoint(&thePoint);
657 }
658 }
659 }
660 if(eGType!=wkbPoint){
661 theRing.addPoint(&firstPoint);//close the ring
662 thePolygon.addRing(&theRing);
663 // SetSpatialFilter(&thePolygon)
664 polyFeature->SetGeometry( &thePolygon );
665 if(createFeature(polyFeature)!=OGRERR_NONE){
666 std::string errorString="Failed to create feature";
667 throw(errorString);
668 OGRFeature::DestroyFeature( polyFeature );
669 }
670 }
671 return getFeatureCount();
672}
673
674int ImgWriterOgr::addData(ImgReaderGdal& imgReader, int theLayer, bool verbose)
675{
676 OGRLayer *poLayer;
677 assert(m_datasource->GetLayerCount()>theLayer);
678 if(verbose)
679 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
680 if(verbose)
681 std::cout << "get layer " << theLayer << std::endl;
682 poLayer = m_datasource->GetLayer(theLayer);
683 //start reading features from the layer
684 OGRFeature *poFeature;
685 if(verbose)
686 std::cout << "reset reading" << std::endl;
687 poLayer->ResetReading();
688 for(int iband=0;iband<imgReader.nrOfBand();++iband){
689 std::ostringstream fs;
690 fs << "band" << iband;
691 createField(fs.str(),OFTReal,theLayer);
692 }
693 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
694 int nfield=poFDefn->GetFieldCount();
695 if(verbose)
696 std::cout << "new number of fields: " << nfield << std::endl;
697 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
698 OGRGeometry *poGeometry;
699 poGeometry = poFeature->GetGeometryRef();
700 assert(poGeometry != NULL
701 && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
702 OGRPoint *poPoint = (OGRPoint *) poGeometry;
703 double x=poPoint->getX();
704 double y=poPoint->getY();
705 for(int iband=0;iband<imgReader.nrOfBand();++iband){
706 double imgData;
707 imgReader.readData(imgData,x,y,iband);
708 //todo: put imgdata in field
709 std::ostringstream fs;
710 fs << "band" << iband;
711 poFeature->SetField(fs.str().c_str(),imgData);
712 }
713 }
714 return(nfield);
715}
716
int nrOfBand(void) const
Get the number of bands of this dataset.
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