pktools 2.6.7
Processing Kernel for geospatial data
pkextractogr.cc
1/**********************************************************************
2pkextractogr.cc: extract pixel values from raster image from a (vector or raster) sample
3Copyright (C) 2008-2017 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 <assert.h>
21#include <math.h>
22#include <stdlib.h>
23#include <sstream>
24#include <string>
25#include <algorithm>
26#include <ctime>
27#include <vector>
28#include "imageclasses/ImgReaderGdal.h"
29#include "imageclasses/ImgWriterOgr.h"
30#include "base/Optionpk.h"
31#include "algorithms/StatFactory.h"
32
33#ifndef PI
34#define PI 3.1415926535897932384626433832795
35#endif
36
37/******************************************************************************/
116namespace rule{
117 enum RULE_TYPE {point=0, mean=1, proportion=2, custom=3, min=4, max=5, mode=6, centroid=7, sum=8, median=9, stdev=10, percentile=11, count=12, allpoints=13};
118}
119
120using namespace std;
121
122int main(int argc, char *argv[])
123{
124 Optionpk<string> image_opt("i", "input", "Raster input dataset containing band information");
125 Optionpk<string> sample_opt("s", "sample", "OGR vector dataset with features to be extracted from input data. Output will contain features with input band information included. Sample image can also be GDAL raster dataset.");
126 Optionpk<string> layer_opt("ln", "ln", "Layer name(s) in sample (leave empty to select all)");
127 Optionpk<unsigned int> random_opt("rand", "random", "Create simple random sample of points. Provide number of points to generate");
128 Optionpk<double> grid_opt("grid", "grid", "Create systematic grid of points. Provide cell grid size (in projected units, e.g,. m)");
129 Optionpk<string> output_opt("o", "output", "Output sample dataset");
130 Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Leave empty to extract all valid data pixels from sample dataset. Make sure to set classes if rule is set to mode, proportion or count");
131 Optionpk<float> threshold_opt("t", "threshold", "Probability threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold per vector sample layer. If using raster land cover maps as a sample dataset, you can provide a threshold value for each class (e.g. -t 80 -t 60). Use value 100 to select all pixels for selected class(es)", 100);
132 Optionpk<double> percentile_opt("perc","perc","Percentile value used for rule percentile",95);
133 Optionpk<string> ogrformat_opt("f", "f", "Output sample dataset format","SQLite");
134 Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or Integer)", "Real");
135 Optionpk<int> band_opt("b", "band", "Band index(es) to extract (0 based). Leave empty to use all bands");
136 Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
137 Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
138 Optionpk<string> rule_opt("r", "rule", "Rule how to report image information per feature (only for vector sample). point (single point within polygon), allpoints (all points within polygon), centroid, mean, stdev, median, proportion, count, min, max, mode, sum, percentile.","centroid");
139 Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "Invalid value(s) for input image");
140 Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Band(s) in input image to check if pixel is valid (used for srcnodata)", 0);
141 Optionpk<float> polythreshold_opt("tp", "thresholdPolygon", "(absolute) threshold for selecting samples in each polygon");
142 Optionpk<short> buffer_opt("buf", "buffer", "Buffer for calculating statistics for point features (in number of pixels) ",0);
143 Optionpk<bool> disc_opt("circ", "circular", "Use a circular disc kernel buffer (for vector point sample datasets only, use in combination with buffer option)", false);
144 Optionpk<short> verbose_opt("v", "verbose", "Verbose mode if > 0", 0,2);
145
146 bstart_opt.setHide(1);
147 bend_opt.setHide(1);
148 bndnodata_opt.setHide(1);
149 srcnodata_opt.setHide(1);
150 polythreshold_opt.setHide(1);
151 percentile_opt.setHide(1);
152 buffer_opt.setHide(1);
153 disc_opt.setHide(1);
154
155 bool doProcess;//stop process when program was invoked with help option (-h --help)
156 try{
157 doProcess=image_opt.retrieveOption(argc,argv);
158 sample_opt.retrieveOption(argc,argv);
159 layer_opt.retrieveOption(argc,argv);
160 random_opt.retrieveOption(argc,argv);
161 grid_opt.retrieveOption(argc,argv);
162 output_opt.retrieveOption(argc,argv);
163 class_opt.retrieveOption(argc,argv);
164 threshold_opt.retrieveOption(argc,argv);
165 percentile_opt.retrieveOption(argc,argv);
166 ogrformat_opt.retrieveOption(argc,argv);
167 ftype_opt.retrieveOption(argc,argv);
168 band_opt.retrieveOption(argc,argv);
169 bstart_opt.retrieveOption(argc,argv);
170 bend_opt.retrieveOption(argc,argv);
171 rule_opt.retrieveOption(argc,argv);
172 bndnodata_opt.retrieveOption(argc,argv);
173 srcnodata_opt.retrieveOption(argc,argv);
174 polythreshold_opt.retrieveOption(argc,argv);
175 buffer_opt.retrieveOption(argc,argv);
176 disc_opt.retrieveOption(argc,argv);
177 verbose_opt.retrieveOption(argc,argv);
178 }
179 catch(string predefinedString){
180 std::cout << predefinedString << std::endl;
181 exit(0);
182 }
183 if(!doProcess){
184 cout << endl;
185 cout << "Usage: pkextractogr -i input [-s sample | -rand number | -grid size] -o output" << endl;
186 cout << endl;
187 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
188 exit(0);//help was invoked, stop processing
189 }
190
191 std::map<std::string, rule::RULE_TYPE> ruleMap;
192 std::map<std::string, std::string> fieldMap;
193 //initialize ruleMap
194 ruleMap["point"]=rule::point;
195 ruleMap["centroid"]=rule::centroid;
196 ruleMap["mean"]=rule::mean;
197 ruleMap["stdev"]=rule::stdev;
198 ruleMap["median"]=rule::median;
199 ruleMap["proportion"]=rule::proportion;
200 ruleMap["count"]=rule::count;
201 ruleMap["min"]=rule::min;
202 ruleMap["max"]=rule::max;
203 ruleMap["custom"]=rule::custom;
204 ruleMap["mode"]=rule::mode;
205 ruleMap["sum"]=rule::sum;
206 ruleMap["percentile"]=rule::percentile;
207 ruleMap["allpoints"]=rule::allpoints;
208
209 fieldMap["point"]="point";
210 fieldMap["centroid"]="cntrd";
211 fieldMap["mean"]="mean";
212 fieldMap["stdev"]="stdev";
213 fieldMap["median"]="median";
214 fieldMap["proportion"]="prop";
215 fieldMap["count"]="count";
216 fieldMap["min"]="min";
217 fieldMap["max"]="max";
218 fieldMap["custom"]="custom";
219 fieldMap["mode"]="mode";
220 fieldMap["sum"]="sum";
221 fieldMap["percentile"]="perc";
222 fieldMap["allpoints"]="allp";
224 if(srcnodata_opt.size()){
225 while(srcnodata_opt.size()<bndnodata_opt.size())
226 srcnodata_opt.push_back(srcnodata_opt[0]);
227 stat.setNoDataValues(srcnodata_opt);
228 }
229
230 if(verbose_opt[0])
231 std::cout << class_opt << std::endl;
233 unsigned long int nsample=0;
234 unsigned long int ntotalvalid=0;
235 unsigned long int ntotalinvalid=0;
236
237 ImgReaderGdal imgReader;
238 // ImgReaderGdal imgReader;
239 if(image_opt.empty()){
240 std::cerr << "No image dataset provided (use option -i). Use --help for help information";
241 exit(1);
242 }
243 if(output_opt.empty()){
244 std::cerr << "No output dataset provided (use option -o). Use --help for help information";
245 exit(1);
246 }
247 try{
248 imgReader.open(image_opt[0],GA_ReadOnly);
249 }
250 catch(std::string errorstring){
251 std::cout << errorstring << std::endl;
252 exit(1);
253 }
254
255 //check if rule contains allpoints
256 if(find(rule_opt.begin(),rule_opt.end(),"allpoints")!=rule_opt.end()){
257 //allpoints should be the only rule
258 rule_opt.clear();
259 rule_opt.push_back("allpoints");
260 }
261
262 //convert start and end band options to vector of band indexes
263 try{
264 if(bstart_opt.size()){
265 if(bend_opt.size()!=bstart_opt.size()){
266 string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
267 throw(errorstring);
268 }
269 band_opt.clear();
270 for(int ipair=0;ipair<bstart_opt.size();++ipair){
271 if(bend_opt[ipair]<=bstart_opt[ipair]){
272 string errorstring="Error: index for end band must be smaller then start band";
273 throw(errorstring);
274 }
275 for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
276 band_opt.push_back(iband);
277 }
278 }
279 }
280 catch(string error){
281 cerr << error << std::endl;
282 exit(1);
283 }
284
285 int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand();
286 if(class_opt.size()){
287 if(nband>1){
288 cerr << "Warning: using only first band or multiband image" << endl;
289 nband=1;
290 band_opt.clear();
291 band_opt.push_back(0);
292 }
293 }
294
295 if(verbose_opt[0]>1)
296 std::cout << "Number of bands in input image: " << imgReader.nrOfBand() << std::endl;
297
298 OGRFieldType fieldType;
299 int ogr_typecount=11;//hard coded for now!
300 if(verbose_opt[0]>1)
301 std::cout << "field and label types can be: ";
302 for(int iType = 0; iType < ogr_typecount; ++iType){
303 if(verbose_opt[0]>1)
304 std::cout << " " << OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType);
305 if( OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType) != NULL
306 && EQUAL(OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType),
307 ftype_opt[0].c_str()))
308 fieldType=(OGRFieldType) iType;
309 }
310 switch( fieldType ){
311 case OFTInteger:
312 case OFTReal:
313 case OFTRealList:
314 case OFTString:
315 if(verbose_opt[0]>1)
316 std::cout << std::endl << "field type is: " << OGRFieldDefn::GetFieldTypeName(fieldType) << std::endl;
317 break;
318 default:
319 cerr << "field type " << OGRFieldDefn::GetFieldTypeName(fieldType) << " not supported" << std::endl;
320 exit(1);
321 break;
322 }
323
324 const char* pszMessage;
325 void* pProgressArg=NULL;
326 GDALProgressFunc pfnProgress=GDALTermProgress;
327 double progress=0;
328 srand(time(NULL));
329
330 bool sampleIsRaster=false;
331 bool sampleIsVirtual=false;
332
333 ImgReaderOgr sampleReaderOgr;
334 ImgWriterOgr sampleWriterOgr;
335
336 if(sample_opt.size()){
337 try{
338 sampleReaderOgr.open(sample_opt[0]);
339 }
340 catch(string errorString){
341 //todo: sampleIsRaster will not work from GDAL 2.0!!?? (unification of driver for raster and vector datasets)
342 sampleIsRaster=true;
343 }
344 }
345 else{
346 try{
347 sampleWriterOgr.open("/vsimem/virtual",ogrformat_opt[0]);
348 char **papszOptions=NULL;
349 sampleIsVirtual=true;
350
351 if(random_opt.size()){
352 //create simple random sampling within boundary
353 double ulx,uly,lrx,lry;
354 imgReader.getBoundingBox(ulx,uly,lrx,lry);
355 if(random_opt[0]>0)
356 sampleWriterOgr.createLayer("points", imgReader.getProjection(), wkbPoint, papszOptions);
357 OGRPoint pt;
358 unsigned int ipoint;
359 for(ipoint=0;ipoint<random_opt[0];++ipoint){
360 OGRFeature *pointFeature;
361 pointFeature=sampleWriterOgr.createFeature();
362 double theX=ulx+static_cast<double>(rand())/(RAND_MAX)*(lrx-ulx);
363 double theY=uly-static_cast<double>(rand())/(RAND_MAX)*(uly-lry);
364 pt.setX(theX);
365 pt.setY(theY);
366 pointFeature->SetGeometry( &pt );
367 if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
368 string errorString="Failed to create feature in vector dataset";
369 throw(errorString);
370 }
371 // OGRFeature::DestroyFeature(pointFeature);
372 }
373 if(!ipoint){
374 string errorString="Error: no random point created";
375 throw(errorString);
376 }
377 }
378 else if(grid_opt.size()){
379 //create systematic grid of points
380 double ulx,uly,lrx,lry;
381 imgReader.getBoundingBox(ulx,uly,lrx,lry);
382 if(uly-grid_opt[0]/2<lry&&ulx+grid_opt[0]/2>lrx){
383 string errorString="Error: grid distance too large";
384 throw(errorString);
385 }
386 else if(grid_opt[0]>0)
387 sampleWriterOgr.createLayer("points", imgReader.getProjection(), wkbPoint, papszOptions);
388 else{
389 string errorString="Error: grid distance must be strictly positive number";
390 throw(errorString);
391 }
392 OGRPoint pt;
393 unsigned int ipoint=0;
394 for(double theY=uly-grid_opt[0]/2;theY>lry;theY-=grid_opt[0]){
395 for(double theX=ulx+grid_opt[0]/2;theX<lrx;theX+=grid_opt[0]){
396 if(verbose_opt[0]>1)
397 cout << "position: " << theX << " " << theY << endl;
398 OGRFeature *pointFeature;
399 pointFeature=sampleWriterOgr.createFeature();
400 pt.setX(theX);
401 pt.setY(theY);
402 pointFeature->SetGeometry( &pt );
403 if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
404 string errorString="Failed to create feature in vector dataset";
405 throw(errorString);
406 }
407 ++ipoint;
408 // OGRFeature::DestroyFeature(pointFeature);
409 }
410 }
411 if(!ipoint){
412 string errorString="Error: no points created in grid";
413 throw(errorString);
414 }
415 }
416 else{
417 std::cerr << "Error: no sample dataset provided (use option -s). Use --help for help information";
418 exit(1);
419 }
420 sampleWriterOgr.close();
421 sampleReaderOgr.open("/vsimem/virtual");
422 }
423 catch(string errorString){
424 cerr << errorString << endl;
425 exit(1);
426 }
427 }
428
429 if(sampleIsRaster){
430 cerr << "Error: sample must be vector dataset in OGR format";
431 exit(1);
432 }
433 else{//vector dataset
434 if(verbose_opt[0]>1)
435 std::cout << "creating image sample writer " << output_opt[0] << std::endl;
436
437 ImgWriterOgr ogrWriter;
438 double vectords_ulx;
439 double vectords_uly;
440 double vectords_lrx;
441 double vectords_lry;
442 bool calculateSpatialStatistics=false;
443 try{
444 sampleReaderOgr.getExtent(vectords_ulx,vectords_uly,vectords_lrx,vectords_lry);
445 bool hasCoverage=((vectords_ulx < imgReader.getLrx())&&(vectords_lrx > imgReader.getUlx())&&(vectords_lry < imgReader.getUly())&&(vectords_uly > imgReader.getLry()));
446 if(!hasCoverage){
447 ostringstream ess;
448 ess << "No coverage in " << image_opt[0] << " for any layer in " << sample_opt[0] << endl;
449 throw(ess.str());
450 }
451 ogrWriter.open(output_opt[0],ogrformat_opt[0]);
452 //if class_opt not set, get number of classes from input image for these rules
453 for(int irule=0;irule<rule_opt.size();++irule){
454 switch(ruleMap[rule_opt[irule]]){
455 case(rule::point):
456 case(rule::centroid):
457 case(rule::allpoints):
458 break;
459 case(rule::proportion):
460 case(rule::count):
461 case(rule::custom):
462 case(rule::mode):{
463 if(class_opt.empty()){
464 int theBand=0;
465 double minValue=0;
466 double maxValue=0;
467 if(band_opt.size())
468 theBand=band_opt[0];
469 imgReader.getMinMax(minValue,maxValue,theBand);
470 int nclass=maxValue-minValue+1;
471 if(nclass<0&&nclass<256){
472 string errorString="Could not automatically define classes, please set class option";
473 throw(errorString);
474 }
475 for(int iclass=minValue;iclass<=maxValue;++iclass)
476 class_opt.push_back(iclass);
477 }
478 }//deliberate fall through: calculate spatial statistics for all non-point like rules
479 default:
480 calculateSpatialStatistics=true;
481 break;
482 }
483 }
484 }
485 catch(string errorString){
486 cerr << errorString << endl;
487 exit(1);
488 }
489
490 //support multiple layers
491 int nlayerRead=sampleReaderOgr.getDataSource()->GetLayerCount();
492 int ilayerWrite=0;
493 unsigned long int ntotalvalid=0;
494
495 if(verbose_opt[0])
496 std::cout << "number of layers: " << nlayerRead << endl;
497
498 for(int ilayer=0;ilayer<nlayerRead;++ilayer){
499 OGRLayer *readLayer=sampleReaderOgr.getLayer(ilayer);
500 string currentLayername=readLayer->GetName();
501 int layerIndex=ilayer;
502 if(layer_opt.size()){
503 vector<string>::const_iterator it=find(layer_opt.begin(),layer_opt.end(),currentLayername);
504 if(it==layer_opt.end())
505 continue;
506 else
507 layerIndex=it-layer_opt.begin();
508 }
509 double layer_ulx;
510 double layer_uly;
511 double layer_lrx;
512 double layer_lry;
513 sampleReaderOgr.getExtent(layer_ulx,layer_uly,layer_lrx,layer_lry,ilayer);
514 bool hasCoverage=((layer_ulx < imgReader.getLrx())&&(layer_lrx > imgReader.getUlx())&&(layer_lry < imgReader.getUly())&&(layer_uly > imgReader.getLry()));
515 if(!hasCoverage)
516 continue;
517
518 //align bounding box to input image
519 layer_ulx-=fmod(layer_ulx-imgReader.getUlx(),imgReader.getDeltaX());
520 layer_lrx+=fmod(imgReader.getLrx()-layer_lrx,imgReader.getDeltaX());
521 layer_uly+=fmod(imgReader.getUly()-layer_uly,imgReader.getDeltaY());
522 layer_lry-=fmod(layer_lry-imgReader.getLry(),imgReader.getDeltaY());
523
524 //do not read outside input image
525 if(layer_ulx<imgReader.getUlx())
526 layer_ulx=imgReader.getUlx();
527 if(layer_lrx>imgReader.getLrx())
528 layer_lrx=imgReader.getLrx();
529 if(layer_uly>imgReader.getUly())
530 layer_uly=imgReader.getUly();
531 if(layer_lry<imgReader.getLry())
532 layer_lry=imgReader.getLry();
533
534 //read entire block for coverage in memory
535 //todo: use different data types
536 vector< Vector2d<float> > readValuesReal(nband);
537 vector< Vector2d<int> > readValuesInt(nband);
538
539 double layer_uli;
540 double layer_ulj;
541 double layer_lri;
542 double layer_lrj;
543 imgReader.geo2image(layer_ulx,layer_uly,layer_uli,layer_ulj);
544 imgReader.geo2image(layer_lrx,layer_lry,layer_lri,layer_lrj);
545
546 OGRwkbGeometryType layerGeometry=readLayer->GetLayerDefn()->GetGeomType();
547
548 if(layerGeometry==wkbPoint){
549 if(calculateSpatialStatistics){
550 if(buffer_opt[0]<1)
551 buffer_opt[0]=1;
552 }
553 }
554
555 //extend bounding box with buffer
556 if(buffer_opt.size()){
557 layer_uli-=buffer_opt[0];
558 layer_ulj-=buffer_opt[0];
559 layer_lri+=buffer_opt[0];
560 layer_lrj+=buffer_opt[0];
561 }
562
563 //we already checked there is coverage
564 layer_uli=(layer_uli<0)? 0 : static_cast<int>(layer_uli);
565 layer_ulj=(layer_ulj<0)? 0 : static_cast<int>(layer_ulj);
566 layer_lri=(layer_lri>=imgReader.nrOfCol())? imgReader.nrOfCol()-1 : static_cast<int>(layer_lri);
567 layer_lrj=(layer_lrj>=imgReader.nrOfRow())? imgReader.nrOfRow()-1 : static_cast<int>(layer_lrj);
568
569 try{
570 for(int iband=0;iband<nband;++iband){
571 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
572 if(theBand<0){
573 string errorString="Error: illegal band (must be positive and starting from 0)";
574 throw(errorString);
575 }
576 if(theBand>=imgReader.nrOfBand()){
577 string errorString="Error: illegal band (must be lower than number of bands in input raster dataset)";
578 throw(errorString);
579 }
580 if(verbose_opt[0])
581 cout << "reading image band " << theBand << " block rows " << layer_ulj << "-" << layer_lrj << ", cols " << layer_uli << "-" << layer_lri << endl;
582 switch( fieldType ){
583 case OFTInteger:
584 imgReader.readDataBlock(readValuesInt[iband],layer_uli,layer_lri,layer_ulj,layer_lrj,theBand);
585 break;
586 case OFTReal:
587 default:
588 imgReader.readDataBlock(readValuesReal[iband],layer_uli,layer_lri,layer_ulj,layer_lrj,theBand);
589 break;
590 }
591 }
592 }
593 catch(std::string e){
594 std::cout << e << std::endl;
595 exit(1);
596 }
597
598
599 float theThreshold=(threshold_opt.size()==layer_opt.size())? threshold_opt[layerIndex]: threshold_opt[0];
600 cout << "processing layer " << currentLayername << endl;
601
602 bool createPolygon=true;
603 if(find(rule_opt.begin(),rule_opt.end(),"allpoints")!=rule_opt.end())
604 createPolygon=false;
605
606 OGRLayer *writeLayer;
607 if(createPolygon){
608 //create polygon
609 if(verbose_opt[0])
610 std::cout << "create polygons" << std::endl;
611 char **papszOptions=NULL;
612 writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPolygon, papszOptions);
613 }
614 else{
615 if(verbose_opt[0])
616 std::cout << "create points in layer " << readLayer->GetName() << std::endl;
617 char **papszOptions=NULL;
618
619 writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPoint, papszOptions);
620 }
621 if(verbose_opt[0])
622 std::cout << "copy fields from layer " << ilayer << std::flush << std::endl;
623 ogrWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
624
625 for(int irule=0;irule<rule_opt.size();++irule){
626 for(int iband=0;iband<nband;++iband){
627 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
628 ostringstream fs;
629 if(rule_opt.size()>1||nband==1)
630 fs << fieldMap[rule_opt[irule]];
631 if(nband>1)
632 fs << "b" << theBand;
633 switch(ruleMap[rule_opt[irule]]){
634 case(rule::proportion):
635 case(rule::count):{//count for each class
636 for(int iclass=0;iclass<class_opt.size();++iclass){
637 ostringstream fsclass;
638 fsclass << fs.str() << class_opt[iclass];
639 ogrWriter.createField(fsclass.str(),fieldType,ilayerWrite);
640 }
641 break;
642 }
643 case(rule::percentile):{//for each percentile
644 for(int iperc=0;iperc<percentile_opt.size();++iperc){
645 ostringstream fsperc;
646 fsperc << fs.str() << percentile_opt[iperc];
647 ogrWriter.createField(fsperc.str(),fieldType,ilayerWrite);
648 }
649 break;
650 }
651 default:
652 ogrWriter.createField(fs.str(),fieldType,ilayerWrite);
653 break;
654 }
655 }
656 }
657 OGRFeature *readFeature;
658 unsigned long int ifeature=0;
659 unsigned long int nfeatureLayer=sampleReaderOgr.getFeatureCount(ilayer);
660 unsigned long int ntotalvalidLayer=0;
661
662 if(nfeatureLayer<=0)
663 continue;
664 progress=0;
665 pfnProgress(progress,pszMessage,pProgressArg);
666 readLayer->ResetReading();
667 while( (readFeature = readLayer->GetNextFeature()) != NULL ){
668 bool validFeature=false;
669 if(verbose_opt[0]>2)
670 std::cout << "reading feature " << readFeature->GetFID() << std::endl;
671 if(theThreshold>0){//percentual value
672 double p=static_cast<double>(rand())/(RAND_MAX);
673 p*=100.0;
674 if(p>theThreshold){
675 continue;//do not select for now, go to next feature
676 }
677 }
678 else{//absolute value
679 if(threshold_opt.size()==layer_opt.size()){
680 if(ntotalvalidLayer>=-theThreshold){
681 continue;//do not select any more pixels, go to next column feature
682 }
683 }
684 else{
685 if(ntotalvalid>=-theThreshold){
686 continue;//do not select any more pixels, go to next column feature
687 }
688 }
689 }
690 if(verbose_opt[0]>2)
691 std::cout << "processing feature " << readFeature->GetFID() << std::endl;
692 //get x and y from readFeature
693 // double x,y;
694 OGRGeometry *poGeometry;
695 poGeometry = readFeature->GetGeometryRef();
696 assert(poGeometry!=NULL);
697 try{
698 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ){
699 OGRPoint readPoint = *((OGRPoint *) poGeometry);
700
701 double i_centre,j_centre;
702 imgReader.geo2image(readPoint.getX(),readPoint.getY(),i_centre,j_centre);
703 //nearest neighbour
704 j_centre=static_cast<int>(j_centre);
705 i_centre=static_cast<int>(i_centre);
706
707 double uli=i_centre-buffer_opt[0];
708 double ulj=j_centre-buffer_opt[0];
709 double lri=i_centre+buffer_opt[0];
710 double lrj=j_centre+buffer_opt[0];
711
712 //nearest neighbour
713 ulj=static_cast<int>(ulj);
714 uli=static_cast<int>(uli);
715 lrj=static_cast<int>(lrj);
716 lri=static_cast<int>(lri);
717
718 //check if j is out of bounds
719 if(static_cast<int>(ulj)<0||static_cast<int>(ulj)>=imgReader.nrOfRow())
720 continue;
721 //check if j is out of bounds
722 if(static_cast<int>(uli)<0||static_cast<int>(lri)>=imgReader.nrOfCol())
723 continue;
724
725 OGRPoint ulPoint,urPoint,llPoint,lrPoint;
726 double ulx;
727 double uly;
728 double lrx;
729 double lry;
730
731 OGRPolygon writePolygon;
732 OGRPoint writePoint;
733 OGRLinearRing writeRing;
734 OGRFeature *writePolygonFeature;
735
736 int nPointPolygon=0;
737 if(createPolygon){
738 if(disc_opt[0]){
739 double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
740 double radius=buffer_opt[0]*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
741 unsigned short nstep = 25;
742 for(int i=0;i<nstep;++i){
743 OGRPoint aPoint;
744 aPoint.setX(readPoint.getX()+imgReader.getDeltaX()/2.0+radius*cos(2*PI*i/nstep));
745 aPoint.setY(readPoint.getY()-imgReader.getDeltaY()/2.0+radius*sin(2*PI*i/nstep));
746 writeRing.addPoint(&aPoint);
747 }
748 writePolygon.addRing(&writeRing);
749 writePolygon.closeRings();
750 }
751 else{
752 double ulx,uly,lrx,lry;
753 imgReader.image2geo(uli,ulj,ulx,uly);
754 imgReader.image2geo(lri,lrj,lrx,lry);
755 ulPoint.setX(ulx-imgReader.getDeltaX()/2.0);
756 ulPoint.setY(uly+imgReader.getDeltaY()/2.0);
757 lrPoint.setX(lrx+imgReader.getDeltaX()/2.0);
758 lrPoint.setY(lry-imgReader.getDeltaY()/2.0);
759 urPoint.setX(lrx+imgReader.getDeltaX()/2.0);
760 urPoint.setY(uly+imgReader.getDeltaY()/2.0);
761 llPoint.setX(ulx-imgReader.getDeltaX()/2.0);
762 llPoint.setY(lry-imgReader.getDeltaY()/2.0);
763
764 writeRing.addPoint(&ulPoint);
765 writeRing.addPoint(&urPoint);
766 writeRing.addPoint(&lrPoint);
767 writeRing.addPoint(&llPoint);
768 writePolygon.addRing(&writeRing);
769 writePolygon.closeRings();
770 }
771 writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
772 if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
773 cerr << "writing feature failed" << std::endl;
774 writePolygonFeature->SetGeometry(&writePolygon);
775 if(verbose_opt[0]>1)
776 std::cout << "copying new fields write polygon " << std::endl;
777 if(verbose_opt[0]>1)
778 std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
779
780 OGRPoint readPoint;
781 if(find(rule_opt.begin(),rule_opt.end(),"centroid")!=rule_opt.end()){
782 if(verbose_opt[0]>1)
783 std::cout << "get centroid" << std::endl;
784 writePolygon.Centroid(&readPoint);
785 double i,j;
786 imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
787 int indexJ=static_cast<int>(j-layer_ulj);
788 int indexI=static_cast<int>(i-layer_uli);
789 bool valid=true;
790 valid=valid&&(indexJ>=0);
791 valid=valid&&(indexJ<imgReader.nrOfRow());
792 valid=valid&&(indexI>=0);
793 valid=valid&&(indexI<imgReader.nrOfCol());
794
795 if(valid){
796 if(srcnodata_opt.empty())
797 validFeature=true;
798 else{
799 for(int vband=0;vband<bndnodata_opt.size();++vband){
800 switch( fieldType ){
801 case OFTInteger:{
802 int value;
803 value=((readValuesInt[vband])[indexJ])[indexI];
804 if(value==srcnodata_opt[vband])
805 valid=false;
806 break;
807 }
808 case OFTReal:{
809 double value;
810 value=((readValuesReal[vband])[indexJ])[indexI];
811 if(value==srcnodata_opt[vband])
812 valid=false;
813 break;
814 }
815 }
816 if(!valid)
817 continue;
818 else
819 validFeature=true;
820 }
821 }
822 }
823 if(valid){
824 assert(readValuesReal.size()==nband);
825 for(int iband=0;iband<nband;++iband){
826 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
827 //write fields for point on surface and centroid
828 string fieldname;
829 ostringstream fs;
830 if(rule_opt.size()>1||nband==1)
831 fs << fieldMap["centroid"];
832 if(nband>1)
833 fs << "b" << theBand;
834 fieldname=fs.str();
835 switch( fieldType ){
836 case OFTInteger:
837 writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(((readValuesInt[iband])[indexJ])[indexI]));
838 break;
839 case OFTReal:
840 writePolygonFeature->SetField(fieldname.c_str(),((readValuesReal[iband])[indexJ])[indexI]);
841 break;
842 default://not supported
843 std::string errorString="field type not supported";
844 throw(errorString);
845 break;
846 }
847 }
848 }
849 }//if centroid
850 if(find(rule_opt.begin(),rule_opt.end(),"point")!=rule_opt.end()){
851 if(verbose_opt[0]>1)
852 std::cout << "get point on surface" << std::endl;
853 if(writePolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
854 writePolygon.Centroid(&readPoint);
855 double i,j;
856 imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
857 int indexJ=static_cast<int>(j-layer_ulj);
858 int indexI=static_cast<int>(i-layer_uli);
859 bool valid=true;
860 valid=valid&&(indexJ>=0);
861 valid=valid&&(indexJ<imgReader.nrOfRow());
862 valid=valid&&(indexI>=0);
863 valid=valid&&(indexI<imgReader.nrOfCol());
864 if(valid&&srcnodata_opt.size()){
865 for(int vband=0;vband<bndnodata_opt.size();++vband){
866 switch( fieldType ){
867 case OFTInteger:{
868 int value;
869 value=((readValuesInt[vband])[indexJ])[indexI];
870 if(value==srcnodata_opt[vband])
871 valid=false;
872 break;
873 }
874 case OFTReal:{
875 double value;
876 value=((readValuesReal[vband])[indexJ])[indexI];
877 if(value==srcnodata_opt[vband])
878 valid=false;
879 break;
880 }
881 }
882 if(!valid)
883 continue;
884 else
885 validFeature=true;
886 }
887 }
888 if(valid){
889 for(int iband=0;iband<nband;++iband){
890 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
891 //write fields for point on surface and centroid
892 string fieldname;
893 ostringstream fs;
894 if(rule_opt.size()>1||nband==1)
895 fs << fieldMap["point"];
896 if(nband>1)
897 fs << "b" << theBand;
898 fieldname=fs.str();
899 switch( fieldType ){
900 case OFTInteger:
901 writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(((readValuesInt[iband])[indexJ])[indexI]));
902 break;
903 case OFTReal:
904 writePolygonFeature->SetField(fieldname.c_str(),((readValuesReal[iband])[indexJ])[indexI]);
905 break;
906 default://not supported
907 std::string errorString="field type not supported";
908 throw(errorString);
909 break;
910 }
911 }
912 }
913 }//if point
914 }//if createPolygon
915
916 if(calculateSpatialStatistics||!createPolygon){
917 Vector2d<double> polyValues;
918 vector<double> polyClassValues;
919
920 if(class_opt.size()){
921 polyClassValues.resize(class_opt.size());
922 //initialize
923 for(int iclass=0;iclass<class_opt.size();++iclass)
924 polyClassValues[iclass]=0;
925 }
926 polyValues.resize(nband);
927
928 OGRPoint thePoint;
929 for(int j=ulj;j<=lrj;++j){
930 for(int i=uli;i<=lri;++i){
931 //check if within raster image
932 if(i<0||i>=imgReader.nrOfCol())
933 continue;
934 if(j<0||j>=imgReader.nrOfRow())
935 continue;
936 int indexJ=j-layer_ulj;
937 int indexI=i-layer_uli;
938 if(indexJ<0)
939 indexJ=0;
940 if(indexI<0)
941 indexI=0;
942 if(indexJ>=imgReader.nrOfRow())
943 indexJ=imgReader.nrOfRow()-1;
944 if(indexI>=imgReader.nrOfCol())
945 indexI=imgReader.nrOfCol()-1;
946
947 double theX=0;
948 double theY=0;
949 imgReader.image2geo(i,j,theX,theY);
950 thePoint.setX(theX);
951 thePoint.setY(theY);
952 if(disc_opt[0]&&buffer_opt[0]>0){
953 double radius=buffer_opt[0]*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
954 if((theX-readPoint.getX())*(theX-readPoint.getX())+(theY-readPoint.getY())*(theY-readPoint.getY())>radius*radius)
955 continue;
956 }
957 bool valid=true;
958
959 if(srcnodata_opt.size()){
960 for(int vband=0;vband<bndnodata_opt.size();++vband){
961 switch( fieldType ){
962 case OFTInteger:{
963 int value=((readValuesInt[vband])[indexJ])[indexI];
964 if(value==srcnodata_opt[vband]){
965 valid=false;
966 }
967 break;
968 }
969 default:{
970 float value=((readValuesReal[vband])[indexJ])[indexI];
971 if(value==srcnodata_opt[vband]){
972 valid=false;
973 }
974 break;
975 }
976 }
977 }
978 }
979 if(!valid)
980 continue;
981 else
982 validFeature=true;
983
984 ++nPointPolygon;
985 OGRFeature *writePointFeature;
986 if(valid&&!createPolygon){//write all points
987 if(polythreshold_opt.size())
988 if(nPointPolygon>polythreshold_opt[0])
989 break;
990 //create feature
991 writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
992 if(verbose_opt[0]>1)
993 std::cout << "copying fields from point feature " << std::endl;
994 if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
995 cerr << "writing feature failed" << std::endl;
996 if(verbose_opt[0]>1)
997 std::cout << "set geometry as point " << std::endl;
998 writePointFeature->SetGeometry(&thePoint);
999 assert(wkbFlatten(writePointFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
1000 if(verbose_opt[0]>1){
1001 std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields:" << std::endl;
1002 for(int iField=0;iField<writePointFeature->GetFieldCount();++iField){
1003 std::string fieldname=writeLayer->GetLayerDefn()->GetFieldDefn(iField)->GetNameRef();
1004 cout << fieldname << endl;
1005 }
1006 }
1007 }
1008 // if(valid&&class_opt.size()){
1009 // short value=0;
1010 // switch( fieldType ){
1011 // case OFTInteger:
1012 // value=((readValuesInt[0])[indexJ])[indexI];
1013 // break;
1014 // case OFTReal:
1015 // value=((readValuesReal[0])[indexJ])[indexI];
1016 // break;
1017 // }
1018 // for(int iclass=0;iclass<class_opt.size();++iclass){
1019 // if(value==class_opt[iclass])
1020 // polyClassValues[iclass]+=1;
1021 // }
1022 // }
1023 if(valid){
1024 for(int iband=0;iband<nband;++iband){
1025 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1026 double value=0;
1027 switch( fieldType ){
1028 case OFTInteger:
1029 value=((readValuesInt[iband])[indexJ])[indexI];
1030 break;
1031 case OFTReal:
1032 value=((readValuesReal[iband])[indexJ])[indexI];
1033 break;
1034 }
1035 if(!iband&&class_opt.size()){
1036 for(int iclass=0;iclass<class_opt.size();++iclass){
1037 if(value==class_opt[iclass])
1038 polyClassValues[iclass]+=1;
1039 }
1040 }
1041 if(verbose_opt[0]>1)
1042 std::cout << ": " << value << std::endl;
1043 if(!createPolygon){//write all points within polygon
1044 string fieldname;
1045 ostringstream fs;
1046 if(rule_opt.size()>1||nband==1)
1047 fs << fieldMap["allpoints"];
1048 if(nband>1)
1049 fs << "b" << theBand;
1050 fieldname=fs.str();
1051 int fieldIndex=writePointFeature->GetFieldIndex(fieldname.c_str());
1052 if(fieldIndex<0){
1053 cerr << "field " << fieldname << " was not found" << endl;
1054 exit(1);
1055 }
1056 if(verbose_opt[0]>1)
1057 std::cout << "set field " << fieldname << " to " << value << std::endl;
1058 switch( fieldType ){
1059 case OFTInteger:
1060 case OFTReal:
1061 writePointFeature->SetField(fieldname.c_str(),value);
1062 break;
1063 default://not supported
1064 assert(0);
1065 break;
1066 }
1067 }
1068 else{
1069 polyValues[iband].push_back(value);
1070 }
1071 }//iband
1072 }
1073 if(valid&&!createPolygon){
1074 //write feature
1075 if(verbose_opt[0]>1)
1076 std::cout << "creating point feature" << std::endl;
1077 if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1078 std::string errorString="Failed to create feature in ogr vector dataset";
1079 throw(errorString);
1080 }
1081 //destroy feature
1082 // OGRFeature::DestroyFeature( writePointFeature );
1083 ++ntotalvalid;
1084 ++ntotalvalidLayer;
1085 }
1086 }//for in i
1087 }//for int j
1088
1089 if(createPolygon){
1090 //do not create if no points found within polygon
1091 if(!nPointPolygon){
1092 if(verbose_opt[0])
1093 cout << "no points found in polygon, continuing" << endl;
1094 continue;
1095 }
1096 //write field attributes to polygon feature
1097 for(int irule=0;irule<rule_opt.size();++irule){
1098 //skip centroid and point
1099 if(ruleMap[rule_opt[irule]]==rule::centroid||ruleMap[rule_opt[irule]]==rule::point)
1100 continue;
1101 for(int iband=0;iband<nband;++iband){
1102 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1103 vector<double> theValue;
1104 vector<string> fieldname;
1105 ostringstream fs;
1106 if(rule_opt.size()>1||nband==1)
1107 fs << fieldMap[rule_opt[irule]];
1108 if(nband>1)
1109 fs << "b" << theBand;
1110 switch(ruleMap[rule_opt[irule]]){
1111 case(rule::proportion):
1112 stat.normalize_pct(polyClassValues);
1113 case(rule::count):{//count for each class
1114 for(int index=0;index<polyClassValues.size();++index){
1115 theValue.push_back(polyClassValues[index]);
1116 ostringstream fsclass;
1117 fsclass << fs.str() << class_opt[index];
1118 fieldname.push_back(fsclass.str());
1119 }
1120 break;
1121 }
1122 case(rule::mode):{
1123 //maximum votes in polygon
1124 if(verbose_opt[0])
1125 std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1126 //search for class with maximum votes
1127 int maxClass=stat.mymin(class_opt);
1128 vector<double>::iterator maxit;
1129 maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
1130 int maxIndex=distance(polyClassValues.begin(),maxit);
1131 maxClass=class_opt[maxIndex];
1132 if(verbose_opt[0]>0)
1133 std::cout << "maxClass: " << maxClass << std::endl;
1134 theValue.push_back(maxClass);
1135 fieldname.push_back(fs.str());
1136 }
1137 case(rule::mean):
1138 theValue.push_back(stat.mean(polyValues[iband]));
1139 fieldname.push_back(fs.str());
1140 break;
1141 case(rule::median):
1142 theValue.push_back(stat.median(polyValues[iband]));
1143 fieldname.push_back(fs.str());
1144 break;
1145 case(rule::stdev):
1146 theValue.push_back(sqrt(stat.var(polyValues[iband])));
1147 fieldname.push_back(fs.str());
1148 break;
1149 case(rule::percentile):{
1150 for(int iperc=0;iperc<percentile_opt.size();++iperc){
1151 theValue.push_back(stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[iperc]));
1152 ostringstream fsperc;
1153 fsperc << fs.str() << percentile_opt[iperc];
1154 fieldname.push_back(fsperc.str());
1155 }
1156 break;
1157 }
1158 case(rule::sum):
1159 theValue.push_back(stat.sum(polyValues[iband]));
1160 fieldname.push_back(fs.str());
1161 break;
1162 case(rule::max):
1163 theValue.push_back(stat.mymax(polyValues[iband]));
1164 fieldname.push_back(fs.str());
1165 break;
1166 case(rule::min):
1167 theValue.push_back(stat.mymin(polyValues[iband]));
1168 fieldname.push_back(fs.str());
1169 break;
1170 case(rule::centroid):
1171 theValue.push_back(polyValues[iband].back());
1172 fieldname.push_back(fs.str());
1173 break;
1174 default://not supported
1175 break;
1176 }
1177 for(int ivalue=0;ivalue<theValue.size();++ivalue){
1178 switch( fieldType ){
1179 case OFTInteger:
1180 writePolygonFeature->SetField(fieldname[ivalue].c_str(),static_cast<int>(theValue[ivalue]));
1181 break;
1182 case OFTReal:
1183 writePolygonFeature->SetField(fieldname[ivalue].c_str(),theValue[ivalue]);
1184 break;
1185 case OFTString:
1186 writePolygonFeature->SetField(fieldname[ivalue].c_str(),type2string<double>(theValue[ivalue]).c_str());
1187 break;
1188 default://not supported
1189 std::string errorString="field type not supported";
1190 throw(errorString);
1191 break;
1192 }
1193 }
1194 }
1195 }
1196 }
1197 }
1198 //test
1199 if(createPolygon&&validFeature){
1200 // if(createPolygon){
1201 //write polygon feature
1202 //todo: create only in case of valid feature
1203 if(verbose_opt[0]>1)
1204 std::cout << "creating polygon feature (1)" << std::endl;
1205 if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1206 std::string errorString="Failed to create polygon feature in ogr vector dataset";
1207 throw(errorString);
1208 }
1209 //test: no need to destroy anymore?
1210 // OGRFeature::DestroyFeature( writePolygonFeature );
1211 ++ntotalvalid;
1212 ++ntotalvalidLayer;
1213 }
1214 }
1215 else{
1216 OGRPolygon readPolygon;
1217 OGRMultiPolygon readMultiPolygon;
1218
1219 //get envelope
1220 OGREnvelope* psEnvelope=new OGREnvelope();
1221
1222 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
1223 readPolygon = *((OGRPolygon *) poGeometry);
1224 readPolygon.closeRings();
1225 readPolygon.getEnvelope(psEnvelope);
1226 }
1227 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
1228 readMultiPolygon = *((OGRMultiPolygon *) poGeometry);
1229 readMultiPolygon.closeRings();
1230 readMultiPolygon.getEnvelope(psEnvelope);
1231 }
1232 else{
1233 std::string test;
1234 test=poGeometry->getGeometryName();
1235 ostringstream oss;
1236 oss << "geometry " << test << " not supported";
1237 throw(oss.str());
1238 }
1239
1240 double ulx,uly,lrx,lry;
1241 double uli,ulj,lri,lrj;
1242 ulx=psEnvelope->MinX;
1243 uly=psEnvelope->MaxY;
1244 lrx=psEnvelope->MaxX;
1245 lry=psEnvelope->MinY;
1246 delete psEnvelope;
1247
1248 //check if feature is covered by input raster dataset
1249 if(!imgReader.covers(ulx,uly,lrx,lry))
1250 continue;
1251
1252 OGRFeature *writePolygonFeature;
1253 int nPointPolygon=0;
1254 if(createPolygon){
1255 writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1256 writePolygonFeature->SetGeometry(poGeometry);
1257 //writePolygonFeature and readFeature are both of type wkbPolygon
1258 if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
1259 cerr << "writing feature failed" << std::endl;
1260 if(verbose_opt[0]>1)
1261 std::cout << "copying new fields write polygon " << std::endl;
1262 if(verbose_opt[0]>1)
1263 std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
1264 }
1265
1266 OGRPoint readPoint;
1267 if(find(rule_opt.begin(),rule_opt.end(),"centroid")!=rule_opt.end()){
1268 if(verbose_opt[0]>1)
1269 std::cout << "get centroid" << std::endl;
1270 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)
1271 readPolygon.Centroid(&readPoint);
1272 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon)
1273 readMultiPolygon.Centroid(&readPoint);
1274
1275 double i,j;
1276 imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
1277 int indexJ=static_cast<int>(j-layer_ulj);
1278 int indexI=static_cast<int>(i-layer_uli);
1279 bool valid=true;
1280 valid=valid&&(indexJ>=0);
1281 valid=valid&&(indexJ<imgReader.nrOfRow());
1282 valid=valid&&(indexI>=0);
1283 valid=valid&&(indexI<imgReader.nrOfCol());
1284 if(valid){
1285 if(srcnodata_opt.empty())
1286 validFeature=true;
1287 else{
1288 for(int vband=0;vband<bndnodata_opt.size();++vband){
1289 switch( fieldType ){
1290 case OFTInteger:{
1291 int value;
1292 value=((readValuesInt[vband])[indexJ])[indexI];
1293 if(value==srcnodata_opt[vband])
1294 valid=false;
1295 break;
1296 }
1297 case OFTReal:{
1298 double value;
1299 value=((readValuesReal[vband])[indexJ])[indexI];
1300 if(value==srcnodata_opt[vband])
1301 valid=false;
1302 break;
1303 }
1304 }
1305 if(!valid)
1306 continue;
1307 else
1308 validFeature=true;
1309 }
1310 }
1311 }
1312 if(valid){
1313 assert(readValuesReal.size()==nband);
1314 for(int iband=0;iband<nband;++iband){
1315 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1316 //write fields for point on surface and centroid
1317 string fieldname;
1318 ostringstream fs;
1319 if(rule_opt.size()>1||nband==1)
1320 fs << fieldMap["centroid"];
1321 if(nband>1)
1322 fs << "b" << theBand;
1323 fieldname=fs.str();
1324 switch( fieldType ){
1325 case OFTInteger:
1326 writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(((readValuesInt[iband])[indexJ])[indexI]));
1327 break;
1328 case OFTReal:
1329 writePolygonFeature->SetField(fieldname.c_str(),((readValuesReal[iband])[indexJ])[indexI]);
1330 break;
1331 default://not supported
1332 std::string errorString="field type not supported";
1333 throw(errorString);
1334 break;
1335 }
1336 }
1337 }
1338 }
1339 if(find(rule_opt.begin(),rule_opt.end(),"point")!=rule_opt.end()){
1340 if(verbose_opt[0]>1)
1341 std::cout << "get point on surface" << std::endl;
1342 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
1343 if(readPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
1344 readPolygon.Centroid(&readPoint);
1345 }
1346 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
1347 // if(readMultiPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
1348 readMultiPolygon.Centroid(&readPoint);
1349 }
1350 double i,j;
1351 imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
1352 int indexJ=static_cast<int>(j-layer_ulj);
1353 int indexI=static_cast<int>(i-layer_uli);
1354 bool valid=true;
1355 valid=valid&&(indexJ>=0);
1356 valid=valid&&(indexJ<imgReader.nrOfRow());
1357 valid=valid&&(indexI>=0);
1358 valid=valid&&(indexI<imgReader.nrOfCol());
1359 if(valid&&srcnodata_opt.size()){
1360 for(int vband=0;vband<bndnodata_opt.size();++vband){
1361 switch( fieldType ){
1362 case OFTInteger:{
1363 int value;
1364 value=((readValuesInt[vband])[indexJ])[indexI];
1365 if(value==srcnodata_opt[vband])
1366 valid=false;
1367 break;
1368 }
1369 case OFTReal:{
1370 double value;
1371 value=((readValuesReal[vband])[indexJ])[indexI];
1372 if(value==srcnodata_opt[vband])
1373 valid=false;
1374 break;
1375 }
1376 }
1377 if(!valid)
1378 continue;
1379 else
1380 validFeature=true;
1381 }
1382 }
1383 if(valid){
1384 for(int iband=0;iband<nband;++iband){
1385 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1386 //write fields for point on surface and centroid
1387 string fieldname;
1388 ostringstream fs;
1389 if(rule_opt.size()>1||nband==1)
1390 fs << fieldMap["point"];
1391 if(nband>1)
1392 fs << "b" << theBand;
1393 fieldname=fs.str();
1394 switch( fieldType ){
1395 case OFTInteger:
1396 writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(((readValuesInt[iband])[indexJ])[indexI]));
1397 break;
1398 case OFTReal:
1399 writePolygonFeature->SetField(fieldname.c_str(),((readValuesReal[iband])[indexJ])[indexI]);
1400 break;
1401 default://not supported
1402 std::string errorString="field type not supported";
1403 throw(errorString);
1404 break;
1405 }
1406 }
1407 }
1408 }
1409 if(calculateSpatialStatistics||ruleMap[rule_opt[0]]==rule::allpoints){
1410 imgReader.geo2image(ulx,uly,uli,ulj);
1411 imgReader.geo2image(lrx,lry,lri,lrj);
1412 //nearest neighbour
1413 ulj=static_cast<int>(ulj);
1414 uli=static_cast<int>(uli);
1415 lrj=static_cast<int>(lrj);
1416 lri=static_cast<int>(lri);
1417 //iterate through all pixels
1418 if(verbose_opt[0]>1)
1419 std::cout << "bounding box for polygon feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
1420
1421 if(uli<0)
1422 uli=0;
1423 if(lri<0)
1424 lri=0;
1425 if(uli>=imgReader.nrOfCol())
1426 uli=imgReader.nrOfCol()-1;
1427 if(lri>=imgReader.nrOfCol())
1428 lri=imgReader.nrOfCol()-1;
1429 if(ulj<0)
1430 ulj=0;
1431 if(lrj<0)
1432 lrj=0;
1433 if(ulj>=imgReader.nrOfRow())
1434 ulj=imgReader.nrOfRow()-1;
1435 if(lrj>=imgReader.nrOfRow())
1436 lrj=imgReader.nrOfRow()-1;
1437
1438 Vector2d<double> polyValues;
1439 vector<double> polyClassValues;
1440
1441 if(class_opt.size()){
1442 polyClassValues.resize(class_opt.size());
1443 //initialize
1444 for(int iclass=0;iclass<class_opt.size();++iclass)
1445 polyClassValues[iclass]=0;
1446 }
1447 polyValues.resize(nband);
1448
1449 OGRPoint thePoint;
1450 for(int j=ulj;j<=lrj;++j){
1451 for(int i=uli;i<=lri;++i){
1452 //check if within raster image
1453 if(i<0||i>=imgReader.nrOfCol())
1454 continue;
1455 if(j<0||j>=imgReader.nrOfRow())
1456 continue;
1457 int indexJ=j-layer_ulj;
1458 int indexI=i-layer_uli;
1459
1460 double theX=0;
1461 double theY=0;
1462 imgReader.image2geo(i,j,theX,theY);
1463 thePoint.setX(theX);
1464 thePoint.setY(theY);
1465 //check if point is on surface
1466 if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
1467 if(!readPolygon.Contains(&thePoint))
1468 continue;
1469 }
1470 else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){
1471 if(!readMultiPolygon.Contains(&thePoint))
1472 continue;
1473 }
1474
1475 bool valid=true;
1476 if(srcnodata_opt.size()){
1477 for(int vband=0;vband<bndnodata_opt.size();++vband){
1478 switch( fieldType ){
1479 case OFTInteger:{
1480 int value=((readValuesInt[vband])[indexJ])[indexI];
1481 if(value==srcnodata_opt[vband]){
1482 valid=false;
1483 }
1484 break;
1485 }
1486 default:{
1487 float value=((readValuesReal[vband])[indexJ])[indexI];
1488 if(value==srcnodata_opt[vband]){
1489 valid=false;
1490 }
1491 break;
1492 }
1493 }
1494 }
1495 }
1496 if(!valid)
1497 continue;
1498 else
1499 validFeature=true;
1500
1501 if(verbose_opt[0]>1)
1502 std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
1503 ++nPointPolygon;
1504
1505 OGRFeature *writePointFeature;
1506 if(!createPolygon){//write all points within polygon
1507 if(polythreshold_opt.size())
1508 if(nPointPolygon>polythreshold_opt[0])
1509 break;
1510 //create feature
1511 writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1512 if(verbose_opt[0]>1)
1513 std::cout << "copying fields from polygons " << std::endl;
1514 if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
1515 cerr << "writing feature failed" << std::endl;
1516 if(verbose_opt[0]>1)
1517 std::cout << "set geometry as point " << std::endl;
1518 writePointFeature->SetGeometry(&thePoint);
1519 assert(wkbFlatten(writePointFeature->GetGeometryRef()->getGeometryType()) == wkbPoint);
1520 if(verbose_opt[0]>1){
1521 std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields:" << std::endl;
1522 for(int iField=0;iField<writePointFeature->GetFieldCount();++iField){
1523 std::string fieldname=writeLayer->GetLayerDefn()->GetFieldDefn(iField)->GetNameRef();
1524 cout << fieldname << endl;
1525 }
1526 }
1527 }
1528 // if(class_opt.size()){
1529 // short value=0;
1530 // switch( fieldType ){
1531 // case OFTInteger:
1532 // value=((readValuesInt[0])[indexJ])[indexI];
1533 // break;
1534 // case OFTReal:
1535 // value=((readValuesReal[0])[indexJ])[indexI];
1536 // break;
1537 // }
1538 // for(int iclass=0;iclass<class_opt.size();++iclass){
1539 // if(value==class_opt[iclass])
1540 // polyClassValues[iclass]+=1;
1541 // }
1542 // }
1543 // else{
1544 for(int iband=0;iband<nband;++iband){
1545 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1546 double value=0;
1547 switch( fieldType ){
1548 case OFTInteger:
1549 value=((readValuesInt[iband])[indexJ])[indexI];
1550 break;
1551 case OFTReal:
1552 value=((readValuesReal[iband])[indexJ])[indexI];
1553 break;
1554 }
1555
1556 if(!iband&&class_opt.size()){
1557 for(int iclass=0;iclass<class_opt.size();++iclass){
1558 if(value==class_opt[iclass])
1559 polyClassValues[iclass]+=1;
1560 }
1561 }
1562 if(verbose_opt[0]>1)
1563 std::cout << ": " << value << std::endl;
1564 if(!createPolygon){//write all points within polygon
1565 string fieldname;
1566 ostringstream fs;
1567 if(rule_opt.size()>1||nband==1)
1568 fs << fieldMap["allpoints"];
1569 if(nband>1)
1570 fs << "b" << theBand;
1571 fieldname=fs.str();
1572 int fieldIndex=writePointFeature->GetFieldIndex(fieldname.c_str());
1573 if(fieldIndex<0){
1574 cerr << "field " << fieldname << " was not found" << endl;
1575 exit(1);
1576 }
1577 if(verbose_opt[0]>1)
1578 std::cout << "set field " << fieldname << " to " << value << std::endl;
1579 switch( fieldType ){
1580 case OFTInteger:
1581 case OFTReal:
1582 writePointFeature->SetField(fieldname.c_str(),value);
1583 break;
1584 default://not supported
1585 assert(0);
1586 break;
1587 }
1588 }
1589 else{
1590 polyValues[iband].push_back(value);
1591 }
1592 }//iband
1593 if(!createPolygon){
1594 //todo: only if valid feature?
1595 //write feature
1596 if(verbose_opt[0]>1)
1597 std::cout << "creating point feature" << std::endl;
1598 if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1599 std::string errorString="Failed to create feature in ogr vector dataset";
1600 throw(errorString);
1601 }
1602 //destroy feature
1603 // OGRFeature::DestroyFeature( writePointFeature );
1604 ++ntotalvalid;
1605 ++ntotalvalidLayer;
1606 }
1607 }//for in i
1608 }//for int j
1609 if(createPolygon){
1610 //do not create if no points found within polygon
1611 if(!nPointPolygon){
1612 if(verbose_opt[0])
1613 cout << "no points found in polygon, continuing" << endl;
1614 continue;
1615 }
1616 //write field attributes to polygon feature
1617 for(int irule=0;irule<rule_opt.size();++irule){
1618 //skip centroid and point
1619 if(ruleMap[rule_opt[irule]]==rule::centroid||ruleMap[rule_opt[irule]]==rule::point)
1620 continue;
1621 for(int iband=0;iband<nband;++iband){
1622 int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1623 vector<double> theValue;
1624 vector<string> fieldname;
1625 ostringstream fs;
1626 if(rule_opt.size()>1||nband==1)
1627 fs << fieldMap[rule_opt[irule]];
1628 if(nband>1)
1629 fs << "b" << theBand;
1630 switch(ruleMap[rule_opt[irule]]){
1631 case(rule::proportion):
1632 stat.normalize_pct(polyClassValues);
1633 case(rule::count):{//count for each class
1634 for(int index=0;index<polyClassValues.size();++index){
1635 theValue.push_back(polyClassValues[index]);
1636 ostringstream fsclass;
1637 fsclass << fs.str() << class_opt[index];
1638 fieldname.push_back(fsclass.str());
1639 }
1640 break;
1641 }
1642 case(rule::mode):{
1643 //maximum votes in polygon
1644 if(verbose_opt[0])
1645 std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1646 //search for class with maximum votes
1647 int maxClass=stat.mymin(class_opt);
1648 vector<double>::iterator maxit;
1649 maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
1650 int maxIndex=distance(polyClassValues.begin(),maxit);
1651 maxClass=class_opt[maxIndex];
1652 if(verbose_opt[0]>0)
1653 std::cout << "maxClass: " << maxClass << std::endl;
1654 theValue.push_back(maxClass);
1655 fieldname.push_back(fs.str());
1656 break;
1657 }
1658 case(rule::mean):
1659 theValue.push_back(stat.mean(polyValues[iband]));
1660 fieldname.push_back(fs.str());
1661 break;
1662 case(rule::median):
1663 theValue.push_back(stat.median(polyValues[iband]));
1664 fieldname.push_back(fs.str());
1665 break;
1666 case(rule::stdev):
1667 theValue.push_back(sqrt(stat.var(polyValues[iband])));
1668 fieldname.push_back(fs.str());
1669 break;
1670 case(rule::percentile):{
1671 for(int iperc=0;iperc<percentile_opt.size();++iperc){
1672 theValue.push_back(stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[iperc]));
1673 ostringstream fsperc;
1674 fsperc << fs.str() << percentile_opt[iperc];
1675 fieldname.push_back(fsperc.str());
1676 }
1677 break;
1678 }
1679 case(rule::sum):
1680 theValue.push_back(stat.sum(polyValues[iband]));
1681 fieldname.push_back(fs.str());
1682 break;
1683 case(rule::max):
1684 theValue.push_back(stat.mymax(polyValues[iband]));
1685 fieldname.push_back(fs.str());
1686 break;
1687 case(rule::min):
1688 theValue.push_back(stat.mymin(polyValues[iband]));
1689 fieldname.push_back(fs.str());
1690 break;
1691 case(rule::centroid):
1692 theValue.push_back(polyValues[iband].back());
1693 fieldname.push_back(fs.str());
1694 break;
1695 default://not supported
1696 break;
1697 }
1698 for(int ivalue=0;ivalue<theValue.size();++ivalue){
1699 switch( fieldType ){
1700 case OFTInteger:
1701 writePolygonFeature->SetField(fieldname[ivalue].c_str(),static_cast<int>(theValue[ivalue]));
1702 break;
1703 case OFTReal:
1704 writePolygonFeature->SetField(fieldname[ivalue].c_str(),theValue[ivalue]);
1705 break;
1706 case OFTString:
1707 writePolygonFeature->SetField(fieldname[ivalue].c_str(),type2string<double>(theValue[ivalue]).c_str());
1708 break;
1709 default://not supported
1710 std::string errorString="field type not supported";
1711 throw(errorString);
1712 break;
1713 }
1714 }
1715 }
1716 }
1717 }
1718 }
1719 if(createPolygon&&validFeature){
1720 //todo: only create if valid feature?
1721 //write polygon feature
1722 if(verbose_opt[0]>1)
1723 std::cout << "creating polygon feature (2)" << std::endl;
1724 if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1725 std::string errorString="Failed to create polygon feature in ogr vector dataset";
1726 throw(errorString);
1727 }
1728 //test: no need to destroy anymore?
1729 // OGRFeature::DestroyFeature( writePolygonFeature );
1730 ++ntotalvalid;
1731 ++ntotalvalidLayer;
1732 }
1733 }
1734 ++ifeature;
1735 if(theThreshold>0){
1736 if(threshold_opt.size()==layer_opt.size())
1737 progress=(100.0/theThreshold)*static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
1738 else
1739 progress=static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
1740 }
1741 else
1742 progress=static_cast<float>(ifeature+1)/(-theThreshold);
1743 pfnProgress(progress,pszMessage,pProgressArg);
1744 }
1745 catch(std::string e){
1746 std::cout << e << std::endl;
1747 continue;
1748 }
1749 catch(int npoint){
1750 if(verbose_opt[0])
1751 std::cout << "number of points read in polygon: " << npoint << std::endl;
1752 continue;
1753 }
1754 }
1755 // if(rbox_opt[0]>0||cbox_opt[0]>0)
1756 // boxWriter.close();
1757 progress=1.0;
1758 pfnProgress(progress,pszMessage,pProgressArg);
1759 ++ilayerWrite;
1760 }//for ilayer
1761 sampleReaderOgr.close();
1762 ogrWriter.close();
1763 }
1764 progress=1.0;
1765 pfnProgress(progress,pszMessage,pProgressArg);
1766 imgReader.close();
1767}
double getDeltaY(void) const
Get the pixel cell spacing in y.
int nrOfRow(void) const
Get the number of rows of this dataset.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
double getUly() const
Get the upper left corner y (georeferenced) coordinate of this dataset.
double getLrx() const
Get the lower right corner x (georeferenced) coordinate of this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
double getDeltaX(void) const
Get the pixel cell spacing in x.
double getLry() const
Get the lower right corner y (georeferenced) coordinate of this dataset.
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked,...
double getUlx() const
Get the upper left corner x (georeferenced) coordinate of this dataset.
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double &minValue, double &maxValue)
Get the minimum and maximum cell values for a specific band in a region of interest defined by startC...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...