Processing Kernel for remote sensing data
pkextract.cc
1 /**********************************************************************
2 pkextract.cc: extract pixel values from raster image from a (vector or raster) sample
3 Copyright (C) 2008-2014 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along 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 namespace rule{
38  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};
39 }
40 
41 using namespace std;
42 
43 int main(int argc, char *argv[])
44 {
45  Optionpk<string> image_opt("i", "input", "Raster input dataset containing band information");
46  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.");
47  Optionpk<string> layer_opt("ln", "ln", "Layer name(s) in sample (leave empty to select all)");
48  Optionpk<unsigned int> random_opt("rand", "random", "Create simple random sample of points. Provide number of points to generate");
49  Optionpk<double> grid_opt("grid", "grid", "Create systematic grid of points. Provide cell grid size (in projected units, e.g,. m)");
50  Optionpk<string> output_opt("o", "output", "Output sample dataset");
51  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 or proportion");
52  Optionpk<float> threshold_opt("t", "threshold", "Probability threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold for vector sample datasets. 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);
53  Optionpk<string> ogrformat_opt("f", "f", "Output sample dataset format","SQLite");
54  Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or Integer)", "Real");
55  Optionpk<string> ltype_opt("lt", "ltype", "Label type: In16 or String", "Integer");
56  Optionpk<bool> polygon_opt("polygon", "polygon", "Create OGRPolygon as geometry instead of OGRPoint.", false);
57  Optionpk<int> band_opt("b", "band", "Band index(es) to extract. Leave empty to use all bands");
58  Optionpk<string> rule_opt("r", "rule", "Rule how to report image information per feature (only for vector sample). point (value at each point or at centroid if polygon), centroid, mean, stdev, median, proportion, min, max, mode, sum.", "centroid");
59  Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "Invalid value(s) for input image");
60  Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Band(s) in input image to check if pixel is valid (used for srcnodata)", 0);
61  Optionpk<float> polythreshold_opt("tp", "thresholdPolygon", "(absolute) threshold for selecting samples in each polygon");
62  Optionpk<string> test_opt("test", "test", "Test sample dataset (use this option in combination with threshold<100 to create a training (output) and test set");
63  Optionpk<string> fieldname_opt("bn", "bname", "For single band input data, this extra attribute name will correspond to the raster values. For multi-band input data, multiple attributes with this prefix will be added (e.g. b0, b1, b2, etc.)", "b");
64  Optionpk<string> label_opt("cn", "cname", "Name of the class label in the output vector dataset", "label");
65  Optionpk<short> geo_opt("geo", "geo", "Use geo coordinates (set to 0 to use image coordinates)", 1);
66  Optionpk<short> down_opt("down", "down", "Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
67  Optionpk<short> buffer_opt("buf", "buffer", "Buffer for calculating statistics for point features ");
68  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);
69  Optionpk<short> verbose_opt("v", "verbose", "Verbose mode if > 0", 0);
70 
71  bool doProcess;//stop process when program was invoked with help option (-h --help)
72  try{
73  doProcess=image_opt.retrieveOption(argc,argv);
74  sample_opt.retrieveOption(argc,argv);
75  layer_opt.retrieveOption(argc,argv);
76  random_opt.retrieveOption(argc,argv);
77  grid_opt.retrieveOption(argc,argv);
78  output_opt.retrieveOption(argc,argv);
79  class_opt.retrieveOption(argc,argv);
80  threshold_opt.retrieveOption(argc,argv);
81  ogrformat_opt.retrieveOption(argc,argv);
82  ftype_opt.retrieveOption(argc,argv);
83  ltype_opt.retrieveOption(argc,argv);
84  polygon_opt.retrieveOption(argc,argv);
85  band_opt.retrieveOption(argc,argv);
86  rule_opt.retrieveOption(argc,argv);
87  bndnodata_opt.retrieveOption(argc,argv);
88  srcnodata_opt.retrieveOption(argc,argv);
89  polythreshold_opt.retrieveOption(argc,argv);
90  test_opt.retrieveOption(argc,argv);
91  fieldname_opt.retrieveOption(argc,argv);
92  label_opt.retrieveOption(argc,argv);
93  geo_opt.retrieveOption(argc,argv);
94  down_opt.retrieveOption(argc,argv);
95  buffer_opt.retrieveOption(argc,argv);
96  disc_opt.retrieveOption(argc,argv);
97  // rbox_opt.retrieveOption(argc,argv);
98  // cbox_opt.retrieveOption(argc,argv);
99  verbose_opt.retrieveOption(argc,argv);
100  }
101  catch(string predefinedString){
102  std::cout << predefinedString << std::endl;
103  exit(0);
104  }
105  if(!doProcess){
106  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
107  exit(0);//help was invoked, stop processing
108  }
109 
110  std::map<std::string, rule::RULE_TYPE> ruleMap;
111  //initialize ruleMap
112  ruleMap["point"]=rule::point;
113  ruleMap["centroid"]=rule::centroid;
114  ruleMap["mean"]=rule::mean;
115  ruleMap["stdev"]=rule::stdev;
116  ruleMap["median"]=rule::median;
117  ruleMap["proportion"]=rule::proportion;
118  ruleMap["min"]=rule::min;
119  ruleMap["max"]=rule::max;
120  ruleMap["custom"]=rule::custom;
121  ruleMap["mode"]=rule::mode;
122  ruleMap["sum"]=rule::sum;
123 
124  if(srcnodata_opt.size()){
125  while(srcnodata_opt.size()<bndnodata_opt.size())
126  srcnodata_opt.push_back(srcnodata_opt[0]);
127  while(bndnodata_opt.size()<srcnodata_opt.size())
128  bndnodata_opt.push_back(bndnodata_opt[0]);
129  }
130 
131  if(verbose_opt[0])
132  std::cout << class_opt << std::endl;
133  statfactory::StatFactory stat;
134  Vector2d<unsigned int> posdata;
135  unsigned long int nsample=0;
136  unsigned long int ntotalvalid=0;
137  unsigned long int ntotalinvalid=0;
138  vector<unsigned long int> nvalid(class_opt.size());
139  vector<unsigned long int> ninvalid(class_opt.size());
140  for(int it=0;it<nvalid.size();++it){
141  nvalid[it]=0;
142  ninvalid[it]=0;
143  }
144 
145  ImgReaderGdal imgReader;
146  if(image_opt.empty()){
147  std::cerr << "No image dataset provided (use option -i). Use --help for help information";
148  exit(0);
149  }
150  if(output_opt.empty()){
151  std::cerr << "No output dataset provided (use option -o). Use --help for help information";
152  exit(0);
153  }
154  try{
155  imgReader.open(image_opt[0]);
156  }
157  catch(std::string errorstring){
158  std::cout << errorstring << std::endl;
159  exit(0);
160  }
161 
162  int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand();
163 
164  if(fieldname_opt.size()<nband){
165  std::string bandString=fieldname_opt[0];
166  fieldname_opt.clear();
167  fieldname_opt.resize(nband);
168  for(int iband=0;iband<nband;++iband){
169  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
170  ostringstream fs;
171  fs << bandString << theBand;
172  fieldname_opt[iband]=fs.str();
173  }
174  }
175 
176  if(verbose_opt[0])
177  std::cout << fieldname_opt << std::endl;
178 
179  if(verbose_opt[0]>1)
180  std::cout << "Number of bands in input image: " << imgReader.nrOfBand() << std::endl;
181 
182  OGRFieldType fieldType;
183  OGRFieldType labelType;
184  int ogr_typecount=11;//hard coded for now!
185  if(verbose_opt[0]>1)
186  std::cout << "field and label types can be: ";
187  for(int iType = 0; iType < ogr_typecount; ++iType){
188  if(verbose_opt[0]>1)
189  std::cout << " " << OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType);
190  if( OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType) != NULL
191  && EQUAL(OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType),
192  ftype_opt[0].c_str()))
193  fieldType=(OGRFieldType) iType;
194  if( OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType) != NULL
195  && EQUAL(OGRFieldDefn::GetFieldTypeName((OGRFieldType)iType),
196  ltype_opt[0].c_str()))
197  labelType=(OGRFieldType) iType;
198  }
199  switch( fieldType ){
200  case OFTInteger:
201  case OFTReal:
202  case OFTRealList:
203  case OFTString:
204  if(verbose_opt[0]>1)
205  std::cout << std::endl << "field type is: " << OGRFieldDefn::GetFieldTypeName(fieldType) << std::endl;
206  break;
207  default:
208  cerr << "field type " << OGRFieldDefn::GetFieldTypeName(fieldType) << " not supported" << std::endl;
209  exit(0);
210  break;
211  }
212  switch( labelType ){
213  case OFTInteger:
214  case OFTReal:
215  case OFTRealList:
216  case OFTString:
217  if(verbose_opt[0]>1)
218  std::cout << std::endl << "label type is: " << OGRFieldDefn::GetFieldTypeName(labelType) << std::endl;
219  break;
220  default:
221  cerr << "label type " << OGRFieldDefn::GetFieldTypeName(labelType) << " not supported" << std::endl;
222  exit(0);
223  break;
224  }
225 
226  const char* pszMessage;
227  void* pProgressArg=NULL;
228  GDALProgressFunc pfnProgress=GDALTermProgress;
229  double progress=0;
230  srand(time(NULL));
231 
232  bool sampleIsRaster=false;
233  bool sampleIsVirtual=false;
234 
235  ImgReaderOgr sampleReaderOgr;
236  ImgWriterOgr sampleWriterOgr;
237 
238  if(sample_opt.size()){
239  try{
240  sampleReaderOgr.open(sample_opt[0]);
241  }
242  catch(string errorString){
243  sampleIsRaster=true;
244  }
245  }
246  else{
247  sampleWriterOgr.open("/vsimem/sample",ogrformat_opt[0]);
248  char **papszOptions=NULL;
249  sampleWriterOgr.createLayer("points", imgReader.getProjection(), wkbPoint, papszOptions);
250  sampleIsVirtual=true;
251 
252  // string fieldName="label";
253  // string fieldValue="class";
254  // sampleWriterOgr.createField(fieldName,OFTString);
255  if(random_opt.size()){
256  //create simple random sampling within boundary
257  OGRPoint pt;
258  double ulx,uly,lrx,lry;
259  imgReader.getBoundingBox(ulx,uly,lrx,lry);
260  for(unsigned int ipoint=1;ipoint<=random_opt[0];++ipoint){
261  OGRFeature *pointFeature;
262  pointFeature=sampleWriterOgr.createFeature();
263  // pointFeature->SetField(fieldName.c_str(),fieldValue.c_str());
264  double theX=ulx+static_cast<double>(rand())/(RAND_MAX)*(lrx-ulx);
265  double theY=uly-static_cast<double>(rand())/(RAND_MAX)*(uly-lry);
266  pt.setX(theX);
267  pt.setY(theY);
268  pointFeature->SetGeometry( &pt );
269  if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
270  cerr << "Failed to create feature in shapefile" << endl;
271  exit( 1 );
272  }
273  OGRFeature::DestroyFeature(pointFeature);
274  }
275  }
276  else if(grid_opt.size()){
277  //create systematic grid of points
278  OGRPoint pt;
279  double ulx,uly,lrx,lry;
280  imgReader.getBoundingBox(ulx,uly,lrx,lry);
281  unsigned int ipoint=0;
282  for(double theY=uly-grid_opt[0]/2;theY>lry;theY-=grid_opt[0]){
283  for(double theX=ulx+grid_opt[0]/2;theX<lrx;theX+=grid_opt[0]){
284  if(verbose_opt[0]>1)
285  cout << "position: " << theX << " " << theY << endl;
286  OGRFeature *pointFeature;
287  pointFeature=sampleWriterOgr.createFeature();
288  // pointFeature->SetField(fieldName.c_str(),fieldValue.c_str());
289  pt.setX(theX);
290  pt.setY(theY);
291  pointFeature->SetGeometry( &pt );
292  if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
293  cerr << "Failed to create feature in shapefile" << endl;
294  exit( 1 );
295  }
296  OGRFeature::DestroyFeature(pointFeature);
297  }
298  }
299  }
300  else{
301  std::cerr << "No sample dataset provided (use option -s). Use --help for help information";
302  exit(0);
303  }
304  sampleReaderOgr.open("/vsimem/sample");
305  }
306 
307  if(sampleIsRaster){
308  if(class_opt.empty()){
309  // std::cout << "Warning: no classes selected, if a classes must be extracted, set to -1 for all classes using option -c -1" << std::endl;
310  ImgReaderGdal classReader;
311  ImgWriterOgr ogrWriter;
312  assert(sample_opt.size());
313  classReader.open(sample_opt[0]);
314  // vector<int> classBuffer(classReader.nrOfCol());
315  vector<double> classBuffer(classReader.nrOfCol());
316  Vector2d<double> imgBuffer(nband);//[band][col]
317  vector<double> sample(2+nband);//x,y,band values
318  Vector2d<double> writeBuffer;
319  // vector<int> writeBufferClass;
320  vector<double> writeBufferClass;
321  vector<int> selectedClass;
322  Vector2d<double> selectedBuffer;
323  double oldimgrow=-1;
324  int irow=0;
325  int icol=0;
326  if(verbose_opt[0]>1)
327  std::cout << "extracting sample from image..." << std::endl;
328  progress=0;
329  pfnProgress(progress,pszMessage,pProgressArg);
330  for(irow=0;irow<classReader.nrOfRow();++irow){
331  if(irow%down_opt[0])
332  continue;
333  // classReader.readData(classBuffer,GDT_Int32,irow);
334  classReader.readData(classBuffer,GDT_Float64,irow);
335  double x,y;//geo coordinates
336  double iimg,jimg;//image coordinates in img image
337  for(icol=0;icol<classReader.nrOfCol();++icol){
338  if(icol%down_opt[0])
339  continue;
340  // int theClass=0;
341  double theClass=classBuffer[icol];
342  // int processClass=-1;
343  int processClass=0;
344  // if(class_opt[0]<0){//process every class except 0
345  // if(classBuffer[icol]){
346  // processClass=0;
347  // theClass=classBuffer[icol];
348  // }
349  // }
350  // else{
351  // for(int iclass=0;iclass<class_opt.size();++iclass){
352  // if(classBuffer[icol]==class_opt[iclass]){
353  // processClass=iclass;
354  // theClass=class_opt[iclass];
355  // }
356  // }
357  // }
358  // if(processClass>=0){
359  bool valid=true;
360  if(valid){
361  if(geo_opt[0]){
362  classReader.image2geo(icol,irow,x,y);
363  sample[0]=x;
364  sample[1]=y;
365  if(verbose_opt[0]>1){
366  std::cout.precision(12);
367  std::cout << theClass << " " << x << " " << y << std::endl;
368  }
369  //find col in img
370  imgReader.geo2image(x,y,iimg,jimg);
371  //nearest neighbour
372  jimg=static_cast<int>(jimg);
373  iimg=static_cast<int>(iimg);
374  if(static_cast<int>(iimg)<0||static_cast<int>(iimg)>=imgReader.nrOfCol())
375  continue;
376  }
377  else{
378  iimg=icol;
379  jimg=irow;
380  sample[0]=iimg;
381  sample[1]=jimg;
382  }
383  if(static_cast<int>(jimg)<0||static_cast<int>(jimg)>=imgReader.nrOfRow())
384  continue;
385 
386  bool valid=true;
387 
388  if(static_cast<int>(jimg)!=static_cast<int>(oldimgrow)){
389  assert(imgBuffer.size()==nband);
390  for(int iband=0;iband<nband;++iband){
391  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
392  imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
393  assert(imgBuffer[iband].size()==imgReader.nrOfCol());
394  if(srcnodata_opt.size()){
395  vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
396  if(bndit!=bndnodata_opt.end()){
397  vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
398  if(bndit!=bndnodata_opt.end()){
399  if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
400  valid=false;
401  }
402  }
403  }
404  }
405  oldimgrow=jimg;
406  }
407 
408  if(valid){
409  for(int iband=0;iband<imgBuffer.size();++iband){
410  if(imgBuffer[iband].size()!=imgReader.nrOfCol()){
411  std::cout << "Error in band " << iband << ": " << imgBuffer[iband].size() << "!=" << imgReader.nrOfCol() << std::endl;
412  assert(imgBuffer[iband].size()==imgReader.nrOfCol());
413  }
414  sample[iband+2]=imgBuffer[iband][static_cast<int>(iimg)];
415  }
416  float theThreshold=(threshold_opt.size()>1)?threshold_opt[processClass]:threshold_opt[0];
417  if(theThreshold>0){//percentual value
418  double p=static_cast<double>(rand())/(RAND_MAX);
419  p*=100.0;
420  if(p>theThreshold)
421  continue;//do not select for now, go to next column
422  }
423  else if(nvalid.size()>processClass){//absolute value
424  if(nvalid[processClass]>-theThreshold)
425  continue;//do not select any more pixels for this class, go to next column to search for other classes
426  }
427  writeBuffer.push_back(sample);
428  writeBufferClass.push_back(theClass);
429  ++ntotalvalid;
430  if(nvalid.size()>processClass)
431  ++(nvalid[processClass]);
432  }
433  else{
434  ++ntotalinvalid;
435  if(ninvalid.size()>processClass)
436  ++(ninvalid[processClass]);
437  }
438  }//processClass
439  }//icol
440  progress=static_cast<float>(irow+1.0)/classReader.nrOfRow();
441  pfnProgress(progress,pszMessage,pProgressArg);
442  }//irow
443  progress=100;
444  pfnProgress(progress,pszMessage,pProgressArg);
445  if(writeBuffer.size()>0){
446  assert(ntotalvalid==writeBuffer.size());
447  if(verbose_opt[0]>0)
448  std::cout << "creating image sample writer " << output_opt[0] << " with " << writeBuffer.size() << " samples (" << ntotalinvalid << " invalid)" << std::endl;
449  ogrWriter.open(output_opt[0],ogrformat_opt[0]);
450  char **papszOptions=NULL;
451  ostringstream slayer;
452  slayer << "training data";
453  std::string layername=slayer.str();
454  ogrWriter.createLayer(layername, imgReader.getProjection(), wkbPoint, papszOptions);
455  std::string fieldname="fid";//number of the point
456  ogrWriter.createField(fieldname,OFTInteger);
457  map<std::string,double> pointAttributes;
458  ogrWriter.createField(label_opt[0],labelType);
459  for(int iband=0;iband<nband;++iband){
460  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
461  ogrWriter.createField(fieldname_opt[iband],fieldType);
462  }
463  std::cout << "writing sample to " << output_opt[0] << "..." << std::endl;
464  progress=0;
465  pfnProgress(progress,pszMessage,pProgressArg);
466  for(int isample=0;isample<writeBuffer.size();++isample){
467  if(verbose_opt[0]>1)
468  std::cout << "writing sample " << isample << std::endl;
469  pointAttributes[label_opt[0]]=writeBufferClass[isample];
470  for(int iband=0;iband<writeBuffer[0].size()-2;++iband){
471  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
472  // ostringstream fs;
473  // if(nband==1)
474  // fs << fieldname_opt[0];
475  // else
476  // fs << fieldname_opt[0] << theBand;
477  // pointAttributes[fs.str()]=writeBuffer[isample][iband+2];
478  pointAttributes[fieldname_opt[iband]]=writeBuffer[isample][iband+2];
479  }
480  if(verbose_opt[0]>1)
481  std::cout << "all bands written" << std::endl;
482  ogrWriter.addPoint(writeBuffer[isample][0],writeBuffer[isample][1],pointAttributes,fieldname,isample);
483  progress=static_cast<float>(isample+1.0)/writeBuffer.size();
484  pfnProgress(progress,pszMessage,pProgressArg);
485  }
486  ogrWriter.close();
487  }
488  else{
489  std::cout << "No data found for any class " << std::endl;
490  }
491  classReader.close();
492  nsample=writeBuffer.size();
493  if(verbose_opt[0])
494  std::cout << "total number of samples written: " << nsample << std::endl;
495  }
496  else{//class_opt.size()!=0
497  assert(class_opt[0]);
498  // if(class_opt[0]){
499  assert(threshold_opt.size()==1||threshold_opt.size()==class_opt.size());
500  ImgReaderGdal classReader;
501  ImgWriterOgr ogrWriter;
502  if(verbose_opt[0]>1){
503  std::cout << "reading position from sample dataset " << std::endl;
504  std::cout << "class thresholds: " << std::endl;
505  for(int iclass=0;iclass<class_opt.size();++iclass){
506  if(threshold_opt.size()>1)
507  std::cout << class_opt[iclass] << ": " << threshold_opt[iclass] << std::endl;
508  else
509  std::cout << class_opt[iclass] << ": " << threshold_opt[0] << std::endl;
510  }
511  }
512  classReader.open(sample_opt[0]);
513  vector<int> classBuffer(classReader.nrOfCol());
514  // vector<double> classBuffer(classReader.nrOfCol());
515  Vector2d<double> imgBuffer(nband);//[band][col]
516  vector<double> sample(2+nband);//x,y,band values
517  Vector2d<double> writeBuffer;
518  vector<int> writeBufferClass;
519  // vector<double> writeBufferClass;
520  vector<int> selectedClass;
521  Vector2d<double> selectedBuffer;
522  double oldimgrow=-1;
523  int irow=0;
524  int icol=0;
525  if(verbose_opt[0]>1)
526  std::cout << "extracting sample from image..." << std::endl;
527  progress=0;
528  pfnProgress(progress,pszMessage,pProgressArg);
529  for(irow=0;irow<classReader.nrOfRow();++irow){
530  if(irow%down_opt[0])
531  continue;
532  classReader.readData(classBuffer,GDT_Int32,irow);
533  // classReader.readData(classBuffer,GDT_Float64,irow);
534  double x,y;//geo coordinates
535  double iimg,jimg;//image coordinates in img image
536  for(icol=0;icol<classReader.nrOfCol();++icol){
537  if(icol%down_opt[0])
538  continue;
539  int theClass=0;
540  // double theClass=0;
541  int processClass=-1;
542  if(class_opt.empty()){//process every class
543  if(classBuffer[icol]){
544  processClass=0;
545  theClass=classBuffer[icol];
546  }
547  }
548  else{
549  for(int iclass=0;iclass<class_opt.size();++iclass){
550  if(classBuffer[icol]==class_opt[iclass]){
551  processClass=iclass;
552  theClass=class_opt[iclass];
553  }
554  }
555  }
556  if(processClass>=0){
557  // if(classBuffer[icol]==class_opt[0]){
558  if(geo_opt[0]){
559  classReader.image2geo(icol,irow,x,y);
560  sample[0]=x;
561  sample[1]=y;
562  if(verbose_opt[0]>1){
563  std::cout.precision(12);
564  std::cout << theClass << " " << x << " " << y << std::endl;
565  }
566  //find col in img
567  imgReader.geo2image(x,y,iimg,jimg);
568  //nearest neighbour
569  jimg=static_cast<int>(jimg);
570  iimg=static_cast<int>(iimg);
571  if(static_cast<int>(iimg)<0||static_cast<int>(iimg)>=imgReader.nrOfCol())
572  continue;
573  }
574  else{
575  iimg=icol;
576  jimg=irow;
577  sample[0]=iimg;
578  sample[1]=jimg;
579  }
580  if(static_cast<int>(jimg)<0||static_cast<int>(jimg)>=imgReader.nrOfRow())
581  continue;
582 
583  bool valid=true;
584 
585  if(static_cast<int>(jimg)!=static_cast<int>(oldimgrow)){
586  assert(imgBuffer.size()==nband);
587  for(int iband=0;iband<nband;++iband){
588  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
589  imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
590  assert(imgBuffer[iband].size()==imgReader.nrOfCol());
591  if(srcnodata_opt.size()){
592  vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
593  if(bndit!=bndnodata_opt.end()){
594  if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
595  valid=false;
596  }
597  }
598  }
599  oldimgrow=jimg;
600  }
601  if(valid){
602  for(int iband=0;iband<imgBuffer.size();++iband){
603  if(imgBuffer[iband].size()!=imgReader.nrOfCol()){
604  std::cout << "Error in band " << iband << ": " << imgBuffer[iband].size() << "!=" << imgReader.nrOfCol() << std::endl;
605  assert(imgBuffer[iband].size()==imgReader.nrOfCol());
606  }
607  sample[iband+2]=imgBuffer[iband][static_cast<int>(iimg)];
608  }
609  float theThreshold=(threshold_opt.size()>1)?threshold_opt[processClass]:threshold_opt[0];
610  if(theThreshold>0){//percentual value
611  double p=static_cast<double>(rand())/(RAND_MAX);
612  p*=100.0;
613  if(p>theThreshold)
614  continue;//do not select for now, go to next column
615  }
616  else if(nvalid.size()>processClass){//absolute value
617  if(nvalid[processClass]>-theThreshold)
618  continue;//do not select any more pixels for this class, go to next column to search for other classes
619  }
620  writeBuffer.push_back(sample);
621  // writeBufferClass.push_back(class_opt[processClass]);
622  writeBufferClass.push_back(theClass);
623  ++ntotalvalid;
624  if(nvalid.size()>processClass)
625  ++(nvalid[processClass]);
626  }
627  else{
628  ++ntotalinvalid;
629  if(ninvalid.size()>processClass)
630  ++(ninvalid[processClass]);
631  }
632  }//processClass
633  }//icol
634  progress=static_cast<float>(irow+1.0)/classReader.nrOfRow();
635  pfnProgress(progress,pszMessage,pProgressArg);
636  }//irow
637  if(writeBuffer.size()>0){
638  assert(ntotalvalid==writeBuffer.size());
639  if(verbose_opt[0]>0)
640  std::cout << "creating image sample writer " << output_opt[0] << " with " << writeBuffer.size() << " samples (" << ntotalinvalid << " invalid)" << std::endl;
641  ogrWriter.open(output_opt[0],ogrformat_opt[0]);
642  char **papszOptions=NULL;
643  ostringstream slayer;
644  slayer << "training data";
645  std::string layername=slayer.str();
646  ogrWriter.createLayer(layername, imgReader.getProjection(), wkbPoint, papszOptions);
647  std::string fieldname="fid";//number of the point
648  ogrWriter.createField(fieldname,OFTInteger);
649  map<std::string,double> pointAttributes;
650  // ogrWriter.createField(label_opt[0],OFTInteger);
651  ogrWriter.createField(label_opt[0],labelType);
652  for(int iband=0;iband<nband;++iband){
653  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
654  // ostringstream fs;
655  // if(nband==1)
656  // fs << fieldname_opt[0];
657  // else
658  // fs << fieldname_opt[0] << theBand;
659  // ogrWriter.createField(fs.str(),fieldType);
660  ogrWriter.createField(fieldname_opt[iband],fieldType);
661  }
662  pfnProgress(progress,pszMessage,pProgressArg);
663  std::cout << "writing sample to " << output_opt[0] << "..." << std::endl;
664  progress=0;
665  pfnProgress(progress,pszMessage,pProgressArg);
666  for(int isample=0;isample<writeBuffer.size();++isample){
667  pointAttributes[label_opt[0]]=writeBufferClass[isample];
668  for(int iband=0;iband<writeBuffer[0].size()-2;++iband){
669  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
670  // ostringstream fs;
671  // if(nband==1)
672  // fs << fieldname_opt[0];
673  // else
674  // fs << fieldname_opt[0] << theBand;
675  // pointAttributes[fs.str()]=writeBuffer[isample][iband+2];
676  pointAttributes[fieldname_opt[iband]]=writeBuffer[isample][iband+2];
677  }
678  ogrWriter.addPoint(writeBuffer[isample][0],writeBuffer[isample][1],pointAttributes,fieldname,isample);
679  progress=static_cast<float>(isample+1.0)/writeBuffer.size();
680  pfnProgress(progress,pszMessage,pProgressArg);
681  }
682  ogrWriter.close();
683  }
684  else{
685  std::cout << "No data found for any class " << std::endl;
686  }
687  classReader.close();
688  nsample=writeBuffer.size();
689  if(verbose_opt[0]){
690  std::cout << "total number of samples written: " << nsample << std::endl;
691  if(nvalid.size()==class_opt.size()){
692  for(int iclass=0;iclass<class_opt.size();++iclass)
693  std::cout << "class " << class_opt[iclass] << " has " << nvalid[iclass] << " samples" << std::endl;
694  }
695  }
696  }
697  }
698  else{//vector dataset
699  if(verbose_opt[0]>1)
700  std::cout << "creating image sample writer " << output_opt[0] << std::endl;
701  ImgWriterOgr ogrWriter;
702  ImgWriterOgr ogrTestWriter;
703  ogrWriter.open(output_opt[0],ogrformat_opt[0]);
704  if(test_opt.size()){
705  if(verbose_opt[0]>1)
706  std::cout << "creating image test writer " << test_opt[0] << std::endl;
707  ogrTestWriter.open(test_opt[0],ogrformat_opt[0]);
708  }
709 
710  //support multiple layers
711  int nlayerRead=sampleReaderOgr.getDataSource()->GetLayerCount();
712  int ilayerWrite=0;
713  if(verbose_opt[0])
714  std::cout << "number of layers: " << nlayerRead << endl;
715 
716  for(int ilayer=0;ilayer<nlayerRead;++ilayer){
717  OGRLayer *readLayer=sampleReaderOgr.getLayer(ilayer);
718  string currentLayername=readLayer->GetName();
719  if(layer_opt.size())
720  if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
721  continue;
722  cout << "processing layer " << currentLayername << endl;
723 
724  readLayer->ResetReading();
725  OGRLayer *writeLayer;
726  OGRLayer *writeTestLayer;
727 
728  if(polygon_opt[0]){
729  if(verbose_opt[0])
730  std::cout << "create polygons" << std::endl;
731  char **papszOptions=NULL;
732  writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPolygon, papszOptions);
733  if(test_opt.size())
734  writeTestLayer=ogrTestWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPolygon, papszOptions);
735  }
736  else{
737  if(verbose_opt[0])
738  std::cout << "create points in layer " << readLayer->GetName() << std::endl;
739  char **papszOptions=NULL;
740 
741  writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPoint, papszOptions);
742  if(test_opt.size()){
743  char **papszOptions=NULL;
744  writeTestLayer=ogrTestWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPoint, papszOptions);
745  }
746  }
747  if(verbose_opt[0])
748  std::cout << "copy fields from layer " << ilayer << std::flush << std::endl;
749  ogrWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
750 
751  if(test_opt.size()){
752  if(verbose_opt[0])
753  std::cout << "copy fields test writer" << std::flush << std::endl;
754  ogrTestWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
755  }
756  // vector<std::string> fieldnames;
757  // if(verbose_opt[0])
758  // std::cout << "get fields" << std::flush << std::endl;
759  // sampleReaderOgr.getFields(fieldnames);
760  // assert(fieldnames.size()==ogrWriter.getFieldCount(ilayerWrite));
761  // map<std::string,double> pointAttributes;
762 
763  if(class_opt.size()){
764  switch(ruleMap[rule_opt[0]]){
765  case(rule::proportion):{//proportion for each class
766  for(int iclass=0;iclass<class_opt.size();++iclass){
767  ostringstream cs;
768  cs << class_opt[iclass];
769  ogrWriter.createField(cs.str(),fieldType,ilayerWrite);
770  }
771  break;
772  }
773  case(rule::custom):
774  case(rule::mode):
775  ogrWriter.createField(label_opt[0],fieldType,ilayerWrite);
776  if(test_opt.size())
777  ogrTestWriter.createField(label_opt[0],fieldType,ilayerWrite);
778  break;
779  }
780  }
781  else{
782  for(int iband=0;iband<nband;++iband){
783  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
784  ostringstream fs;
785  fs << fieldname_opt[iband];
786  if(verbose_opt[0]>1)
787  std::cout << "creating field " << fs.str() << std::endl;
788 
789  ogrWriter.createField(fs.str(),fieldType,ilayerWrite);
790  if(test_opt.size())
791  ogrTestWriter.createField(fs.str(),fieldType,ilayerWrite);
792  }
793  }
794  OGRFeature *readFeature;
795  unsigned long int ifeature=0;
796  unsigned long int nfeature=sampleReaderOgr.getFeatureCount();
797  progress=0;
798  pfnProgress(progress,pszMessage,pProgressArg);
799  while( (readFeature = readLayer->GetNextFeature()) != NULL ){
800  bool validFeature=false;
801  bool writeTest=false;//write this feature to test_opt[0] instead of output_opt
802  if(verbose_opt[0]>0)
803  std::cout << "reading feature " << readFeature->GetFID() << std::endl;
804  if(threshold_opt[0]>0){//percentual value
805  double p=static_cast<double>(rand())/(RAND_MAX);
806  p*=100.0;
807  if(p>threshold_opt[0]){
808  if(test_opt.size())
809  writeTest=true;
810  else
811  continue;//do not select for now, go to next feature
812  }
813  }
814  else{//absolute value
815  if(ntotalvalid>-threshold_opt[0]){
816  if(test_opt.size())
817  writeTest=true;
818  else
819  continue;//do not select any more pixels, go to next column feature
820  }
821  }
822  if(verbose_opt[0]>0)
823  std::cout << "processing feature " << readFeature->GetFID() << std::endl;
824  //get x and y from readFeature
825  double x,y;
826  OGRGeometry *poGeometry;
827  poGeometry = readFeature->GetGeometryRef();
828  assert(poGeometry!=NULL);
829  try{
830  if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ){
831 
832  if(!buffer_opt.size()){
833  switch(ruleMap[rule_opt[0]]){
834  case(rule::point):
835  case(rule::centroid):
836  buffer_opt.push_back(1);//default
837  break;
838  default:
839  buffer_opt.push_back(3);//default
840  }
841  }
842 
843  if(verbose_opt[0]>1)
844  std::cout << "boundary: " << buffer_opt[0] << std::endl;
845 
846  OGRPolygon writePolygon;
847  OGRLinearRing writeRing;
848  OGRPoint writeCentroidPoint;
849  OGRFeature *writePolygonFeature;
850  OGRFeature *writeCentroidFeature;
851 
852  OGRPoint *poPoint = (OGRPoint *) poGeometry;
853  writeCentroidPoint=*poPoint;
854 
855  x=poPoint->getX();
856  y=poPoint->getY();
857 
858  double i_centre,j_centre;
859  if(geo_opt[0])
860  imgReader.geo2image(x,y,i_centre,j_centre);
861  else{
862  i_centre=x;
863  j_centre=y;
864  }
865  //nearest neighbour
866  j_centre=static_cast<int>(j_centre);
867  i_centre=static_cast<int>(i_centre);
868 
869  double uli=i_centre-buffer_opt[0]/2;
870  double ulj=j_centre-buffer_opt[0]/2;
871  double lri=i_centre+buffer_opt[0]/2;
872  double lrj=j_centre+buffer_opt[0]/2;
873  // double uli=i_centre-buffer_opt[0]/2-0.5;
874  // double ulj=j_centre-buffer_opt[0]/2-0.5;
875  // double lri=i_centre+buffer_opt[0]/2+0.5;
876  // double lrj=j_centre+buffer_opt[0]/2+0.5;
877 
878  //nearest neighbour
879  ulj=static_cast<int>(ulj);
880  uli=static_cast<int>(uli);
881  lrj=static_cast<int>(lrj);
882  lri=static_cast<int>(lri);
883 
884  if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
885  uli=i_centre;
886  ulj=j_centre;
887  lri=i_centre;
888  lrj=j_centre;
889  }
890 
891  //check if j is out of bounds
892  if(static_cast<int>(ulj)<0||static_cast<int>(ulj)>=imgReader.nrOfRow())
893  continue;
894  //check if j is out of bounds
895  if(static_cast<int>(uli)<0||static_cast<int>(lri)>=imgReader.nrOfCol())
896  continue;
897 
898  OGRPoint ulPoint,urPoint,llPoint,lrPoint;
899  double ulx,uly;
900  double urx,ury;
901 
902  if(polygon_opt[0]){
903  if(disc_opt[0]){
904  double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
905  double radius=buffer_opt[0]/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
906  unsigned short nstep = 25;
907  for(int i=0;i<nstep;++i){
908  OGRPoint aPoint;
909  aPoint.setX(x+radius*cos(2*PI*i/nstep));
910  aPoint.setY(y+radius*sin(2*PI*i/nstep));
911  writeRing.addPoint(&aPoint);
912  }
913  writePolygon.addRing(&writeRing);
914  writePolygon.closeRings();
915  }
916  else{
917  double llx,lly;
918  double lrx,lry;
919  imgReader.image2geo(uli,ulj,ulx,uly);
920  imgReader.image2geo(lri,lrj,lrx,lry);
921  ulPoint.setX(ulx-imgReader.getDeltaX()/2.0);
922  ulPoint.setY(uly+imgReader.getDeltaY()/2.0);
923  lrPoint.setX(lrx+imgReader.getDeltaX()/2.0);
924  lrPoint.setY(lry-imgReader.getDeltaY()/2.0);
925  urPoint.setX(lrx+imgReader.getDeltaX()/2.0);
926  urPoint.setY(uly+imgReader.getDeltaY()/2.0);
927  llPoint.setX(ulx-imgReader.getDeltaX()/2.0);
928  llPoint.setY(lry-imgReader.getDeltaY()/2.0);
929 
930  writeRing.addPoint(&ulPoint);
931  writeRing.addPoint(&urPoint);
932  writeRing.addPoint(&lrPoint);
933  writeRing.addPoint(&llPoint);
934  writePolygon.addRing(&writeRing);
935  writePolygon.closeRings();
936  }
937  }
938 
939  int nPointWindow=0;//similar to nPointPolygon for polygons
940  if(polygon_opt[0]){
941  if(writeTest)
942  writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
943  else
944  writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
945  }
946  else if(ruleMap[rule_opt[0]]!=rule::point){
947  if(writeTest)
948  writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
949  else
950  writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
951  }
952  Vector2d<double> windowValues;
953  vector<double> windowClassValues;
954 
955  if(class_opt.size()){
956  windowClassValues.resize(class_opt.size());
957  //initialize
958  for(int iclass=0;iclass<class_opt.size();++iclass)
959  windowClassValues[iclass]=0;
960  }
961  else
962  windowValues.resize(nband);
963  vector< Vector2d<double> > readValues(nband);
964  for(int iband=0;iband<nband;++iband){
965  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
966  //test
967  assert(uli>=0);
968  assert(uli<imgReader.nrOfCol());
969  assert(lri>=0);
970  assert(lri<imgReader.nrOfCol());
971  assert(ulj>=0);
972  assert(ulj<imgReader.nrOfRow());
973  assert(lrj>=0);
974  assert(lrj<imgReader.nrOfRow());
975  imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
976  }
977 
978  OGRPoint thePoint;
979  for(int j=ulj;j<=lrj;++j){
980  for(int i=uli;i<=lri;++i){
981  //check if within raster image
982  if(i<0||i>=imgReader.nrOfCol())
983  continue;
984  if(j<0||j>=imgReader.nrOfRow())
985  continue;
986  //no need to check if point is on surface
987  double theX=0;
988  double theY=0;
989  imgReader.image2geo(i,j,theX,theY);
990  thePoint.setX(theX);
991  thePoint.setY(theY);
992 
993  if(disc_opt[0]&&buffer_opt[0]>1){
994  double radius=buffer_opt[0]/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
995  if((theX-x)*(theX-x)+(theY-y)*(theY-y)>radius*radius)
996  continue;
997  }
998  bool valid=true;
999 
1000  if(srcnodata_opt.size()){
1001  for(int vband=0;vband<bndnodata_opt.size();++vband){
1002  double value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
1003  if(value==srcnodata_opt[vband]){
1004  valid=false;
1005  break;
1006  }
1007  }
1008  }
1009 
1010  if(!valid)
1011  continue;
1012  else
1013  validFeature=true;
1014 
1015  // writeRing.addPoint(&thePoint);//already done
1016 
1017  ++nPointWindow;
1018  OGRFeature *writePointFeature;
1019  if(!polygon_opt[0]){
1020  //create feature
1021  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, stdev, median, sum or centroid (only create point at centroid)
1022  if(writeTest)
1023  writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1024  else
1025  writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1026  if(verbose_opt[0]>1)
1027  std::cout << "copying fields from polygons " << std::endl;
1028  if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
1029  cerr << "writing feature failed" << std::endl;
1030  writePointFeature->SetGeometry(&thePoint);
1031  OGRGeometry *updateGeometry;
1032  updateGeometry = writePointFeature->GetGeometryRef();
1033  OGRPoint *poPoint = (OGRPoint *) updateGeometry;
1034  if(verbose_opt[0]>1)
1035  std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
1036  }
1037  }
1038  if(class_opt.size()){
1039  short value=((readValues[0])[j-ulj])[i-uli];
1040  for(int iclass=0;iclass<class_opt.size();++iclass){
1041  if(value==class_opt[iclass])
1042  windowClassValues[iclass]+=1;
1043  }
1044  }
1045  else{
1046  for(int iband=0;iband<nband;++iband){
1047  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1048  assert(j-ulj>=0);
1049  assert(j-ulj<readValues[iband].size());
1050  assert(i-uli>=0);
1051  assert(i-uli<((readValues[iband])[j-ulj]).size());
1052  double value=((readValues[iband])[j-ulj])[i-uli];
1053  // imgReader.readData(value,GDT_Float64,i,j,theBand);
1054  if(verbose_opt[0]>1)
1055  std::cout << ": " << value << std::endl;
1056  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
1057  windowValues[iband].push_back(value);
1058  }
1059  else{
1060  if(verbose_opt[0]>1)
1061  std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
1062  switch( fieldType ){
1063  case OFTInteger:
1064  case OFTReal:
1065  writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
1066  break;
1067  case OFTString:
1068  writePointFeature->SetField(fieldname_opt[iband].c_str(),type2string<double>(value).c_str());
1069  break;
1070  default://not supported
1071  assert(0);
1072  break;
1073  }
1074  }//else
1075  }//iband
1076  }//else (class_opt.size())
1077  if(!polygon_opt[0]){
1078  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean or median value (only at centroid)
1079  //write feature
1080  if(verbose_opt[0]>1)
1081  std::cout << "creating point feature" << std::endl;
1082  if(writeTest){
1083  if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1084  std::string errorString="Failed to create feature in test ogr vector dataset";
1085  throw(errorString);
1086  }
1087  }
1088  else{
1089  if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1090  std::string errorString="Failed to create feature in ogr vector dataset";
1091  throw(errorString);
1092  }
1093  }
1094  //destroy feature
1095  OGRFeature::DestroyFeature( writePointFeature );
1096  ++ntotalvalid;
1097  if(verbose_opt[0])
1098  std::cout << "ntotalvalid(2): " << ntotalvalid << std::endl;
1099  }
1100  }
1101  }
1102  }
1103  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
1104  //do not create if no points found within polygon
1105  if(!nPointWindow){
1106  if(verbose_opt[0])
1107  cout << "no points found in window, continuing" << endl;
1108  continue;
1109  }
1110  //add ring to polygon
1111  if(polygon_opt[0]){
1112  // writePolygon.addRing(&writeRing);//already done
1113  // writePolygon.closeRings();//already done
1114  //write geometry of writePolygon
1115  if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
1116  cerr << "writing feature failed" << std::endl;
1117  writePolygonFeature->SetGeometry(&writePolygon);
1118  if(verbose_opt[0]>1)
1119  std::cout << "copying new fields write polygon " << std::endl;
1120  if(verbose_opt[0]>1)
1121  std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
1122  //write polygon feature
1123  }
1124  else{//write value of polygon to centroid point
1125  //create feature
1126  if(verbose_opt[0]>1)
1127  std::cout << "copying fields from polygons " << std::endl;
1128  if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
1129  cerr << "writing feature failed" << std::endl;
1130  writeCentroidFeature->SetGeometry(&writeCentroidPoint);
1131  OGRGeometry *updateGeometry;
1132  updateGeometry = writeCentroidFeature->GetGeometryRef();
1133  assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
1134  if(verbose_opt[0]>1)
1135  std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
1136  }
1137  if(class_opt.empty()){
1138  if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
1139  if(verbose_opt[0])
1140  std::cout << "number of points in window: " << nPointWindow << std::endl;
1141  for(int index=0;index<windowValues.size();++index){
1142  //test
1143  if(windowValues[index].size()!=1){
1144  cerr << "Error: windowValues[index].size()=" << windowValues[index].size() << endl;
1145  assert(windowValues[index].size()==1);
1146  }
1147  double theValue=windowValues[index].back();
1148 
1149  if(verbose_opt[0])
1150  std::cout << "number of points in window: " << nPointWindow << std::endl;
1151  int theBand=(band_opt.size()) ? band_opt[index] : index;
1152 
1153  if(verbose_opt[0]>1)
1154  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
1155  switch( fieldType ){
1156  case OFTInteger:
1157  case OFTReal:
1158  if(polygon_opt[0])
1159  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
1160  else
1161  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
1162  break;
1163  case OFTString:
1164  if(polygon_opt[0])
1165  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1166  else
1167  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1168  break;
1169  default://not supported
1170  std::cout << "field type not supported yet..." << std::endl;
1171  break;
1172  }
1173  }
1174  }
1175  else{//ruleMap[rule_opt[0]] is not rule::point
1176  double theValue=0;
1177  for(int index=0;index<windowValues.size();++index){
1178  if(windowValues[index].size()!=buffer_opt[0]*buffer_opt[0]){
1179  cerr << "Error: windowValues[index].size()=" << windowValues[index].size() << endl;
1180  }
1181  if(ruleMap[rule_opt[0]]==rule::mean)
1182  theValue=stat.mean(windowValues[index]);
1183  else if(ruleMap[rule_opt[0]]==rule::stdev)
1184  theValue=sqrt(stat.var(windowValues[index]));
1185  else if(ruleMap[rule_opt[0]]==rule::median)
1186  theValue=stat.median(windowValues[index]);
1187  else if(ruleMap[rule_opt[0]]==rule::sum)
1188  theValue=stat.sum(windowValues[index]);
1189  else if(ruleMap[rule_opt[0]]==rule::max)
1190  theValue=stat.mymax(windowValues[index]);
1191  else if(ruleMap[rule_opt[0]]==rule::min)
1192  theValue=stat.mymin(windowValues[index]);
1193  else{//rule::centroid
1194  if(verbose_opt[0])
1195  std::cout << "number of points in polygon: " << nPointWindow << std::endl;
1196  assert(nPointWindow<=1);
1197  assert(nPointWindow==windowValues[index].size());
1198  theValue=windowValues[index].back();
1199  }
1200  if(verbose_opt[0]>1)
1201  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
1202  switch( fieldType ){
1203  case OFTInteger:
1204  case OFTReal:
1205  if(polygon_opt[0])
1206  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
1207  else
1208  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
1209  break;
1210  case OFTString:
1211  if(polygon_opt[0])
1212  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1213  else
1214  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1215  break;
1216  default://not supported
1217  std::cout << "field type not supported yet..." << std::endl;
1218  break;
1219  }
1220  }
1221  }
1222  }
1223  else{//class_opt is set
1224  if(ruleMap[rule_opt[0]]==rule::proportion){
1225  if(verbose_opt[0])
1226  std::cout << "number of points in polygon: " << nPointWindow << std::endl;
1227  stat.normalize_pct(windowClassValues);
1228  for(int index=0;index<windowClassValues.size();++index){
1229  double theValue=windowClassValues[index];
1230  ostringstream fs;
1231  fs << class_opt[index];
1232  if(polygon_opt[0])
1233  writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
1234  else
1235  writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
1236  }
1237  }
1238  else if(ruleMap[rule_opt[0]]==rule::custom){
1239  assert(polygon_opt[0]);//not implemented for points
1240  if(verbose_opt[0])
1241  std::cout << "number of points in polygon: " << nPointWindow << std::endl;
1242  stat.normalize_pct(windowClassValues);
1243  assert(windowClassValues.size()==2);//11:broadleaved, 12:coniferous
1244  if(windowClassValues[0]>=75)//broadleaved
1245  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
1246  else if(windowClassValues[1]>=75)//coniferous
1247  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
1248  else if(windowClassValues[0]>25&&windowClassValues[1]>25)//mixed
1249  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
1250  else{
1251  if(verbose_opt[0]){
1252  std::cout << "No valid value in windowClassValues..." << std::endl;
1253  for(int index=0;index<windowClassValues.size();++index){
1254  double theValue=windowClassValues[index];
1255  std::cout << theValue << " ";
1256  }
1257  std::cout << std::endl;
1258  }
1259  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
1260  }
1261  }
1262  else if(ruleMap[rule_opt[0]]==rule::mode){
1263  //maximum votes in polygon
1264  if(verbose_opt[0])
1265  std::cout << "number of points in window: " << nPointWindow << std::endl;
1266  //search for class with maximum votes
1267  int maxClass=stat.mymin(class_opt);
1268  vector<double>::iterator maxit;
1269  maxit=stat.mymax(windowClassValues,windowClassValues.begin(),windowClassValues.end());
1270  int maxIndex=distance(windowClassValues.begin(),maxit);
1271  maxClass=class_opt[maxIndex];
1272  if(verbose_opt[0]>0)
1273  std::cout << "maxClass: " << maxClass << std::endl;
1274  if(polygon_opt[0])
1275  writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
1276  else
1277  writeCentroidFeature->SetField(label_opt[0].c_str(),maxClass);
1278  }
1279  }
1280  if(polygon_opt[0]){
1281  if(verbose_opt[0]>1)
1282  std::cout << "creating polygon feature" << std::endl;
1283  if(writeTest){
1284  if(writeTestLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1285  std::string errorString="Failed to create polygon feature in ogr vector dataset";
1286  throw(errorString);
1287  }
1288  }
1289  else{
1290  if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1291  std::string errorString="Failed to create polygon feature in ogr vector dataset";
1292  throw(errorString);
1293  }
1294  }
1295  OGRFeature::DestroyFeature( writePolygonFeature );
1296  ++ntotalvalid;
1297  if(verbose_opt[0])
1298  std::cout << "ntotalvalid(1): " << ntotalvalid << std::endl;
1299  }
1300  else{
1301  if(verbose_opt[0]>1)
1302  std::cout << "creating point feature in centroid" << std::endl;
1303  if(writeTest){
1304  if(writeTestLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
1305  std::string errorString="Failed to create point feature in ogr vector dataset";
1306  throw(errorString);
1307  }
1308  }
1309  else{
1310  //test
1311  assert(validFeature);
1312  if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
1313  std::string errorString="Failed to create point feature in ogr vector dataset";
1314  throw(errorString);
1315  }
1316  }
1317  OGRFeature::DestroyFeature( writeCentroidFeature );
1318  ++ntotalvalid;
1319  if(verbose_opt[0])
1320  std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
1321  }
1322  }
1323  }//if wkbPoint
1324  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
1325 
1326  OGRPolygon readPolygon = *((OGRPolygon *) poGeometry);
1327  OGRPolygon writePolygon;
1328  OGRLinearRing writeRing;
1329  OGRPoint writeCentroidPoint;
1330  OGRFeature *writePolygonFeature;
1331  OGRFeature *writeCentroidFeature;
1332 
1333  readPolygon.closeRings();
1334 
1335  if(verbose_opt[0]>1)
1336  std::cout << "get point on polygon" << std::endl;
1337  if(ruleMap[rule_opt[0]]==rule::centroid)
1338  readPolygon.Centroid(&writeCentroidPoint);
1339  else
1340  readPolygon.PointOnSurface(&writeCentroidPoint);
1341 
1342  double ulx,uly,lrx,lry;
1343  double uli,ulj,lri,lrj;
1344  if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
1345  ulx=writeCentroidPoint.getX();
1346  uly=writeCentroidPoint.getY();
1347  lrx=ulx;
1348  lry=uly;
1349  }
1350  else{
1351  //get envelope
1352  if(verbose_opt[0])
1353  std::cout << "reading envelope for polygon " << ifeature << std::endl;
1354  OGREnvelope* psEnvelope=new OGREnvelope();
1355  readPolygon.getEnvelope(psEnvelope);
1356  ulx=psEnvelope->MinX;
1357  uly=psEnvelope->MaxY;
1358  lrx=psEnvelope->MaxX;
1359  lry=psEnvelope->MinY;
1360  delete psEnvelope;
1361  }
1362  if(geo_opt[0]){
1363  imgReader.geo2image(ulx,uly,uli,ulj);
1364  imgReader.geo2image(lrx,lry,lri,lrj);
1365  }
1366  else{
1367  uli=ulx;
1368  ulj=uly;
1369  lri=lrx;
1370  lrj=lry;
1371  }
1372  //nearest neighbour
1373  ulj=static_cast<int>(ulj);
1374  uli=static_cast<int>(uli);
1375  lrj=static_cast<int>(lrj);
1376  lri=static_cast<int>(lri);
1377  //iterate through all pixels
1378  if(verbose_opt[0]>1)
1379  std::cout << "bounding box for polygon feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
1380 
1381  if(uli<0)
1382  uli=0;
1383  if(lri<0)
1384  lri=0;
1385  if(uli>=imgReader.nrOfCol())
1386  uli=imgReader.nrOfCol()-1;
1387  if(lri>=imgReader.nrOfCol())
1388  lri=imgReader.nrOfCol()-1;
1389  if(ulj<0)
1390  ulj=0;
1391  if(lrj<0)
1392  lrj=0;
1393  if(ulj>=imgReader.nrOfRow())
1394  ulj=imgReader.nrOfRow()-1;
1395  if(lrj>=imgReader.nrOfRow())
1396  lrj=imgReader.nrOfRow()-1;
1397  // if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||lrj>=imgReader.nrOfRow())
1398  // continue;
1399 
1400  int nPointPolygon=0;
1401 
1402  if(polygon_opt[0]){
1403  if(writeTest)
1404  writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1405  else
1406  writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1407  }
1408  else if(ruleMap[rule_opt[0]]!=rule::point){
1409  if(writeTest)
1410  writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1411  else
1412  writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1413  }
1414  // vector<double> polyValues;
1415  Vector2d<double> polyValues;
1416  vector<double> polyClassValues;
1417 
1418  if(class_opt.size()){
1419 
1420  polyClassValues.resize(class_opt.size());
1421  //initialize
1422  for(int iclass=0;iclass<class_opt.size();++iclass)
1423  polyClassValues[iclass]=0;
1424  }
1425  else
1426  polyValues.resize(nband);
1427  vector< Vector2d<double> > readValues(nband);
1428  for(int iband=0;iband<nband;++iband){
1429  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1430  //test
1431  assert(uli>=0);
1432  assert(uli<imgReader.nrOfCol());
1433  assert(lri>=0);
1434  assert(lri<imgReader.nrOfCol());
1435  assert(ulj>=0);
1436  assert(ulj<imgReader.nrOfRow());
1437  assert(lrj>=0);
1438  assert(lrj<imgReader.nrOfRow());
1439  imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
1440  }
1441 
1442  OGRPoint thePoint;
1443  for(int j=ulj;j<=lrj;++j){
1444  for(int i=uli;i<=lri;++i){
1445  //check if within raster image
1446  if(i<0||i>=imgReader.nrOfCol())
1447  continue;
1448  if(j<0||j>=imgReader.nrOfRow())
1449  continue;
1450  //check if point is on surface
1451  double theX=0;
1452  double theY=0;
1453  imgReader.image2geo(i,j,theX,theY);
1454  thePoint.setX(theX);
1455  thePoint.setY(theY);
1456 
1457  if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
1458  continue;
1459 
1460  bool valid=true;
1461 
1462  if(srcnodata_opt.size()){
1463  for(int vband=0;vband<bndnodata_opt.size();++vband){
1464  double value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
1465  if(value==srcnodata_opt[vband]){
1466  valid=false;
1467  break;
1468  }
1469  }
1470  }
1471 
1472  if(!valid)
1473  continue;
1474  else
1475  validFeature=true;
1476 
1477  writeRing.addPoint(&thePoint);//todo: check if I need to add all interior points to ring or do I need to check if point is on ring first?
1478  if(verbose_opt[0]>1)
1479  std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
1480  ++nPointPolygon;
1481 
1482  if(polythreshold_opt.size())
1483  if(nPointPolygon>polythreshold_opt[0])
1484  continue;
1485  // throw(nPointPolygon);
1486  OGRFeature *writePointFeature;
1487  if(!polygon_opt[0]){
1488  //create feature
1489  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, stdev, median, sum or centroid (only create point at centroid)
1490  if(writeTest)
1491  writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1492  else
1493  writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1494  if(verbose_opt[0]>1)
1495  std::cout << "copying fields from polygons " << std::endl;
1496  if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
1497  cerr << "writing feature failed" << std::endl;
1498  writePointFeature->SetGeometry(&thePoint);
1499  OGRGeometry *updateGeometry;
1500  updateGeometry = writePointFeature->GetGeometryRef();
1501  OGRPoint *poPoint = (OGRPoint *) updateGeometry;
1502  if(verbose_opt[0]>1)
1503  std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
1504  }
1505  }
1506  if(class_opt.size()){
1507  short value=((readValues[0])[j-ulj])[i-uli];
1508  for(int iclass=0;iclass<class_opt.size();++iclass){
1509  if(value==class_opt[iclass])
1510  polyClassValues[iclass]+=1;
1511  }
1512  }
1513  else{
1514  for(int iband=0;iband<nband;++iband){
1515  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1516  assert(j-ulj>=0);
1517  assert(j-ulj<readValues[iband].size());
1518  assert(i-uli>=0);
1519  assert(i-uli<((readValues[iband])[j-ulj]).size());
1520  double value=((readValues[iband])[j-ulj])[i-uli];
1521  // imgReader.readData(value,GDT_Float64,i,j,theBand);
1522  if(verbose_opt[0]>1)
1523  std::cout << ": " << value << std::endl;
1524  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point)
1525  polyValues[iband].push_back(value);
1526  else{
1527  if(verbose_opt[0]>1)
1528  std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
1529  switch( fieldType ){
1530  case OFTInteger:
1531  case OFTReal:
1532  writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
1533  break;
1534  case OFTString:
1535  writePointFeature->SetField(fieldname_opt[iband].c_str(),type2string<double>(value).c_str());
1536  break;
1537  default://not supported
1538  assert(0);
1539  break;
1540  }
1541  }//else
1542  }//iband
1543  }//else (class_opt.size())
1544  if(!polygon_opt[0]){
1545  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean or median value (only at centroid)
1546  //write feature
1547  if(verbose_opt[0]>1)
1548  std::cout << "creating point feature" << std::endl;
1549  if(writeTest){
1550  if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1551  std::string errorString="Failed to create feature in test ogr vector dataset";
1552  throw(errorString);
1553  }
1554  }
1555  else{
1556  if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
1557  std::string errorString="Failed to create feature in ogr vector dataset";
1558  throw(errorString);
1559  }
1560  }
1561  //destroy feature
1562  OGRFeature::DestroyFeature( writePointFeature );
1563  ++ntotalvalid;
1564  if(verbose_opt[0])
1565  std::cout << "ntotalvalid(2): " << ntotalvalid << std::endl;
1566  }
1567  }
1568  }
1569  }
1570  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
1571  //do not create if no points found within polygon
1572  if(!nPointPolygon){
1573  if(verbose_opt[0])
1574  cout << "no points found in polygon, continuing" << endl;
1575  continue;
1576  }
1577  //add ring to polygon
1578  if(polygon_opt[0]){
1579  writePolygon.addRing(&writeRing);
1580  writePolygon.closeRings();
1581  //write geometry of writePolygon
1582  if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
1583  cerr << "writing feature failed" << std::endl;
1584  writePolygonFeature->SetGeometry(&writePolygon);
1585  if(verbose_opt[0]>1)
1586  std::cout << "copying new fields write polygon " << std::endl;
1587  if(verbose_opt[0]>1)
1588  std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
1589  //write polygon feature
1590  }
1591  else{//write value of polygon to centroid point
1592  //create feature
1593  if(verbose_opt[0]>1)
1594  std::cout << "copying fields from polygons " << std::endl;
1595  if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
1596  cerr << "writing feature failed" << std::endl;
1597  writeCentroidFeature->SetGeometry(&writeCentroidPoint);
1598  OGRGeometry *updateGeometry;
1599  updateGeometry = writeCentroidFeature->GetGeometryRef();
1600  assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
1601  if(verbose_opt[0]>1)
1602  std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
1603  }
1604  if(class_opt.empty()){
1605  if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
1606  if(verbose_opt[0])
1607  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1608  for(int index=0;index<polyValues.size();++index){
1609  //test
1610  assert(polyValues[index].size()==1);
1611  double theValue=polyValues[index].back();
1612 
1613  if(verbose_opt[0])
1614  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1615  int theBand=(band_opt.size()) ? band_opt[index] : index;
1616 
1617  if(verbose_opt[0]>1)
1618  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
1619  switch( fieldType ){
1620  case OFTInteger:
1621  case OFTReal:
1622  if(polygon_opt[0])
1623  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
1624  else
1625  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
1626  break;
1627  case OFTString:
1628  if(polygon_opt[0])
1629  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1630  else
1631  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1632  break;
1633  default://not supported
1634  std::cout << "field type not supported yet..." << std::endl;
1635  break;
1636  }
1637  }
1638  }
1639  else{//ruleMap[rule_opt[0]] is not rule::point
1640  double theValue=0;
1641  for(int index=0;index<polyValues.size();++index){
1642  if(ruleMap[rule_opt[0]]==rule::mean)
1643  theValue=stat.mean(polyValues[index]);
1644  else if(ruleMap[rule_opt[0]]==rule::stdev)
1645  theValue=sqrt(stat.var(polyValues[index]));
1646  else if(ruleMap[rule_opt[0]]==rule::median)
1647  theValue=stat.median(polyValues[index]);
1648  else if(ruleMap[rule_opt[0]]==rule::sum)
1649  theValue=stat.sum(polyValues[index]);
1650  else if(ruleMap[rule_opt[0]]==rule::max)
1651  theValue=stat.mymax(polyValues[index]);
1652  else if(ruleMap[rule_opt[0]]==rule::min)
1653  theValue=stat.mymin(polyValues[index]);
1654  else{//rule::centroid
1655  if(verbose_opt[0])
1656  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1657  assert(nPointPolygon<=1);
1658  assert(nPointPolygon==polyValues[index].size());
1659  theValue=polyValues[index].back();
1660  }
1661  if(verbose_opt[0]>1)
1662  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
1663  switch( fieldType ){
1664  case OFTInteger:
1665  case OFTReal:
1666  if(polygon_opt[0])
1667  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
1668  else
1669  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
1670  break;
1671  case OFTString:
1672  if(polygon_opt[0])
1673  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1674  else
1675  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
1676  break;
1677  default://not supported
1678  std::cout << "field type not supported yet..." << std::endl;
1679  break;
1680  }
1681  }
1682  }
1683  }
1684  else{//class_opt is set
1685  if(ruleMap[rule_opt[0]]==rule::proportion){
1686  if(verbose_opt[0])
1687  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1688  stat.normalize_pct(polyClassValues);
1689  for(int index=0;index<polyClassValues.size();++index){
1690  double theValue=polyClassValues[index];
1691  ostringstream fs;
1692  fs << class_opt[index];
1693  if(polygon_opt[0])
1694  writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
1695  else
1696  writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
1697  }
1698  }
1699  else if(ruleMap[rule_opt[0]]==rule::custom){
1700  assert(polygon_opt[0]);//not implemented for points
1701  if(verbose_opt[0])
1702  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1703  stat.normalize_pct(polyClassValues);
1704  assert(polyClassValues.size()==2);//11:broadleaved, 12:coniferous
1705  if(polyClassValues[0]>=75)//broadleaved
1706  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
1707  else if(polyClassValues[1]>=75)//coniferous
1708  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
1709  else if(polyClassValues[0]>25&&polyClassValues[1]>25)//mixed
1710  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
1711  else{
1712  if(verbose_opt[0]){
1713  std::cout << "No valid value in polyClassValues..." << std::endl;
1714  for(int index=0;index<polyClassValues.size();++index){
1715  double theValue=polyClassValues[index];
1716  std::cout << theValue << " ";
1717  }
1718  std::cout << std::endl;
1719  }
1720  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
1721  }
1722  }
1723  else if(ruleMap[rule_opt[0]]==rule::mode){
1724  //maximum votes in polygon
1725  if(verbose_opt[0])
1726  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
1727  //search for class with maximum votes
1728  int maxClass=stat.mymin(class_opt);
1729  vector<double>::iterator maxit;
1730  maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
1731  int maxIndex=distance(polyClassValues.begin(),maxit);
1732  maxClass=class_opt[maxIndex];
1733  if(verbose_opt[0]>0)
1734  std::cout << "maxClass: " << maxClass << std::endl;
1735  if(polygon_opt[0])
1736  writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
1737  else
1738  writeCentroidFeature->SetField(label_opt[0].c_str(),maxClass);
1739  }
1740  }
1741  if(polygon_opt[0]){
1742  if(verbose_opt[0]>1)
1743  std::cout << "creating polygon feature" << std::endl;
1744  if(writeTest){
1745  if(writeTestLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1746  std::string errorString="Failed to create polygon feature in ogr vector dataset";
1747  throw(errorString);
1748  }
1749  }
1750  else{
1751  if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
1752  std::string errorString="Failed to create polygon feature in ogr vector dataset";
1753  throw(errorString);
1754  }
1755  }
1756  OGRFeature::DestroyFeature( writePolygonFeature );
1757  ++ntotalvalid;
1758  if(verbose_opt[0])
1759  std::cout << "ntotalvalid(1): " << ntotalvalid << std::endl;
1760  }
1761  else{
1762  if(verbose_opt[0]>1)
1763  std::cout << "creating point feature in centroid" << std::endl;
1764  if(writeTest){
1765  if(writeTestLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
1766  std::string errorString="Failed to create point feature in ogr vector dataset";
1767  throw(errorString);
1768  }
1769  }
1770  else{
1771  //test
1772  assert(validFeature);
1773  if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
1774  std::string errorString="Failed to create point feature in ogr vector dataset";
1775  throw(errorString);
1776  }
1777  }
1778  OGRFeature::DestroyFeature( writeCentroidFeature );
1779  ++ntotalvalid;
1780  if(verbose_opt[0])
1781  std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
1782  }
1783  }
1784  }
1785  else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){//todo: try to use virtual OGRGeometry instead of OGRMultiPolygon and OGRPolygon
1786  OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
1787  OGRPolygon writePolygon;
1788  OGRLinearRing writeRing;
1789  OGRPoint writeCentroidPoint;
1790  OGRFeature *writePolygonFeature;
1791  OGRFeature *writeCentroidFeature;
1792 
1793  readPolygon.closeRings();
1794 
1795  if(verbose_opt[0]>1)
1796  std::cout << "get centroid point from polygon" << std::endl;
1797 
1798  readPolygon.Centroid(&writeCentroidPoint);
1799 
1800  double ulx,uly,lrx,lry;
1801  double uli,ulj,lri,lrj;
1802  if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
1803  ulx=writeCentroidPoint.getX();
1804  uly=writeCentroidPoint.getY();
1805  lrx=ulx;
1806  lry=uly;
1807  }
1808  else{
1809  //get envelope
1810  if(verbose_opt[0])
1811  std::cout << "reading envelope for polygon " << ifeature << std::endl;
1812  OGREnvelope* psEnvelope=new OGREnvelope();
1813  readPolygon.getEnvelope(psEnvelope);
1814  ulx=psEnvelope->MinX;
1815  uly=psEnvelope->MaxY;
1816  lrx=psEnvelope->MaxX;
1817  lry=psEnvelope->MinY;
1818  delete psEnvelope;
1819  }
1820  // if(geo_opt[0]){
1821  imgReader.geo2image(ulx,uly,uli,ulj);
1822  imgReader.geo2image(lrx,lry,lri,lrj);
1823  // }
1824  // else{
1825  // uli=ulx;
1826  // ulj=uly;
1827  // lri=lrx;
1828  // lrj=lry;
1829  // }
1830  //nearest neighbour
1831  ulj=static_cast<int>(ulj);
1832  uli=static_cast<int>(uli);
1833  lrj=static_cast<int>(lrj);
1834  lri=static_cast<int>(lri);
1835  //iterate through all pixels
1836  if(verbose_opt[0]>1)
1837  std::cout << "bounding box for multipologon feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
1838 
1839  if(uli<0)
1840  uli=0;
1841  if(lri<0)
1842  lri=0;
1843  if(uli>=imgReader.nrOfCol())
1844  uli=imgReader.nrOfCol()-1;
1845  if(lri>=imgReader.nrOfCol())
1846  lri=imgReader.nrOfCol()-1;
1847  if(ulj<0)
1848  ulj=0;
1849  if(lrj<0)
1850  lrj=0;
1851  if(ulj>=imgReader.nrOfRow())
1852  ulj=imgReader.nrOfRow()-1;
1853  if(lrj>=imgReader.nrOfRow())
1854  lrj=imgReader.nrOfRow()-1;
1855  // if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||lrj>=imgReader.nrOfRow())
1856  // continue;
1857 
1858  int nPointPolygon=0;
1859  if(polygon_opt[0]){
1860  if(writeTest)
1861  writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1862  else
1863  writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1864  }
1865  else if(ruleMap[rule_opt[0]]!=rule::point){
1866  if(writeTest)
1867  writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1868  else
1869  writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1870  }
1871  // vector<double> polyValues;
1872  Vector2d<double> polyValues;
1873  vector<double> polyClassValues;
1874 
1875  if(class_opt.size()){
1876  polyClassValues.resize(class_opt.size());
1877  //initialize
1878  for(int iclass=0;iclass<class_opt.size();++iclass)
1879  polyClassValues[iclass]=0;
1880  }
1881  else
1882  polyValues.resize(nband);
1883  vector< Vector2d<double> > readValues(nband);
1884  for(int iband=0;iband<nband;++iband){
1885  int theBand=(band_opt.size()) ? band_opt[iband] : iband;
1886  //test
1887  assert(uli>=0);
1888  assert(uli<imgReader.nrOfCol());
1889  assert(lri>=0);
1890  assert(lri<imgReader.nrOfCol());
1891  assert(ulj>=0);
1892  assert(ulj<imgReader.nrOfRow());
1893  assert(lrj>=0);
1894  assert(lrj<imgReader.nrOfRow());
1895  imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
1896  }
1897 
1898  OGRPoint thePoint;
1899  for(int j=ulj;j<=lrj;++j){
1900  for(int i=uli;i<=lri;++i){
1901  //check if within raster image
1902  if(i<0||i>=imgReader.nrOfCol())
1903  continue;
1904  if(j<0||j>=imgReader.nrOfRow())
1905  continue;
1906  //check if point is on surface
1907  double theX=0;
1908  double theY=0;
1909  imgReader.image2geo(i,j,theX,theY);
1910  thePoint.setX(theX);
1911  thePoint.setY(theY);
1912 
1913  if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
1914  continue;
1915 
1916  bool valid=true;
1917 
1918  if(srcnodata_opt.size()){
1919  for(int vband=0;vband<bndnodata_opt.size();++vband){
1920  double value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
1921  if(value==srcnodata_opt[vband]){
1922  valid=false;
1923  break;
1924  }
1925  }
1926  }
1927 
1928  if(!valid)
1929  continue;
1930  else
1931  validFeature=true;
1932 
1933  writeRing.addPoint(&thePoint);
1934  if(verbose_opt[0]>1)
1935  std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
1936  ++nPointPolygon;
1937 
1938  if(polythreshold_opt.size())
1939  if(nPointPolygon>polythreshold_opt[0])
1940  continue;
1941  // throw(nPointPolygon);
1942  OGRFeature *writePointFeature;
1943  if(!polygon_opt[0]){
1944  //create feature
1945  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, stdev, median or sum (only create point at centroid)
1946  if(writeTest)
1947  writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
1948  else
1949  writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
1950  if(verbose_opt[0]>1)
1951  std::cout << "copying fields from polygons " << std::endl;
1952  if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
1953  cerr << "writing feature failed" << std::endl;
1954  writePointFeature->SetGeometry(&thePoint);
1955  OGRGeometry *updateGeometry;
1956  updateGeometry = writePointFeature->GetGeometryRef();
1957  OGRPoint *poPoint = (OGRPoint *) updateGeometry;
1958  if(verbose_opt[0]>1)
1959  std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
1960  }
1961  }
1962  if(class_opt.size()){
1963  short value=((readValues[0])[j-ulj])[i-uli];
1964  for(int iclass=0;iclass<class_opt.size();++iclass){
1965  if(value==class_opt[iclass])
1966  polyClassValues[iclass]+=1;
1967  }
1968  }
1969  else{
1970  for(int iband=0;iband<nband;++iband){
1971  //test
1972  assert(j-ulj>=0);
1973  assert(j-ulj<readValues[iband].size());
1974  assert(i-uli>=0);
1975  assert(i-uli<((readValues[iband])[j-ulj]).size());
1976  double value=((readValues[iband])[j-ulj])[i-uli];
1977  // imgReader.readData(value,GDT_Float64,i,j,theBand);
1978  if(verbose_opt[0]>1)
1979  std::cout << ": " << value << std::endl;
1980  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point)
1981  polyValues[iband].push_back(value);
1982  else{
1983  if(verbose_opt[0]>1)
1984  std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
1985  switch( fieldType ){
1986  case OFTInteger:
1987  case OFTReal:
1988  writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
1989  break;
1990  case OFTString:
1991  writePointFeature->SetField(fieldname_opt[iband].c_str(),type2string<double>(value).c_str());
1992  break;
1993  // case OFTRealList:{
1994  // int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
1995  // int nCount;
1996  // const double *theList;
1997  // theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
1998  // vector<double> vectorList(nCount+1);
1999  // for(int index=0;index<nCount;++index)
2000  // vectorList[nCount]=value;
2001  // writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
2002  // break;
2003  // }
2004  default://not supported
2005  assert(0);
2006  break;
2007  }
2008  }//else
2009  }//iband
2010  }//else (class_opt.size())
2011  if(!polygon_opt[0]){
2012  if(ruleMap[rule_opt[0]]==rule::point){//do not create in case of mean, stdev or median value (only at centroid)
2013  //write feature
2014  if(verbose_opt[0]>1)
2015  std::cout << "creating point feature" << std::endl;
2016  if(writeTest){
2017  if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
2018  std::string errorString="Failed to create feature in ogr vector dataset";
2019  throw(errorString);
2020  }
2021  }
2022  else{
2023  if(writeLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
2024  std::string errorString="Failed to create feature in ogr vector dataset";
2025  throw(errorString);
2026  }
2027  }
2028  //destroy feature
2029  OGRFeature::DestroyFeature( writePointFeature );
2030  }
2031  }
2032  // ++isample;
2033  ++ntotalvalid;
2034  if(verbose_opt[0])
2035  std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
2036  }
2037  }
2038 
2039  //test
2040  if(!validFeature)
2041  continue;
2042  if(polygon_opt[0]||ruleMap[rule_opt[0]]!=rule::point){
2043  //do not create if no points found within polygon
2044  if(!nPointPolygon)
2045  continue;
2046  //add ring to polygon
2047  if(polygon_opt[0]){
2048  writePolygon.addRing(&writeRing);
2049  writePolygon.closeRings();
2050  //write geometry of writePolygon
2051  if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
2052  cerr << "writing feature failed" << std::endl;
2053  writePolygonFeature->SetGeometry(&writePolygon);
2054  if(verbose_opt[0]>1)
2055  std::cout << "copying new fields write polygon " << std::endl;
2056  if(verbose_opt[0]>1)
2057  std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
2058  //write polygon feature
2059  }
2060  else{//write mean /median value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean, stdev or median )
2061  //create feature
2062  if(verbose_opt[0]>1)
2063  std::cout << "copying fields from polygons " << std::endl;
2064  if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
2065  cerr << "writing feature failed" << std::endl;
2066  writeCentroidFeature->SetGeometry(&writeCentroidPoint);
2067  OGRGeometry *updateGeometry;
2068  updateGeometry = writeCentroidFeature->GetGeometryRef();
2069  assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
2070  if(verbose_opt[0]>1)
2071  std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
2072  }
2073  if(class_opt.empty()){
2074  if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
2075  if(verbose_opt[0])
2076  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2077  for(int index=0;index<polyValues.size();++index){
2078  //test
2079  assert(polyValues[index].size()==1);
2080  double theValue=polyValues[index].back();
2081  if(verbose_opt[0])
2082  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2083  int theBand=(band_opt.size()) ? band_opt[index] : index;
2084 
2085  if(verbose_opt[0]>1)
2086  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
2087  switch( fieldType ){
2088  case OFTInteger:
2089  case OFTReal:
2090  if(polygon_opt[0])
2091  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
2092  else
2093  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
2094  break;
2095  case OFTString:
2096  if(polygon_opt[0])
2097  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
2098  else
2099  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
2100  break;
2101  // case OFTRealList:{
2102  // int fieldIndex;
2103  // int nCount;
2104  // const double *theList;
2105  // vector<double> vectorList;
2106  // if(polygon_opt[0]){
2107  // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
2108  // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
2109  // vectorList.resize(nCount+1);
2110  // for(int index=0;index<nCount;++index)
2111  // vectorList[index]=theList[index];
2112  // vectorList[nCount]=theValue;
2113  // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
2114  // }
2115  // else{
2116  // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
2117  // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
2118  // vectorList.resize(nCount+1);
2119  // for(int index=0;index<nCount;++index)
2120  // vectorList[index]=theList[index];
2121  // vectorList[nCount]=theValue;
2122  // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
2123  // }
2124  // break;
2125  // }
2126  default://not supported
2127  std::cout << "field type not supported yet..." << std::endl;
2128  break;
2129  }
2130  }
2131  }
2132  else{//ruleMap[rule_opt[0]] is not rule::point
2133  double theValue=0;
2134  for(int index=0;index<polyValues.size();++index){
2135  if(ruleMap[rule_opt[0]]==rule::mean)
2136  theValue=stat.mean(polyValues[index]);
2137  else if(ruleMap[rule_opt[0]]==rule::stdev)
2138  theValue=sqrt(stat.var(polyValues[index]));
2139  else if(ruleMap[rule_opt[0]]==rule::median)
2140  theValue=stat.median(polyValues[index]);
2141  else if(ruleMap[rule_opt[0]]==rule::sum)
2142  theValue=stat.sum(polyValues[index]);
2143  else if(ruleMap[rule_opt[0]]==rule::max)
2144  theValue=stat.mymax(polyValues[index]);
2145  else if(ruleMap[rule_opt[0]]==rule::min)
2146  theValue=stat.mymin(polyValues[index]);
2147  else{//rule::centroid
2148  if(verbose_opt[0])
2149  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2150  assert(nPointPolygon<=1);
2151  assert(nPointPolygon==polyValues[index].size());
2152  theValue=polyValues[index].back();
2153  }
2154  if(verbose_opt[0]>1)
2155  std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
2156  switch( fieldType ){
2157  case OFTInteger:
2158  case OFTReal:
2159  if(polygon_opt[0])
2160  writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
2161  else
2162  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
2163  break;
2164  case OFTString:
2165  if(polygon_opt[0])
2166  writePolygonFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
2167  else
2168  writeCentroidFeature->SetField(fieldname_opt[index].c_str(),type2string<double>(theValue).c_str());
2169  break;
2170  // case OFTRealList:{
2171  // int fieldIndex;
2172  // int nCount;
2173  // const double *theList;
2174  // vector<double> vectorList;
2175  // if(polygon_opt[0]){
2176  // fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
2177  // theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
2178  // vectorList.resize(nCount+1);
2179  // for(int index=0;index<nCount;++index)
2180  // vectorList[index]=theList[index];
2181  // vectorList[nCount]=theValue;
2182  // writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
2183  // }
2184  // else{
2185  // fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
2186  // theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
2187  // vectorList.resize(nCount+1);
2188  // for(int index=0;index<nCount;++index)
2189  // vectorList[index]=theList[index];
2190  // vectorList[nCount]=theValue;
2191  // writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
2192  // }
2193  // break;
2194  // }
2195  default://not supported
2196  std::cout << "field type not supported yet..." << std::endl;
2197  break;
2198  }
2199  }
2200  }
2201  }
2202  else{//class_opt is set
2203  if(ruleMap[rule_opt[0]]==rule::proportion){
2204  if(verbose_opt[0])
2205  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2206  stat.normalize_pct(polyClassValues);
2207  for(int index=0;index<polyClassValues.size();++index){
2208  double theValue=polyClassValues[index];
2209  ostringstream fs;
2210  fs << class_opt[index];
2211  if(polygon_opt[0])
2212  writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
2213  else
2214  writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
2215  }
2216  }
2217  else if(ruleMap[rule_opt[0]]==rule::custom){
2218  assert(polygon_opt[0]);//not implemented for points
2219  if(verbose_opt[0])
2220  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2221  stat.normalize_pct(polyClassValues);
2222  assert(polyClassValues.size()==2);//11:broadleaved, 12:coniferous
2223  if(polyClassValues[0]>=75)//broadleaved
2224  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
2225  else if(polyClassValues[1]>=75)//coniferous
2226  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
2227  else if(polyClassValues[0]>25&&polyClassValues[1]>25)//mixed
2228  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
2229  else{
2230  if(verbose_opt[0]){
2231  std::cout << "No valid value in polyClassValues..." << std::endl;
2232  for(int index=0;index<polyClassValues.size();++index){
2233  double theValue=polyClassValues[index];
2234  std::cout << theValue << " ";
2235  }
2236  std::cout << std::endl;
2237  }
2238  writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
2239  }
2240  }
2241  else if(ruleMap[rule_opt[0]]==rule::mode){
2242  //maximum votes in polygon
2243  if(verbose_opt[0])
2244  std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
2245  //search for class with maximum votes
2246  int maxClass=stat.mymin(class_opt);
2247  vector<double>::iterator maxit;
2248  maxit=stat.mymax(polyClassValues,polyClassValues.begin(),polyClassValues.end());
2249  int maxIndex=distance(polyClassValues.begin(),maxit);
2250  maxClass=class_opt[maxIndex];
2251  if(verbose_opt[0]>0)
2252  std::cout << "maxClass: " << maxClass << std::endl;
2253  if(polygon_opt[0])
2254  writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
2255  else
2256  writeCentroidFeature->SetField(label_opt[0].c_str(),maxClass);
2257  }
2258  }
2259 
2260  if(polygon_opt[0]){
2261  if(verbose_opt[0]>1)
2262  std::cout << "creating polygon feature" << std::endl;
2263  if(writeTest){
2264  if(writeTestLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
2265  std::string errorString="Failed to create polygon feature in ogr vector dataset";
2266  throw(errorString);
2267  }
2268  }
2269  else{
2270  if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
2271  std::string errorString="Failed to create polygon feature in ogr vector dataset";
2272  throw(errorString);
2273  }
2274  }
2275  OGRFeature::DestroyFeature( writePolygonFeature );
2276  ++ntotalvalid;
2277  if(verbose_opt[0])
2278  std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
2279  }
2280  else{
2281  if(verbose_opt[0]>1)
2282  std::cout << "creating point feature in centroid" << std::endl;
2283  if(writeTest){
2284  if(writeTestLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
2285  std::string errorString="Failed to create point feature in ogr vector dataset";
2286  throw(errorString);
2287  }
2288  }
2289  else{
2290  //test
2291  assert(validFeature);
2292  if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
2293  std::string errorString="Failed to create point feature in ogr vector dataset";
2294  throw(errorString);
2295  }
2296  }
2297  OGRFeature::DestroyFeature( writeCentroidFeature );
2298  ++ntotalvalid;
2299  if(verbose_opt[0])
2300  std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
2301  }
2302  }
2303  }
2304  else{
2305  std::string test;
2306  test=poGeometry->getGeometryName();
2307  ostringstream oss;
2308  oss << "geometry " << test << " not supported";
2309  throw(oss.str());
2310  }
2311  ++ifeature;
2312  progress=static_cast<float>(ifeature+1)/nfeature;
2313  pfnProgress(progress,pszMessage,pProgressArg);
2314  }
2315  catch(std::string e){
2316  std::cout << e << std::endl;
2317  continue;
2318  }
2319  catch(int npoint){
2320  if(verbose_opt[0])
2321  std::cout << "number of points read in polygon: " << npoint << std::endl;
2322  continue;
2323  }
2324  }//end of getNextFeature
2325  // if(rbox_opt[0]>0||cbox_opt[0]>0)
2326  // boxWriter.close();
2327  progress=1.0;
2328  pfnProgress(progress,pszMessage,pProgressArg);
2329  ++ilayerWrite;
2330  }//for ilayer
2331  ogrWriter.close();
2332  if(sampleIsVirtual)
2333  sampleWriterOgr.close();
2334  if(test_opt.size())
2335  ogrTestWriter.close();
2336  }//else (vector)
2337  progress=1.0;
2338  pfnProgress(progress,pszMessage,pProgressArg);
2339  imgReader.close();
2340 }
2341 
Generated on Wed Nov 19 2014 12:26:59 by   doxygen 1.8.6