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