pktools 2.6.7
Processing Kernel for geospatial data
ImgReaderOgr.h
1/**********************************************************************
2ImgReaderOgr.h: class to read 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#ifndef _IMGREADEROGR_H_
21#define _IMGREADEROGR_H_
22
23#include <assert.h>
24#include <fstream>
25#include <string>
26#include <sstream>
27#include <map>
28#include <vector>
29#include <iostream>
30#include "ogrsf_frmts.h"
31#include "base/Vector2d.h"
32#include "ImgReaderGdal.h"
33
34//--------------------------------------------------------------------------
36{
37public:
38 ImgReaderOgr(void);
39 ImgReaderOgr(const std::string& filename);
40 ~ImgReaderOgr(void);
41 void open(const std::string& filename);
42 void close(void);
43 // int readData(Vector2d<double>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, int layer=0, bool pos=false);//default layer 0 and no pos information in data
44 template <typename T> int readXY(std::vector<T>& xVector, std::vector<T>& yVector, int layer=0, bool verbose=false);
45 template <typename T> int readY(std::vector<T>& yVector, int layer=0, bool verbose=false);
46 template <typename T> int readData(std::vector<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, OGRFeature *poFeature, int layer=0, bool pos=false, bool verbose=false);
47 template <typename T> int readData(std::vector<T>& data, const OGRFieldType& fieldType, const std::string& theField, int layer=0, bool verbose=false);
48 template <typename T> int readData(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
49 template <typename T> int readData(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
50 template <typename T> int readData(std::map<std::string,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer=0, bool pos=false, bool verbose=false);//default layer 0 and no pos information in data
51 unsigned int readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
52 std::vector<std::string>& fields,
53 const std::vector<unsigned short>& bands,
54 const std::string& label,
55 const std::vector<std::string>& layers,
56 int verbose=false);
57
58 unsigned int readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
59 std::vector<std::string>& fields,
60 double start,
61 double end,
62 const std::string& label,
63 const std::vector<std::string>& layers,
64 int verbose=false);
65
66 // void shape2ascii(std::ostream& theOstream, const std::string& pointname, int layer=0, bool verbose=false);
67 unsigned long int getFeatureCount(int layer=0) const;
68 int getFieldCount(int layer=0) const;
69 OGRLayer* getLayer(int layer=0){return m_datasource->GetLayer(layer);};
70 std::string getProjection(int layer=0) const;
71 OGRwkbGeometryType getGeometryType(int layer=0) const;
72 std::string getLayerName(int layer=0){return m_datasource->GetLayer(layer)->GetLayerDefn()->GetName();};
73 // int getLayer(int layer=0) const;
74 int getFields(std::vector<std::string>& fields, int layer=0) const;
75 int getFields(std::vector<OGRFieldDefn*>& fields, int layer=0) const;
76#if GDAL_VERSION_MAJOR < 2
77 OGRDataSource* getDataSource(void) {return m_datasource;};
78 OGRSFDriver* getDriver(void) const {return m_datasource->GetDriver();};
79#else
80 GDALDataset* getDataSource(void) {return m_datasource;};
81 GDALDriver* getDriver(void) const {return m_datasource->GetDriver();};
82#endif
83 int getLayerCount(void) const {return m_datasource->GetLayerCount();};
84// OGRLayer *executeSql(const std::string& output,const std::string& sqlStatement, OGRGeometry* spatialFilter=NULL);
85 template<typename T> int readSql(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& sqlStatement, OGRGeometry* spatialFilter=NULL, int layer=0, bool pos=false, bool verbose=false);
86 template<typename T> int readSql(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, const std::string& sqlStatement, OGRGeometry* spatialFilter, int layer=0, bool pos=false, bool verbose=false);
87 bool getExtent(double& ulx, double& uly, double& lrx, double& lry, int layer);
88 bool getExtent(double& ulx, double& uly, double& lrx, double& lry);
89
90 void setFieldSeparator(const char fs){ m_fs=fs;};
91 char getFieldSeparator() const { return m_fs;};
92 friend std::ostream& operator<<(std::ostream& theOstream, ImgReaderOgr& theImageReader);
93
94protected:
95 void setCodec(void);
96
97 std::string m_filename;
98#if GDAL_VERSION_MAJOR < 2
99 OGRDataSource *m_datasource;
100#else
101 GDALDataset *m_datasource;
102#endif
103 char m_fs;
104};
105
106//read data from all features in a map, organized by classes
107template <typename T> int ImgReaderOgr::readData(std::map<int,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer, bool pos, bool verbose)
108{
109 assert(m_datasource->GetLayerCount()>layer);
110 OGRLayer *poLayer;
111 if(verbose)
112 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
113 poLayer = m_datasource->GetLayer(layer);
114 if(poLayer!=NULL){
115 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
116 if(fields.empty()){
117 fields.resize(poFDefn->GetFieldCount());
118 if(verbose)
119 std::cout << "resized fields to " << fields.size() << std::endl;
120 }
121 //start reading features from the layer
122 OGRFeature *poFeature;
123 if(verbose)
124 std::cout << "reset reading" << std::endl;
125 poLayer->ResetReading();
126 unsigned long int ifeature=0;
127 int posOffset=(pos)?2:0;
128 if(verbose)
129 std::cout << "going through features" << std::endl << std::flush;
130 int theClass=0;
131 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
132 std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
133 if(verbose)
134 std::cout << "reading feature " << ifeature << std::endl << std::flush;
135 OGRGeometry *poGeometry;
136 poGeometry = poFeature->GetGeometryRef();
137 if(verbose){
138 if(poGeometry == NULL)
139 std::cerr << "no geometry defined" << std::endl << std::flush;
140 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
141 std::cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
142 }
143 assert(poGeometry != NULL );
144 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
145 if(pos){
146 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint){
147 OGRPoint *poPoint;
148 poPoint = (OGRPoint *) poGeometry;
149 theFeature.push_back(poPoint->getX());
150 theFeature.push_back(poPoint->getY());
151 }
152 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
153 OGRPoint thePoint;
154 OGRPolygon * poPolygon = (OGRPolygon *) poGeometry;
155 poPolygon->Centroid(&thePoint);
156 theFeature.push_back(thePoint.getX());
157 theFeature.push_back(thePoint.getY());
158 }
159 else{
160 //Centroid for non polygon geometry not supported until OGR 1.8.0, comment out if version < 1.8.0 is installed...";
161 OGRPoint thePoint;
162 poGeometry->Centroid(&thePoint);
163 theFeature.push_back(thePoint.getX());
164 theFeature.push_back(thePoint.getY());
165 }
166 }
167 // OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
168 std::string featurename;
169 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
170 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
171 std::string fieldname=poFieldDefn->GetNameRef();
172 if(fieldname==label)
173 theClass=poFeature->GetFieldAsInteger(iField);
174 else{
175 switch(fieldType){
176 case(OFTReal):
177 if(fields.size()<poFDefn->GetFieldCount()){
178 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
179 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
180 }
181 else{
182 fields[iField]=fieldname;
183 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
184 }
185 break;
186 case(OFTInteger):
187 if(fields.size()<poFDefn->GetFieldCount()){
188 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
189 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
190 }
191 else{
192 fields[iField]=fieldname;
193 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
194 }
195 break;
196 default:
197 {
198 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
199 throw(errorstring);
200 }
201 break;
202 }
203 }
204 }
205 data[theClass].push_back(theFeature);
206 ++ifeature;
207 ++ifeature;
208 }
209 if(verbose)
210 std::cout << "number of features read: " << ifeature << std::endl << std::flush;
211 typename std::map<int,Vector2d<T> >::const_iterator mit=data.begin();
212 int nband=0;
213 if(verbose)
214 std::cout << "read classes: " << std::flush;
215 while(mit!=data.end()){
216 if(verbose)
217 std::cout << mit->first << " " << std::flush;
218 if(!nband)
219 nband=fields.size();
220 if(pos)
221 assert((mit->second)[0].size()==nband+2);
222 else
223 assert((mit->second)[0].size()==nband);
224 ++mit;
225 }
226 if(verbose)
227 std::cout << std::endl << std::flush;
228 return(nband);
229 }
230 else{
231 std::ostringstream ess;
232 ess << "no layer in " << m_filename;
233 throw(ess.str());
234 }
235}
236
237//read data from all features in a map, organized by class names
238template <typename T> int ImgReaderOgr::readData(std::map<std::string,Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, int layer, bool pos, bool verbose)
239{
240 assert(m_datasource->GetLayerCount()>layer);
241 OGRLayer *poLayer;
242 if(verbose)
243 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
244 poLayer = m_datasource->GetLayer(layer);
245 if(poLayer!=NULL){
246 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
247 assert(poFDefn!=NULL);
248
249 if(fields.empty()){
250 fields.resize(poFDefn->GetFieldCount());
251 if(verbose)
252 std::cout << "resized fields to " << fields.size() << std::endl;
253 }
254
255 //start reading features from the layer
256 OGRFeature *poFeature;
257 if(verbose)
258 std::cout << "reset reading" << std::endl;
259 poLayer->ResetReading();
260 unsigned long int ifeature=0;
261 int posOffset=(pos)?2:0;
262 if(verbose)
263 std::cout << "going through features to fill in string map" << std::endl << std::flush;
264 std::string theClass;
265 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
266 std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
267 if(verbose)
268 std::cout << "reading feature " << ifeature << std::endl << std::flush;
269 OGRGeometry *poGeometry;
270 poGeometry = poFeature->GetGeometryRef();
271 if(verbose){
272 if(poGeometry == NULL)
273 std::cerr << "no geometry defined" << std::endl << std::flush;
274 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
275 std::cerr << "Warning: poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
276 }
277 assert(poGeometry != NULL );
278 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
279 if(pos){
280 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint){
281 OGRPoint *poPoint;
282 poPoint = (OGRPoint *) poGeometry;
283 theFeature.push_back(poPoint->getX());
284 theFeature.push_back(poPoint->getY());
285 }
286 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
287 OGRPoint thePoint;
288 poGeometry->Centroid(&thePoint);
289 theFeature.push_back(thePoint.getX());
290 theFeature.push_back(thePoint.getY());
291 }
292 else{
293 //Centroid for non polygon geometry not supported until OGR 1.8.0, comment out if version < 1.8.0 is installed...";
294 OGRPoint thePoint;
295 poGeometry->Centroid(&thePoint);
296 theFeature.push_back(thePoint.getX());
297 theFeature.push_back(thePoint.getY());
298 }
299 }
300 // OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();//got LayerDefn already...
301 std::string featurename;
302 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
303 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
304 std::string fieldname=poFieldDefn->GetNameRef();
305 if(fieldname==label){
306 theClass=poFeature->GetFieldAsString(iField);
307 if(verbose)
308 std::cout << "read feature for " << theClass << std::endl;
309 }
310 else{
311 switch(fieldType){
312 case(OFTReal):
313 if(fields.size()<poFDefn->GetFieldCount()){
314 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
315 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
316 }
317 else{
318 fields[iField]=fieldname;
319 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
320 }
321 break;
322 case(OFTInteger):
323 if(fields.size()<poFDefn->GetFieldCount()){
324 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
325 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
326 }
327 else{
328 fields[iField]=fieldname;
329 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
330 }
331 break;
332 default:
333 {
334 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
335 throw(errorstring);
336 }
337 break;
338 }
339 }
340 assert(poFDefn!=NULL);
341 }
342 data[theClass].push_back(theFeature);
343 ++ifeature;
344 }
345 if(verbose)
346 std::cout << "number of features read: " << ifeature << std::endl << std::flush;
347 typename std::map<std::string,Vector2d<T> >::const_iterator mit=data.begin();
348 int nband=0;
349 if(verbose)
350 std::cout << "read classes: " << std::flush;
351 while(mit!=data.end()){
352 if(verbose)
353 std::cout << mit->first << " " << std::flush;
354 if(!nband)
355 nband=fields.size();
356 if(pos)
357 assert((mit->second)[0].size()==nband+2);
358 else
359 assert((mit->second)[0].size()==nband);
360 ++mit;
361 }
362 if(verbose)
363 std::cout << std::endl << std::flush;
364 return(nband);
365 }
366 else{
367 std::ostringstream ess;
368 ess << "no layer in " << m_filename;
369 throw(ess.str());
370 }
371}
372
373//read x positions
374template <typename T> int ImgReaderOgr::readXY(std::vector<T>& xVector, std::vector<T>& yVector, int layer, bool verbose){
375 assert(m_datasource->GetLayerCount()>layer);
376 OGRLayer *poLayer;
377 if(verbose)
378 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
379 poLayer = m_datasource->GetLayer(layer);
380 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
381 //start reading features from the layer
382 OGRFeature *poFeature;
383 if(verbose)
384 std::cout << "reset reading" << std::endl;
385 poLayer->ResetReading();
386 unsigned long int ifeature=0;
387 if(verbose)
388 std::cout << "going through features" << std::endl << std::flush;
389 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
390 if(verbose)
391 std::cout << "reading feature " << ifeature << std::endl << std::flush;
392 OGRGeometry *poGeometry;
393 poGeometry = poFeature->GetGeometryRef();
394 if(verbose){
395 if(poGeometry == NULL)
396 std::cerr << "no geometry defined" << std::endl << std::flush;
397 else// if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
398 std::cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl;
399 }
400 // assert(poGeometry != NULL
401 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
402 OGRPoint *poPoint = (OGRPoint *) poGeometry;
403 xVector.push_back(poPoint->getX());
404 yVector.push_back(poPoint->getY());
405 ++ifeature;
406 }
407 assert(xVector.size()==yVector.size());
408 if(xVector.size()){
409 return xVector.size();
410 }
411 else{
412 std::ostringstream ess;
413 ess << "no layer in " << m_filename;
414 throw(ess.str());
415 }
416}
417
418//read data from a single feature
419template <typename T> int ImgReaderOgr::readData(std::vector<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, OGRFeature *poFeature, int layer, bool pos, bool verbose)
420{
421 assert(m_datasource->GetLayerCount()>layer);
422 OGRLayer *poLayer;
423 if(verbose)
424 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
425 poLayer = m_datasource->GetLayer(layer);
426 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
427 if(fields.empty()){
428 fields.resize(poFDefn->GetFieldCount());
429 if(verbose)
430 std::cout << "resized fields to " << fields.size() << std::endl;
431 }
432 OGRGeometry *poGeometry;
433 poGeometry = poFeature->GetGeometryRef();
434 if(verbose){
435 if(poGeometry == NULL)
436 std::cerr << "no geometry defined" << std::endl << std::flush;
437 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
438 std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
439 }
440 assert(poGeometry != NULL);
441 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
442 OGRPoint *poPoint = (OGRPoint *) poGeometry;
443 if(pos){
444 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint){
445 OGRPoint *poPoint;
446 poPoint = (OGRPoint *) poGeometry;
447 data.push_back(poPoint->getX());
448 data.push_back(poPoint->getY());
449 }
450 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
451 OGRPoint thePoint;
452 poGeometry->Centroid(&thePoint);
453 data.push_back(thePoint.getX());
454 data.push_back(thePoint.getY());
455 }
456 else{
457 //Centroid for non polygon geometry not supported until OGR 1.8.0, comment out if version < 1.8.0 is installed...";
458 OGRPoint thePoint;
459 poGeometry->Centroid(&thePoint);
460 data.push_back(thePoint.getX());
461 data.push_back(thePoint.getY());
462 }
463
464 }
465 std::string featurename;
466 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
467 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
468 std::string fieldname=poFieldDefn->GetNameRef();
469 switch(fieldType){
470 case(OFTReal):
471 if(fields.size()<poFDefn->GetFieldCount()){
472 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
473 data.push_back(poFeature->GetFieldAsDouble(iField));
474 }
475 else{
476 fields[iField]=fieldname;
477 data.push_back(poFeature->GetFieldAsDouble(iField));
478 }
479 break;
480 case(OFTInteger):
481 if(fields.size()<poFDefn->GetFieldCount()){
482 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
483 data.push_back(poFeature->GetFieldAsDouble(iField));
484 }
485 else{
486 fields[iField]=fieldname;
487 data.push_back(poFeature->GetFieldAsDouble(iField));
488 }
489 break;
490 default:
491 {
492 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
493 throw(errorstring);
494 }
495 break;
496 }
497 }
498 // assert(data.size()==ifeature);
499 if(data.size()){
500 if(pos)
501 assert(data.size()==fields.size()+2);
502 else
503 assert(data.size()==fields.size());
504 return fields.size();
505 }
506 else{
507 std::ostringstream ess;
508 ess << "no layer in " << m_filename;
509 throw(ess.str());
510 }
511}
512
513//read one field from all features
514template <typename T> inline int ImgReaderOgr::readData(std::vector<T>& data, const OGRFieldType& fieldType, const std::string& theField, int layer, bool verbose)
515{
516 assert(m_datasource->GetLayerCount()>layer);
517 OGRLayer *poLayer;
518 if(verbose)
519 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
520 poLayer = m_datasource->GetLayer(layer);
521 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
522 int nfield=(theField!="")? poFDefn->GetFieldCount() : 1;
523 if(theField==""){
524 //read first field available
525 if(verbose)
526 std::cout << "read first field from total of " << nfield << std::endl;
527 }
528
529 //start reading features from the layer
530 OGRFeature *poFeature;
531 if(verbose)
532 std::cout << "reset reading" << std::endl;
533 poLayer->ResetReading();
534 unsigned long int ifeature=0;
535 if(verbose)
536 std::cout << "going through features" << std::endl << std::flush;
537 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
538 // std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
539 T theFeature;
540 if(verbose)
541 std::cout << "reading feature " << ifeature << std::endl << std::flush;
542 OGRGeometry *poGeometry;
543 poGeometry = poFeature->GetGeometryRef();
544 if(verbose){
545 if(poGeometry == NULL)
546 std::cerr << "no geometry defined" << std::endl << std::flush;
547 else// if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
548 std::cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl;
549 }
550 // assert(poGeometry != NULL
551 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
552 OGRPoint *poPoint = (OGRPoint *) poGeometry;
553
554 for(int iField=0;iField<nfield;++iField){
555 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
556 std::string fieldname=poFieldDefn->GetNameRef();
557 if(fieldname!=theField)
558 continue;
559 switch(fieldType){
560 case(OFTInteger):
561 case(OFTReal):
562 theFeature=poFeature->GetFieldAsDouble(iField);
563 break;
564 default:
565 {
566 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
567 throw(errorstring);
568 }
569 break;
570 }
571 }
572 data.push_back(theFeature);
573 if(verbose)
574 std::cout << "feature is: " << theFeature << std::endl;
575 ++ifeature;
576 }
577 if(data.size()){
578 return data.size();
579 }
580 else{
581 std::ostringstream ess;
582 ess << "no layer in " << m_filename;
583 throw(ess.str());
584 }
585}
586
587//specialization for string: read one field from all features
588template <> inline int ImgReaderOgr::readData(std::vector<std::string>& data, const OGRFieldType& fieldType, const std::string& theField, int layer, bool verbose)
589{
590 assert(m_datasource->GetLayerCount()>layer);
591 OGRLayer *poLayer;
592 if(verbose)
593 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
594 poLayer = m_datasource->GetLayer(layer);
595 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
596 int nfield=(theField!="")? poFDefn->GetFieldCount() : 1;
597 if(theField==""){
598 //read first field available
599 if(verbose)
600 std::cout << "read first field from total of " << nfield << std::endl;
601 }
602
603 //start reading features from the layer
604 OGRFeature *poFeature;
605 if(verbose)
606 std::cout << "reset reading" << std::endl;
607 poLayer->ResetReading();
608 unsigned long int ifeature=0;
609 if(verbose)
610 std::cout << "going through features" << std::endl << std::flush;
611 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
612 std::string theFeature;
613 if(verbose)
614 std::cout << "reading feature " << ifeature << std::endl << std::flush;
615 OGRGeometry *poGeometry;
616 poGeometry = poFeature->GetGeometryRef();
617 if(verbose){
618 if(poGeometry == NULL)
619 std::cerr << "no geometry defined" << std::endl << std::flush;
620 else// if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
621 std::cout << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl;
622 }
623 // assert(poGeometry != NULL
624 // && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
625 OGRPoint *poPoint = (OGRPoint *) poGeometry;
626
627 for(int iField=0;iField<nfield;++iField){
628 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
629 std::string fieldname=poFieldDefn->GetNameRef();
630 if(fieldname!=theField)
631 continue;
632 switch(fieldType){
633 case(OFTInteger):
634 case(OFTReal):
635 case(OFTString):
636 theFeature=poFeature->GetFieldAsString(iField);
637 break;
638 default:
639 {
640 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
641 throw(errorstring);
642 }
643 break;
644 }
645 }
646 data.push_back(theFeature);
647 if(verbose)
648 std::cout << "feature is: " << theFeature << std::endl;
649 ++ifeature;
650 }
651 if(data.size()){
652 return data.size();
653 }
654 else{
655 std::ostringstream ess;
656 ess << "no layer in " << m_filename;
657 throw(ess.str());
658 }
659}
660
661//read data from all features
662template <typename T> int ImgReaderOgr::readData(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, int layer, bool pos, bool verbose)
663{
664 assert(m_datasource->GetLayerCount()>layer);
665 OGRLayer *poLayer;
666 if(verbose)
667 std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
668 poLayer = m_datasource->GetLayer(layer);
669 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
670 if(fields.empty()){
671 fields.resize(poFDefn->GetFieldCount());
672 if(verbose)
673 std::cout << "resized fields to " << fields.size() << std::endl;
674 }
675 //start reading features from the layer
676 OGRFeature *poFeature;
677 if(verbose)
678 std::cout << "reset reading" << std::endl;
679 poLayer->ResetReading();
680 unsigned long int ifeature=0;
681 int posOffset=(pos)?2:0;
682 if(verbose)
683 std::cout << "going through features" << std::endl << std::flush;
684 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
685 std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
686 if(verbose)
687 std::cout << "reading feature " << ifeature << std::endl << std::flush;
688 OGRGeometry *poGeometry;
689 poGeometry = poFeature->GetGeometryRef();
690 if(verbose){
691 if(poGeometry == NULL)
692 std::cerr << "no geometry defined" << std::endl << std::flush;
693 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
694 std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
695 }
696 assert(poGeometry != NULL
697 && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
698 OGRPoint *poPoint = (OGRPoint *) poGeometry;
699 if(pos){
700 theFeature.push_back(poPoint->getX());
701 theFeature.push_back(poPoint->getY());
702 }
703 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
704 std::string featurename;
705 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
706 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
707 std::string fieldname=poFieldDefn->GetNameRef();
708 switch(fieldType){
709 case(OFTReal):
710 if(fields.size()<poFDefn->GetFieldCount()){
711 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
712 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
713 }
714 else{
715 fields[iField]=fieldname;
716 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
717 }
718 break;
719 case(OFTInteger):
720 if(fields.size()<poFDefn->GetFieldCount()){
721 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
722 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
723 }
724 else{
725 fields[iField]=fieldname;
726 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
727 }
728 break;
729 default:
730 {
731 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
732 throw(errorstring);
733 }
734 break;
735 }
736 }
737 data.push_back(theFeature);
738 ++ifeature;
739 }
740// assert(data.size()==ifeature);
741 if(data.size()){
742 if(pos)
743 assert(data[0].size()==fields.size()+2);
744 else
745 assert(data[0].size()==fields.size());
746 return fields.size();
747 }
748 else{
749 std::ostringstream ess;
750 ess << "no layer in " << m_filename;
751 throw(ess.str());
752 }
753}
754
755template<typename T> int ImgReaderOgr::readSql(std::map<int, Vector2d<T> >& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& label, const std::string& sqlStatement, OGRGeometry* spatialFilter, int layer, bool pos, bool verbose)
756{
757 assert(m_datasource->GetLayerCount()>layer);
758 OGRLayer *poLayer;
759 poLayer = m_datasource->ExecuteSQL(sqlStatement.c_str(), spatialFilter,NULL );
760 if(poLayer!=NULL){
761 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
762 if(fields.empty()){
763 fields.resize(poFDefn->GetFieldCount());
764 if(verbose)
765 std::cout << "resized fields to " << fields.size() << std::endl;
766 }
767 //start reading features from the layer
768 OGRFeature *poFeature;
769 if(verbose)
770 std::cout << "reset reading" << std::endl;
771 poLayer->ResetReading();
772 unsigned long int ifeature=0;
773 int posOffset=(pos)?2:0;
774 if(verbose)
775 std::cout << "going through features" << std::endl << std::flush;
776 int theClass=0;
777 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
778 std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
779 if(verbose)
780 std::cout << "reading feature " << ifeature << std::endl << std::flush;
781 OGRGeometry *poGeometry;
782 poGeometry = poFeature->GetGeometryRef();
783 if(verbose){
784 if(poGeometry == NULL)
785 std::cerr << "no geometry defined" << std::endl << std::flush;
786 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
787 std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
788 }
789 assert(poGeometry != NULL
790 && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
791 OGRPoint *poPoint = (OGRPoint *) poGeometry;
792 if(pos){
793 theFeature.push_back(poPoint->getX());
794 theFeature.push_back(poPoint->getY());
795 }
796 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
797 std::string featurename;
798 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
799 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
800 std::string fieldname=poFieldDefn->GetNameRef();
801 if(fieldname==label)
802 theClass=poFeature->GetFieldAsInteger(iField);
803 else{
804 switch(fieldType){
805 case(OFTReal):
806 if(fields.size()<poFDefn->GetFieldCount()){
807 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
808 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
809 }
810 else{
811 fields[iField]=fieldname;
812 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
813 }
814 break;
815 case(OFTInteger):
816 if(fields.size()<poFDefn->GetFieldCount()){
817 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
818 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
819 }
820 else{
821 fields[iField]=fieldname;
822 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
823 }
824 break;
825 default:
826 {
827 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
828 throw(errorstring);
829 }
830 break;
831 }
832 }
833 }
834 data[theClass].push_back(theFeature);
835 ++ifeature;
836 }
837 if(verbose)
838 std::cout << "number of features read: " << ifeature << std::endl << std::flush;
839 typename std::map<int,Vector2d<T> >::const_iterator mit=data.begin();
840 int nband=0;
841 if(verbose)
842 std::cout << "read classes: " << std::flush;
843 while(mit!=data.end()){
844 if(verbose)
845 std::cout << mit->first << " " << std::flush;
846 if(!nband)
847 nband=fields.size();
848 if(pos)
849 assert((mit->second)[0].size()==nband+2);
850 else
851 assert((mit->second)[0].size()==nband);
852 ++mit;
853 }
854 if(verbose)
855 std::cout << std::endl << std::flush;
856 return(nband);
857 }
858 else{
859 std::ostringstream ess;
860 ess << "no layer in " << m_filename;
861 throw(ess.str());
862 }
863}
864
865template<typename T> int ImgReaderOgr::readSql(Vector2d<T>& data, const OGRFieldType& fieldType, std::vector<std::string>& fields, const std::string& sqlStatement, OGRGeometry* spatialFilter, int layer, bool pos, bool verbose)
866{
867 assert(m_datasource->GetLayerCount()>layer);
868 OGRLayer *poLayer;
869 poLayer = m_datasource->ExecuteSQL(sqlStatement.c_str(), spatialFilter,NULL );
870 if(poLayer!=NULL){
871 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
872 if(fields.empty()){
873 fields.resize(poFDefn->GetFieldCount());
874 if(verbose)
875 std::cout << "resized fields to " << fields.size() << std::endl;
876 }
877 //start reading features from the layer
878 OGRFeature *poFeature;
879 if(verbose)
880 std::cout << "reset reading" << std::endl;
881 poLayer->ResetReading();
882 unsigned long int ifeature=0;
883 int posOffset=(pos)?2:0;
884 if(verbose)
885 std::cout << "going through features" << std::endl << std::flush;
886 while( (poFeature = poLayer->GetNextFeature()) != NULL ){
887 std::vector<T> theFeature;//(fields.size()+posOffset);//x,y+selectedfields
888 if(verbose)
889 std::cout << "reading feature " << ifeature << std::endl << std::flush;
890 OGRGeometry *poGeometry;
891 poGeometry = poFeature->GetGeometryRef();
892 if(verbose){
893 if(poGeometry == NULL)
894 std::cerr << "no geometry defined" << std::endl << std::flush;
895 else if(wkbFlatten(poGeometry->getGeometryType()) != wkbPoint)
896 std::cerr << "poGeometry type: " << wkbFlatten(poGeometry->getGeometryType()) << std::endl << std::flush;
897 }
898 assert(poGeometry != NULL
899 && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
900 OGRPoint *poPoint = (OGRPoint *) poGeometry;
901 if(pos){
902 theFeature.push_back(poPoint->getX());
903 theFeature.push_back(poPoint->getY());
904 }
905 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
906 std::string featurename;
907 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
908 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
909 std::string fieldname=poFieldDefn->GetNameRef();
910 switch(fieldType){
911 case(OFTReal):
912 if(fields.size()<poFDefn->GetFieldCount()){
913 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
914 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
915 }
916 else{
917 fields[iField]=fieldname;
918 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
919 }
920 break;
921 case(OFTInteger):
922 if(fields.size()<poFDefn->GetFieldCount()){
923 if(find(fields.begin(),fields.end(),fieldname)!=fields.end())
924 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
925 }
926 else{
927 fields[iField]=fieldname;
928 theFeature.push_back(poFeature->GetFieldAsDouble(iField));
929 }
930 break;
931 default:
932 {
933 std::string errorstring="field type not supported in ImgReaderOgr::ReadData";
934 throw(errorstring);
935 }
936 break;
937 }
938 }
939 data.push_back(theFeature);
940 ++ifeature;
941 }
942 m_datasource->ReleaseResultSet( poLayer );
943 // assert(data.size()==ifeature);
944 if(data.size()){
945 if(pos)
946 assert(data[0].size()==fields.size()+2);
947 else
948 assert(data[0].size()==fields.size());
949 return fields.size();
950 }
951 else
952 return(0);
953 }
954 else{
955 std::ostringstream ess;
956 ess << "no layer in " << m_filename;
957 throw(ess.str());
958 }
959}
960
961#endif // _IMGREADEROGR_H_