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