pktools 2.6.7
Processing Kernel for geospatial data
ImgReaderOgr.cc
1/**********************************************************************
2ImgReaderOgr.cc: 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#include <iostream>
21#include <fstream>
22#include "ImgReaderOgr.h"
23#include "ImgWriterOgr.h"
24#include "cpl_string.h"
25//---------------------------------------------------------------------------
26ImgReaderOgr::ImgReaderOgr(void)
27 : m_fs(' ')
28{}
29
30ImgReaderOgr::ImgReaderOgr(const std::string& filename)
31{
32 open(filename);
33}
34
35ImgReaderOgr::~ImgReaderOgr(void)
36{
37}
38
39//---------------------------------------------------------------------------
40
41void ImgReaderOgr::open(const std::string& filename)
42{
43 m_fs=' ';
44 m_filename = filename;
45 setCodec();
46}
47
48//---------------------------------------------------------------------------
49void ImgReaderOgr::close(void)
50{
51#if GDAL_VERSION_MAJOR < 2
52 OGRDataSource::DestroyDataSource(m_datasource);
53#else
54 GDALClose(m_datasource);
55#endif
56}
57
58//---------------------------------------------------------------------------
59void ImgReaderOgr::setCodec(void){
60#if GDAL_VERSION_MAJOR < 2
61 //register the drivers
62 OGRRegisterAll();
63 //open the input OGR datasource. Datasources can be files, RDBMSes, directories full of files, or even remote web services depending on the driver being used. However, the datasource name is always a single string.
64 m_datasource = OGRSFDriverRegistrar::Open(m_filename.c_str(), FALSE);//FAlSE: do not update
65#else
66 //register the drivers
67 GDALAllRegister();
68 //open the input OGR datasource. Datasources can be files, RDBMSes, directories full of files, or even remote web services depending on the driver being used. However, the datasource name is always a single string.
69 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_VECTOR||GDAL_OF_READONLY, NULL, NULL, NULL);
70 // m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_READONLY||GDAL_OF_VECTOR, NULL, NULL, NULL);
71#endif
72 if( m_datasource == NULL ){
73#if GDAL_VERSION_MAJOR < 2
74 std::string errorString="Open failed";
75 throw(errorString);
76#else
77 m_datasource = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
78 if( m_datasource == NULL ){
79 std::string errorString="Open failed";
80 throw(errorString);
81 }
82#endif
83 }
84}
85
86bool ImgReaderOgr::getExtent(double& ulx, double& uly, double& lrx, double& lry, int layer)
87{
88 OGREnvelope oExt;
89 if(getLayer(layer)->GetExtent(&oExt,TRUE)==OGRERR_NONE){
90 ulx=oExt.MinX;
91 uly=oExt.MaxY;
92 lrx=oExt.MaxX;
93 lry=oExt.MinY;
94 return true;
95 }
96 else
97 return false;
98}
99
100bool ImgReaderOgr::getExtent(double& ulx, double& uly, double& lrx, double& lry)
101{
102 bool success=false;
103 OGREnvelope oExt;
104 for(int ilayer=0;ilayer<getLayerCount();++ilayer){
105 if(getLayer(ilayer)->GetExtent(&oExt,TRUE)==OGRERR_NONE){
106 if(!ilayer){
107 ulx=oExt.MinX;
108 uly=oExt.MaxY;
109 lrx=oExt.MaxX;
110 lry=oExt.MinY;
111 }
112 else{
113 if(ulx>oExt.MinX)
114 ulx=oExt.MinX;
115 if(uly<oExt.MaxY)
116 uly=oExt.MaxY;
117 if(lrx<oExt.MaxX)
118 lrx=oExt.MaxX;
119 if(lry>oExt.MinY)
120 lry=oExt.MinY;
121 }
122 success=true;
123 }
124 else
125 success=false;
126 }
127 return success;
128}
129
130unsigned long int ImgReaderOgr::getFeatureCount(int layer) const
131{
132 return(m_datasource->GetLayer(layer)->GetFeatureCount());
133}
134
135int ImgReaderOgr::getFieldCount(int layer) const
136{
137 if(layer<0)
138 layer=m_datasource->GetLayerCount()-1;
139 assert(m_datasource->GetLayerCount()>layer);
140 OGRLayer *poLayer;
141 if((poLayer = m_datasource->GetLayer(layer))==NULL){
142 std::string errorstring="Could not get layer";
143 throw(errorstring);
144 }
145 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
146 return(poFDefn->GetFieldCount());
147}
148
149int ImgReaderOgr::getFields(std::vector<std::string>& fields, int layer) const
150{
151 if(layer<0)
152 layer=m_datasource->GetLayerCount()-1;
153 assert(m_datasource->GetLayerCount()>layer);
154 OGRLayer *poLayer;
155 if((poLayer = m_datasource->GetLayer(layer))==NULL){
156 std::string errorstring="Could not get layer";
157 throw(errorstring);
158 }
159 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
160 fields.clear();
161 fields.resize(poFDefn->GetFieldCount());
162 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
163 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
164 fields[iField]=poFieldDefn->GetNameRef();
165 }
166 return(fields.size());
167}
168
169int ImgReaderOgr::getFields(std::vector<OGRFieldDefn*>& fields, int layer) const
170{
171 if(layer<0)
172 layer=m_datasource->GetLayerCount()-1;
173 assert(m_datasource->GetLayerCount()>layer);
174 OGRLayer *poLayer;
175 if((poLayer = m_datasource->GetLayer(layer))==NULL){
176 std::string errorstring="Could not get layer";
177 throw(errorstring);
178 }
179 OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
180 fields.clear();
181 fields.resize(poFDefn->GetFieldCount());
182 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
183 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
184 fields[iField]=poFDefn->GetFieldDefn(iField);
185 }
186 assert(fields.size()==getFieldCount(layer));
187 return(fields.size());
188}
189
190std::string ImgReaderOgr::getProjection(int layer) const
191{
192 if(m_datasource->GetLayer(layer)->GetSpatialRef()){
193 char* ppszResult;
194 m_datasource->GetLayer(layer)->GetSpatialRef()->exportToWkt(&ppszResult);
195 return(ppszResult);
196 }
197 else
198 return "";
199}
200
201OGRwkbGeometryType ImgReaderOgr::getGeometryType(int layer) const
202{
203 return m_datasource->GetLayer(layer)->GetLayerDefn()->GetGeomType();
204}
205
206std::ostream& operator<<(std::ostream& theOstream, ImgReaderOgr& theImageReader){
207 //An OGRDataSource can potentially have many layers associated with it. The number of layers available can be queried with OGRDataSource::GetLayerCount() and individual layers fetched by index using OGRDataSource::GetLayer(). However, we wil just fetch the layer by name.
208 //todo: try to open and catch if failure...
209 // ofstream fpoints(filename.c_str(),ios::out);
210
211 int nlayerRead=theImageReader.getDataSource()->GetLayerCount();
212
213 for(int ilayer=0;ilayer<nlayerRead;++ilayer){
214 OGRLayer *readLayer=theImageReader.getLayer(ilayer);
215 OGRFeatureDefn *poFDefn = readLayer->GetLayerDefn();
216
217 theOstream << "#";
218 int iField=0;
219 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
220 OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
221 std::string fieldname=poFieldDefn->GetNameRef();
222 theOstream << fieldname << theImageReader.getFieldSeparator();
223 }
224 theOstream << std::endl;
225
226 readLayer->ResetReading();
227
228 //start reading features from the layer
229 OGRFeature *poFeature;
230 unsigned long int ifeature=0;
231 while( (poFeature = readLayer->GetNextFeature()) != NULL ){
232 OGRGeometry *poGeometry;
233 poGeometry = poFeature->GetGeometryRef();
234 assert(poGeometry != NULL);
235 double x,y;
236 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint){
237 OGRPoint *poPoint = (OGRPoint *) poGeometry;
238 x=poPoint->getX();
239 y=poPoint->getY();
240 }
241 std::vector<std::string> vfields(poFDefn->GetFieldCount());
242 std::string featurename;
243 std::vector<std::string>::iterator fit=vfields.begin();
244 for(int iField=0;iField<poFDefn->GetFieldCount();++iField){
245 *(fit++)=poFeature->GetFieldAsString(iField);
246 }
247 theOstream.precision(12);
248 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint)
249 theOstream << x << theImageReader.getFieldSeparator() << y;
250 for(fit=vfields.begin();fit!=vfields.end();++fit)
251 theOstream << theImageReader.getFieldSeparator() << *fit;
252 theOstream << std::endl;
253 ++ifeature;
254 }
255 }
256 return(theOstream);
257}
258
259// OGRLayer * ImgReaderOgr::executeSql(const std::string& output, const std::string& sqlStatement, OGRGeometry* spatialFilter)
260// {
261// OGRLayer *poResultSet;
262// poResultSet = m_datasource->ExecuteSQL(sqlStatement.c_str(), spatialFilter,NULL );
263
264// if( poResultSet != NULL ){
265// ImgWriterOgr imgWriter;
266// imgWriter.open(output);
267// imgWriter.copyLayer(poResultSet,output);
268// m_datasource->ReleaseResultSet( poResultSet );
269// imgWriter.close();
270// }
271// }
272
273unsigned int ImgReaderOgr::readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
274 std::vector<std::string>& fields,
275 const std::vector<unsigned short>& bands,
276 const std::string& label,
277 const std::vector<std::string>& layers,
278 int verbose)
279{
280 mapPixels.clear();
281 int nsample=0;
282 int totalSamples=0;
283 int nband=0;
284 if(verbose)
285 std::cout << "reading OGR dataset " << m_filename << std::endl;
286 for(int ilayer=0;ilayer<getLayerCount();++ilayer){
287 std::string currentLayername=getLayer(ilayer)->GetName();
288 if(layers.size())
289 if(find(layers.begin(),layers.end(),currentLayername)==layers.end())
290 continue;
291 try{
292 //only retain bands in fields
293 getFields(fields,ilayer);
294 std::vector<std::string>::iterator fit=fields.begin();
295 if(verbose>1)
296 std::cout << "reading fields: ";
297 while(fit!=fields.end()){
298 if(verbose)
299 std::cout << *fit << " ";
300 // size_t pos=(*fit).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ");
301 if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
302 // if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
303 std::size_t digits=(*fit).substr(1,1).find_first_of("0123456789");
304 std::size_t digite=(*fit).substr(1).find_first_not_of("0123456789");
305 if((*fit)=="B" || (*fit)=="b" || (*fit)=="Band")//B is only band
306 ++fit;
307 else if(digits!=std::string::npos&&digite==std::string::npos){
308 std::string digitString=(*fit).substr(digits);
309 // int theBand=atoi((*fit).substr(1).c_str());
310 int theBand=atoi(digitString.c_str());
311 if(bands.size()){
312 bool validBand=false;
313 for(int iband=0;iband<bands.size();++iband){
314 if(theBand==bands[iband])
315 validBand=true;
316 }
317 if(validBand)
318 ++fit;
319 else
320 fields.erase(fit);
321 }
322 else
323 ++fit;
324 }
325 else
326 fields.erase(fit);
327 }
328 else
329 fields.erase(fit);
330 }
331 if(verbose)
332 std::cout << std::endl;
333 if(verbose){
334 std::cout << "fields:";
335 for(std::vector<std::string>::iterator fit=fields.begin();fit!=fields.end();++fit)
336 std::cout << " " << *fit;
337 std::cout << std::endl;
338 }
339 if(!nband){
340 if(verbose)
341 std::cout << "reading data" << std::endl;
342 nband=readData(mapPixels,OFTReal,fields,label,ilayer,true,verbose==2);
343 }
344 else{
345 assert(nband==readData(mapPixels,OFTReal,fields,label,ilayer,true,false));
346 }
347 nsample=getFeatureCount(ilayer);
348 totalSamples+=nsample;
349 if(verbose)
350 std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
351 }
352 catch(std::string e){
353 std::ostringstream estr;
354 estr << e << " " << m_filename;
355 throw(estr.str());
356 }
357 }
358 if(verbose)
359 std::cout << "total number of samples read " << totalSamples << std::endl;
360 return totalSamples;
361}
362
363unsigned int ImgReaderOgr::readDataImageOgr(std::map<std::string,Vector2d<float> > &mapPixels, //[classNr][pixelNr][bandNr],
364 std::vector<std::string>& fields,
365 double start,
366 double end,
367 const std::string& label,
368 const std::vector<std::string>& layers,
369 int verbose)
370{
371 mapPixels.clear();
372 int nsample=0;
373 int totalSamples=0;
374 int nband=0;
375 if(verbose)
376 std::cout << "reading OGR dataset file " << m_filename << std::endl;
377 for(int ilayer=0;ilayer<getLayerCount();++ilayer){
378 std::string currentLayername=getLayer(ilayer)->GetName();
379 if(layers.size())
380 if(find(layers.begin(),layers.end(),currentLayername)==layers.end())
381 continue;
382 try{
383 //only retain bands in fields
384 getFields(fields,ilayer);
385 std::vector<std::string>::iterator fit=fields.begin();
386 if(verbose)
387 std::cout << "reading fields: ";
388 while(fit!=fields.end()){
389 if(verbose)
390 std::cout << *fit << " ";
391 if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
392 // if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
393 std::size_t digits=(*fit).substr(1,1).find_first_of("0123456789");
394 std::size_t digite=(*fit).substr(1).find_first_not_of("0123456789");
395 if(*fit=="B" || *fit=="b"|| *fit=="Band")
396 ++fit;
397 else if(digits!=std::string::npos&&digite==std::string::npos){
398 std::string digitString=(*fit).substr(digits);
399 int iband=atoi(digitString.c_str());
400 // int iband=atoi((*fit).substr(1).c_str());
401 if((start||end)&&(iband<start||iband>end))
402 fields.erase(fit);
403 else
404 ++fit;
405 }
406 else
407 fields.erase(fit);
408 }
409 else
410 fields.erase(fit);
411 }
412 if(verbose)
413 std::cout << std::endl;
414 if(verbose){
415 std::cout << "fields:";
416 for(std::vector<std::string>::iterator fit=fields.begin();fit!=fields.end();++fit)
417 std::cout << " " << *fit;
418 std::cout << std::endl;
419 }
420 if(!nband){
421 if(verbose)
422 std::cout << "reading data" << std::endl;
423 nband=readData(mapPixels,OFTReal,fields,label,ilayer,true,verbose==2);
424 }
425 else{
426 assert(nband==readData(mapPixels,OFTReal,fields,label,ilayer,true,false));
427 }
428 nsample=getFeatureCount(ilayer);
429 totalSamples+=nsample;
430 if(verbose)
431 std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
432 }
433 catch(std::string e){
434 std::ostringstream estr;
435 estr << e << " " << m_filename;
436 throw(estr.str());
437 }
438 if(verbose)
439 std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
440 }
441 if(verbose)
442 std::cout << "total number of samples read " << totalSamples << std::endl;
443 return totalSamples;
444}