[pktools] 179/375: geotransform should be correct for both projected and not-projected raster images. removed dependence of isGeoRef() in applications
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:11 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit d8e4b7ce156b627c356d10a7efc219795dc3c115
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Fri Jan 10 17:51:40 2014 +0100
geotransform should be correct for both projected and not-projected raster images. removed dependence of isGeoRef() in applications
---
ChangeLog | 1 +
src/algorithms/StatFactory.h | 341 +++++++++++++++++++++++++++++++-------
src/apps/pkcrop.cc | 2 +-
src/apps/pkdiff.cc | 4 +-
src/apps/pkdsm2shadow.cc | 11 +-
src/apps/pkdumpimg.cc | 2 +-
src/apps/pkeditogr.cc | 280 ++++++++++++++++---------------
src/apps/pkenhance.cc | 8 +-
src/apps/pkextract.cc | 11 +-
src/apps/pkfilter.cc | 11 +-
src/apps/pkgetmask.cc | 104 ++++++++----
src/apps/pkinfo.cc | 47 ++----
src/apps/pkndvi.cc | 12 +-
src/apps/pkreclass.cc | 28 +---
src/apps/pksetmask.cc | 34 +---
src/apps/pkstatascii.cc | 43 ++---
src/apps/pkstatogr.cc | 59 ++++---
src/imageclasses/ImgReaderGdal.cc | 29 ++--
src/imageclasses/ImgReaderGdal.h | 4 +-
src/imageclasses/ImgReaderOgr.cc | 146 ++++++++--------
20 files changed, 710 insertions(+), 467 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 6232853..4dd2a19 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -224,6 +224,7 @@ version 2.5
support multiple layers
- pkeditogr
support other ogr file formats than ESRI Shape (e.g, using option -f SQLite)
+ support multiple layers
- Filter2d.h
renamed mask to nodata
- pkinfo
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index 0903e77..c110c66 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -104,7 +104,21 @@ public:
gsl_rng_set(r,theSeed);
return r;
}
-
+ bool isNoData(double value) const{
+ if(m_noDataValues.empty())
+ return false;
+ else
+ return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();
+ };
+ unsigned int pushNodDataValue(double noDataValue){
+ if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
+ m_noDataValues.push_back(noDataValue);
+ return m_noDataValues.size();
+ };
+ unsigned int setNoDataValues(std::vector<double> vnodata){
+ m_noDataValues=vnodata;
+ return m_noDataValues.size();
+ };
static double getRandomValue(const gsl_rng* r, const std::string type, double a=0, double b=0){
std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
initDist(m_distMap);
@@ -120,18 +134,25 @@ public:
}
return randValue;
};
- template<class T> T max(const std::vector<T>& v) const;
template<class T> T min(const std::vector<T>& v) const;
+ template<class T> T max(const std::vector<T>& v) const;
+ template<class T> T min(const std::vector<T>& v, T minConstraint) const;
+ template<class T> T max(const std::vector<T>& v, T maxConstraint) const;
// template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
- template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
- template<class T> typename std::vector<T>::iterator max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const;
- template<class T> typename std::vector<T>::const_iterator absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
template<class T> typename std::vector<T>::const_iterator min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
template<class T> typename std::vector<T>::iterator min(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T minConstraint) const;
+ template<class T> typename std::vector<T>::iterator min(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T minConstraint) const;
+ template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::iterator max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T maxConstraint) const;
+ template<class T> typename std::vector<T>::iterator max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T maxConstraint) const;
template<class T> typename std::vector<T>::const_iterator absmin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
+ template<class T> typename std::vector<T>::const_iterator absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const;
template<class T> void minmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T& theMin, T& theMax) const;
template<class T> T sum(const std::vector<T>& v) const;
template<class T> double mean(const std::vector<T>& v) const;
+ template<class T> T eraseNoData(std::vector<T>& v) const;
template<class T> T median(const std::vector<T>& v) const;
template<class T> double var(const std::vector<T>& v) const;
template<class T> double moment(const std::vector<T>& v, int n) const;
@@ -177,74 +198,195 @@ private:
m_distMap["gaussian"]=gaussian;
m_distMap["uniform"]=uniform;
}
-
+ std::vector<double> m_noDataValues;
};
-template<class T> typename std::vector<T>::iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
+{
+ typename std::vector<T>::const_iterator tmpIt=begin;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(!isNoData(*it))
+ if(*tmpIt>*it)
+ tmpIt=it;
+ }
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
{
typename std::vector<T>::iterator tmpIt=begin;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(!isNoData(*it))
+ if(*tmpIt>*it)
+ tmpIt=it;
+ }
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T minConstraint) const
+{
+ typename std::vector<T>::const_iterator tmpIt=v.end();
+ T minValue=minConstraint;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
+ if((minConstraint<=*it)&&(*it<=minValue)){
+ tmpIt=it;
+ minValue=*it;
+ }
+ }
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T minConstraint) const
+{
+ typename std::vector<T>::iterator tmpIt=v.end();
+ T minValue=minConstraint;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
+ if((minConstraint<=*it)&&(*it<=minValue)){
+ tmpIt=it;
+ minValue=*it;
+ }
+ }
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
+{
+ typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
if(*tmpIt<*it)
tmpIt=it;
}
return tmpIt;
}
-template<class T> T StatFactory::max(const std::vector<T>& v) const
+template<class T> inline typename std::vector<T>::iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) const
{
- T maxValue=*(v.begin());
- for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
- if(maxValue<*it)
+ typename std::vector<T>::iterator tmpIt=begin;
+ for (typename std::vector<T>::iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
+ if(*tmpIt<*it)
+ tmpIt=it;
+ }
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T maxConstraint) const
+{
+ typename std::vector<T>::const_iterator tmpIt=v.end();
+ T maxValue=maxConstraint;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
+ if((maxValue<=*it)&&(*it<=maxConstraint)){
+ tmpIt=it;
maxValue=*it;
+ }
}
- return maxValue;
+ return tmpIt;
+}
+
+template<class T> inline typename std::vector<T>::iterator StatFactory::max(const std::vector<T>& v, typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end, T maxConstraint) const
+{
+ typename std::vector<T>::iterator tmpIt=v.end();
+ T maxValue=maxConstraint;
+ for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
+ if((maxValue<=*it)&&(*it<=maxConstraint)){
+ tmpIt=it;
+ maxValue=*it;
+ }
+ }
+ return tmpIt;
}
-template<class T> T StatFactory::min(const std::vector<T>& v) const
+
+
+
+template<class T> inline T StatFactory::min(const std::vector<T>& v) const
{
T minValue=*(v.begin());
for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ if(isNoData(*it))
+ continue;
if(minValue>*it)
minValue=*it;
}
return minValue;
}
-template<class T> typename std::vector<T>::const_iterator StatFactory::absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
+ template<class T> inline T StatFactory::min(const std::vector<T>& v, T minConstraint) const
{
- typename std::vector<T>::const_iterator tmpIt=begin;
- for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
- if(abs(*tmpIt)<abs(*it))
- tmpIt=it;
+ T minValue=minConstraint;
+ for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ if((minConstraint<=*it)&&(*it<=minValue))
+ minValue=*it;
}
- return tmpIt;
+ return minValue;
}
-template<class T> typename std::vector<T>::const_iterator StatFactory::min(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
+template<class T> inline T StatFactory::max(const std::vector<T>& v) const
+{
+ T maxValue=*(v.begin());
+ for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ if(maxValue<*it)
+ maxValue=*it;
+ }
+ return maxValue;
+}
+
+template<class T> inline T StatFactory::max(const std::vector<T>& v, T maxConstraint) const
+{
+ T maxValue=maxConstraint;
+ for (typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ if((maxValue<=*it)&&(*it<=maxConstraint))
+ maxValue=*it;
+ }
+ return maxValue;
+}
+
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::absmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
- if(*tmpIt>*it)
+ if(isNoData(*it))
+ continue;
+ if(abs(*tmpIt)<abs(*it))
tmpIt=it;
}
return tmpIt;
}
-template<class T> typename std::vector<T>::const_iterator StatFactory::absmin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
+template<class T> inline typename std::vector<T>::const_iterator StatFactory::absmin(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end) const
{
typename std::vector<T>::const_iterator tmpIt=begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
if(abs(*tmpIt)>abs(*it))
tmpIt=it;
}
}
-template<class T> void StatFactory::minmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T& theMin, T& theMax) const
+template<class T> inline void StatFactory::minmax(const std::vector<T>& v, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, T& theMin, T& theMax) const
{
theMin=*begin;
theMax=*begin;
for (typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
+ if(isNoData(*it))
+ continue;
if(theMin>*it)
theMin=*it;
if(theMax<*it)
@@ -252,24 +394,51 @@ template<class T> void StatFactory::minmax(const std::vector<T>& v, typename std
}
}
-template<class T> T StatFactory::sum(const std::vector<T>& v) const
+template<class T> inline T StatFactory::sum(const std::vector<T>& v) const
{
typename std::vector<T>::const_iterator it;
T tmpSum=0;
- for (it = v.begin(); it!= v.end(); ++it)
+ for (it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
tmpSum+=*it;
+ }
return tmpSum;
}
-template<class T> double StatFactory::mean(const std::vector<T>& v) const
+template<class T> inline double StatFactory::mean(const std::vector<T>& v) const
{
assert(v.size());
- return static_cast<double>(sum(v))/v.size();
+ typename std::vector<T>::const_iterator it;
+ T tmpSum=0;
+ unsigned int validSize=0;
+ for (it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ ++validSize;
+ tmpSum+=*it;
+ }
+ if(validSize)
+ return static_cast<double>(tmpSum)/validSize;
+ else
+ return 0;
+}
+
+template<class T> inline T StatFactory::eraseNoData(std::vector<T>& v) const
+{
+ typename std::vector<T>::iterator it;
+ while(it!=v.end()){
+ if(isNoData(*it))
+ v.erase(it);
+ else
+ ++it;
+ }
}
template<class T> T StatFactory::median(const std::vector<T>& v) const
{
std::vector<T> tmpV=v;
+ eraseNoData(tmpV);
sort(tmpV.begin(),tmpV.end());
if(tmpV.size()%2)
return tmpV[tmpV.size()/2];
@@ -280,44 +449,69 @@ template<class T> T StatFactory::median(const std::vector<T>& v) const
template<class T> double StatFactory::var(const std::vector<T>& v) const
{
typename std::vector<T>::const_iterator it;
- double v1=0;
- double m1=mean(v);
- double n=v.size();
- assert(n>1);
- for (it = v.begin(); it!= v.end(); ++it)
- v1+=(*it-m1)*(*it-m1);
- v1/=(n-1);
-// if(v1<0){
-// for (it = v.begin(); it!= v.end(); ++it)
-// std::cout << *it << " ";
-// std::cout << std::endl;
-// }
- assert(v1>=0);
- return v1;
+ unsigned int validSize=0;
+ double m1=0;
+ double m2=0;
+ for (it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ m1+=*it;
+ m2+=(*it)*(*it);
+ ++validSize;
+ }
+ if(validSize){
+ m2/=validSize;
+ m1/=validSize;
+ return m2-m1*m1;
+ }
+ else
+ return 0;
+ /* double v1=0; */
+ /* double m1=mean(v); */
+ /* double n=v.size(); */
+ /* assert(n>1); */
+ /* for (it = v.begin(); it!= v.end(); ++it) */
+ /* v1+=(*it-m1)*(*it-m1); */
+ /* v1/=(n-1); */
+ /* assert(v1>=0); */
+ /* return v1; */
}
template<class T> double StatFactory::moment(const std::vector<T>& v, int n) const
{
assert(v.size());
+ unsigned int validSize=0;
typename std::vector<T>::const_iterator it;
double m=0;
// double m1=mean(v);
for(it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
m+=pow((*it),n);
+ ++validSize;
}
- return m/v.size();
+ if(validSize)
+ return m/validSize;
+ else
+ return 0;
+ /* return m/v.size(); */
}
//central moment
template<class T> double StatFactory::cmoment(const std::vector<T>& v, int n) const
{
assert(v.size());
+ unsigned int validSize=0;
typename std::vector<T>::const_iterator it;
double m=0;
double m1=mean(v);
for(it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
m+=pow((*it-m1),n);
+ ++validSize;
}
+ return m/validSize;
return m/v.size();
}
@@ -336,14 +530,31 @@ template<class T> double StatFactory::kurtosis(const std::vector<T>& v) const
template<class T> void StatFactory::meanVar(const std::vector<T>& v, double& m1, double& v1) const
{
typename std::vector<T>::const_iterator it;
+ unsigned int validSize=0;
+ m1=0;
v1=0;
- m1=mean(v);
- double n=v.size();
- assert(n>1);
- for (it = v.begin(); it!= v.end(); ++it)
- v1+=(*(it)-m1)*(*(it)-m1);
- v1/=(n-1);
- assert(v1>=0);
+ double m2=0;
+ for (it = v.begin(); it!= v.end(); ++it){
+ if(isNoData(*it))
+ continue;
+ m1+=*it;
+ m2+=(*it)*(*it);
+ ++validSize;
+ }
+ if(validSize){
+ m2/=validSize;
+ m1/=validSize;
+ v1=m2-m1*m1;
+ }
+ /* typename std::vector<T>::const_iterator it; */
+ /* v1=0; */
+ /* m1=mean(v); */
+ /* double n=v.size(); */
+ /* assert(n>1); */
+ /* for (it = v.begin(); it!= v.end(); ++it) */
+ /* v1+=(*(it)-m1)*(*(it)-m1); */
+ /* v1/=(n-1); */
+ /* assert(v1>=0); */
}
template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>& input, std::vector<T2>& output, unsigned char lbound, unsigned char ubound) const
@@ -359,18 +570,24 @@ template<class T1, class T2> void StatFactory::scale2byte(const std::vector<T1>&
template<class T> void StatFactory::distribution(const std::vector<T>& input, typename std::vector<T>::const_iterator begin, typename std::vector<T>::const_iterator end, std::vector<double>& output, int nbin, T &minimum, T &maximum, double sigma, const std::string &filename)
{
- if(maximum<=minimum)
- minmax(input,begin,end,minimum,maximum);
- // if(!minimum)
- // minimum=*(min(input,begin,end));
- // if(!maximum)
- // maximum=*(max(input,begin,end));
+ double minValue=0;
+ double maxValue=0;
+ minmax(input,begin,end,minValue,maxValue);
+ if(minimum<maximum&&minimum>minValue)
+ minValue=minimum;
+ if(minimum<maximum&&maximum<maxValue)
+ maxValue=maximum;
+
+ //todo: check...
+ minimum=minValue;
+ maximum=maxValue;
+
if(maximum<=minimum){
std::ostringstream s;
s<<"Error: could not calculate distribution (min>=max)";
throw(s.str());
}
- assert(nbin>1);
+ assert(nbin);
assert(input.size());
if(output.size()!=nbin){
output.resize(nbin);
@@ -378,6 +595,12 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
}
typename std::vector<T>::const_iterator it;
for(it=begin;it!=end;++it){
+ if(*it<minimum)
+ continue;
+ if(*it>maximum)
+ continue;
+ if(isNoData(*it))
+ continue;
if(sigma>0){
// minimum-=2*sigma;
// maximum+=2*sigma;
@@ -394,7 +617,7 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
int theBin=0;
if(*it==maximum)
theBin=nbin-1;
- else if(*it>=minimum && *it<maximum)
+ else if(*it>minimum && *it<maximum)
theBin=static_cast<int>(static_cast<double>((*it)-minimum)/(maximum-minimum)*nbin);
++output[theBin];
// if(*it==maximum)
@@ -411,8 +634,8 @@ template<class T> void StatFactory::distribution(const std::vector<T>& input, t
s<<"Error opening distribution file , " << filename;
throw(s.str());
}
- for(int bin=0;bin<nbin;++bin)
- outputfile << (maximum-minimum)*bin/(nbin-1)+minimum << " " << static_cast<double>(output[bin])/input.size() << std::endl;
+ for(int ibin=0;ibin<nbin;++ibin)
+ outputfile << minimum+static_cast<double>(maximum-minimum)*(ibin+0.5)/nbin << " " << static_cast<double>(output[ibin])/input.size() << std::endl;
outputfile.close();
}
}
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index 45a9ffe..45d2a0e 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -375,7 +375,7 @@ int main(int argc, char *argv[])
cout << "projection: " << projection_opt[0] << endl;
imgWriter.setProjectionProj4(projection_opt[0]);
}
- else if(imgReader.isGeoRef())
+ else
imgWriter.setProjection(imgReader.getProjection());
if(imgWriter.getDataType()==GDT_Byte){
if(colorTable_opt.size()){
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index d47a40d..448cc0d 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -577,9 +577,7 @@ int main(int argc, char *argv[])
gdalWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),1,inputReader.getDataType(),inputReader.getImageType(),option_opt);
if(nodata_opt.size())
gdalWriter.GDALSetNoDataValue(nodata_opt[0]);
- if(inputReader.isGeoRef()){
- gdalWriter.copyGeoTransform(inputReader);
- }
+ gdalWriter.copyGeoTransform(inputReader);
if(colorTable_opt.size())
gdalWriter.setColorTable(colorTable_opt[0]);
else if(inputReader.getColorTable()!=NULL){
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
index 15ede4e..df62063 100644
--- a/src/apps/pkdsm2shadow.cc
+++ b/src/apps/pkdsm2shadow.cc
@@ -117,12 +117,11 @@ int main(int argc,char **argv) {
cout << errorstring << endl;
exit(4);
}
- if(input.isGeoRef()){
- output.setProjection(input.getProjection());
- double gt[6];
- input.getGeoTransform(gt);
- output.setGeoTransform(gt);
- }
+ output.setProjection(input.getProjection());
+ double gt[6];
+ input.getGeoTransform(gt);
+ output.setGeoTransform(gt);
+
if(input.getColorTable()!=NULL)
output.setColorTable(input.getColorTable());
diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc
index 22f16ff..671d4e2 100644
--- a/src/apps/pkdumpimg.cc
+++ b/src/apps/pkdumpimg.cc
@@ -217,7 +217,7 @@ int main(int argc, char *argv[])
//test
// vector<double> readBuffer(readncol);
vector<double> readBuffer(readncol+1);
- assert(imgWriter.isGeoRef());
+ // assert(imgWriter.isGeoRef());
if(band_opt.empty()){
for(int iband=0;iband<imgReader.nrOfBand();++iband)
band_opt.push_back(iband);
diff --git a/src/apps/pkeditogr.cc b/src/apps/pkeditogr.cc
index 8de2f3a..b1331be 100644
--- a/src/apps/pkeditogr.cc
+++ b/src/apps/pkeditogr.cc
@@ -100,9 +100,6 @@ int main(int argc, char *argv[])
}
//start reading features from the layer
- if(verbose_opt[0])
- cout << "reset reading" << endl;
- ogrReader.getLayer()->ResetReading();
// if(field_opt.size())
// assert(field_opt.size()==ogrReader.getFieldCount());
unsigned long int ifeature=0;
@@ -110,105 +107,117 @@ int main(int argc, char *argv[])
cout << "going through features" << endl << flush;
ImgWriterOgr ogrWriter(output_opt[0],ogrformat_opt[0]);
- OGRLayer* writeLayer=ogrWriter.createLayer(output_opt[0],ogrReader.getProjection(),ogrReader.getGeometryType(),NULL);
- std::vector<OGRFieldDefn*> readFields;
- std::vector<OGRFieldDefn*> writeFields;
- ogrReader.getFields(readFields);
- writeFields=readFields;
- try{
- for(int ifield=0;ifield<readFields.size();++ifield){
- // if(field_opt.size()>ifield)
- // writeFields[ifield]->SetName(field_opt[ifield].c_str());
- if(verbose_opt[0])
- std::cout << readFields[ifield]->GetNameRef() << " -> " << writeFields[ifield]->GetNameRef() << std::endl;
- if(writeLayer->CreateField(writeFields[ifield]) != OGRERR_NONE ){
- ostringstream es;
- // if(field_opt.size()>ifield)
- // es << "Creating field " << field_opt[ifield] << " failed";
- // else
- es << "Creating field " << readFields[ifield] << " failed";
- string errorString=es.str();
- throw(errorString);
- }
- }
- }
- catch(string errorString){
- std::cerr << errorString << std::endl;
- exit(1);
- }
+
+ //support multiple layers
+ int nlayer=ogrReader.getLayerCount();
if(verbose_opt[0])
- std::cout << "add " << addname_opt.size() << " fields" << std::endl;
- if(addname_opt.size()){
- assert(addname_opt.size()==addtype_opt.size());
- while(addvalue_opt.size()<addname_opt.size())
- addvalue_opt.push_back(addvalue_opt.back());
- }
- for(int iname=0;iname<addname_opt.size();++iname){
- if(verbose_opt[0])
- std::cout << addname_opt[iname] << " " << std::endl;
- ogrWriter.createField(addname_opt[iname],fieldType[iname]);
- }
- if(verbose_opt[0]){
- std::cout << std::endl;
- std::cout << addname_opt.size() << " fields created" << std::endl;
- }
- const char* pszMessage;
- void* pProgressArg=NULL;
- GDALProgressFunc pfnProgress=GDALTermProgress;
- double progress=0;
- OGRFeature *poFeature;
- unsigned long int nfeature=ogrReader.getFeatureCount();
- while((poFeature = ogrReader.getLayer()->GetNextFeature()) != NULL ){
+ std::cout << "number of layers: " << nlayer << endl;
+
+ for(int ilayer=0;ilayer<nlayer;++ilayer){
+ OGRLayer *readLayer=ogrReader.getLayer(ilayer);
if(verbose_opt[0])
- std::cout << "feature " << ifeature << std::endl;
- ++ifeature;
- bool doSelect;
- if(like_opt.empty())
- doSelect=true;
- else{
- assert(selectField_opt.size());
- int fieldIndex=poFeature->GetFieldIndex(selectField_opt[0].c_str());
- string fieldValue=poFeature->GetFieldAsString(fieldIndex);
- if(stringent_opt[0]){
- if(fieldValue==like_opt[0])
- doSelect=true;
- else
- doSelect=false;
- }
- else{
- doSelect=false;
- for(int ilike=0;ilike<like_opt.size();++ilike){
- if(fieldValue.find(like_opt[ilike])!=std::string::npos){
- if(verbose_opt[0])
- std::cout << "found " << like_opt[ilike] << " in " << fieldValue << std::endl;
- doSelect=true;
- break;
- }
- }
+ cout << "reset reading" << endl;
+ readLayer->ResetReading();
+
+ OGRLayer *writeLayer=ogrWriter.createLayer(output_opt[0],ogrReader.getProjection(),ogrReader.getGeometryType(ilayer),NULL);
+ std::vector<OGRFieldDefn*> readFields;
+ std::vector<OGRFieldDefn*> writeFields;
+ ogrReader.getFields(readFields,ilayer);
+ writeFields=readFields;
+ try{
+ for(int ifield=0;ifield<readFields.size();++ifield){
+ // if(field_opt.size()>ifield)
+ // writeFields[ifield]->SetName(field_opt[ifield].c_str());
+ if(verbose_opt[0])
+ std::cout << readFields[ifield]->GetNameRef() << " -> " << writeFields[ifield]->GetNameRef() << std::endl;
+ if(writeLayer->CreateField(writeFields[ifield]) != OGRERR_NONE ){
+ ostringstream es;
+ // if(field_opt.size()>ifield)
+ // es << "Creating field " << field_opt[ifield] << " failed";
+ // else
+ es << "Creating field " << readFields[ifield] << " failed";
+ string errorString=es.str();
+ throw(errorString);
+ }
}
}
- if(!doSelect){
+ catch(string errorString){
+ std::cerr << errorString << std::endl;
+ exit(1);
+ }
+ if(verbose_opt[0])
+ std::cout << "add " << addname_opt.size() << " fields" << std::endl;
+ if(addname_opt.size()){
+ assert(addname_opt.size()==addtype_opt.size());
+ while(addvalue_opt.size()<addname_opt.size())
+ addvalue_opt.push_back(addvalue_opt.back());
+ }
+ for(int iname=0;iname<addname_opt.size();++iname){
if(verbose_opt[0])
- std::cout << "string not found in feature " << ifeature << std::endl;
- continue;
+ std::cout << addname_opt[iname] << " " << std::endl;
+ ogrWriter.createField(addname_opt[iname],fieldType[iname]);
}
- OGRFeature *poDstFeature = NULL;
- poDstFeature=ogrWriter.createFeature();
- if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
- const char* fmt;
- string errorString="Unable to translate feature %d from layer %s.\n";
- fmt=errorString.c_str();
- CPLError( CE_Failure, CPLE_AppDefined,
- fmt,
- poFeature->GetFID(), ogrWriter.getLayerName().c_str() );
- OGRFeature::DestroyFeature( poFeature );
- OGRFeature::DestroyFeature( poDstFeature );
+ if(verbose_opt[0]){
+ std::cout << std::endl;
+ std::cout << addname_opt.size() << " fields created" << std::endl;
}
- long int fid=poFeature->GetFID();
- poDstFeature->SetFID( poFeature->GetFID() );
- for(int ifeature=0;ifeature<setfeature_opt.size();++ifeature){
- if(fid==setfeature_opt[ifeature]){
- switch(poDstFeature->GetFieldDefnRef(fid)->GetType()){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ OGRFeature *poFeature;
+ unsigned long int nfeature=ogrReader.getFeatureCount(ilayer);
+ while((poFeature = ogrReader.getLayer(ilayer)->GetNextFeature()) != NULL ){
+ if(verbose_opt[0])
+ std::cout << "feature " << ifeature << std::endl;
+ ++ifeature;
+ bool doSelect;
+ if(like_opt.empty())
+ doSelect=true;
+ else{
+ assert(selectField_opt.size());
+ int fieldIndex=poFeature->GetFieldIndex(selectField_opt[0].c_str());
+ string fieldValue=poFeature->GetFieldAsString(fieldIndex);
+ if(stringent_opt[0]){
+ if(fieldValue==like_opt[0])
+ doSelect=true;
+ else
+ doSelect=false;
+ }
+ else{
+ doSelect=false;
+ for(int ilike=0;ilike<like_opt.size();++ilike){
+ if(fieldValue.find(like_opt[ilike])!=std::string::npos){
+ if(verbose_opt[0])
+ std::cout << "found " << like_opt[ilike] << " in " << fieldValue << std::endl;
+ doSelect=true;
+ break;
+ }
+ }
+ }
+ }
+ if(!doSelect){
+ if(verbose_opt[0])
+ std::cout << "string not found in feature " << ifeature << std::endl;
+ continue;
+ }
+ OGRFeature *poDstFeature = NULL;
+ poDstFeature=ogrWriter.createFeature(ilayer);
+ if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
+ const char* fmt;
+ string errorString="Unable to translate feature %d from layer %s.\n";
+ fmt=errorString.c_str();
+ CPLError( CE_Failure, CPLE_AppDefined,
+ fmt,
+ poFeature->GetFID(), ogrWriter.getLayerName().c_str() );
+ OGRFeature::DestroyFeature( poFeature );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ long int fid=poFeature->GetFID();
+ poDstFeature->SetFID( poFeature->GetFID() );
+ for(int ifeature=0;ifeature<setfeature_opt.size();++ifeature){
+ if(fid==setfeature_opt[ifeature]){
+ switch(poDstFeature->GetFieldDefnRef(fid)->GetType()){
case(OFTReal):
poDstFeature->SetField(setname_opt[ifeature].c_str(),string2type<float>(setvalue_opt[ifeature]));
break;
@@ -222,52 +231,53 @@ int main(int argc, char *argv[])
std::cerr << "Error: field type not supported" << std::endl;
exit(1);
break;
+ }
}
}
- }
- //set default values for new fields
- if(verbose_opt[0])
- std::cout << "set default values for new fields in feature " << ifeature << std::endl;
- for(int iname=0;iname<addname_opt.size();++iname){
- switch(fieldType[iname]){
- case(OFTReal):
- if(verbose_opt[0])
- std::cout << "set field " << addname_opt[iname] << " to default " << string2type<float>(addvalue_opt[iname]) << std::endl;
- poDstFeature->SetField(addname_opt[iname].c_str(),string2type<float>(addvalue_opt[iname]));
- break;
- case(OFTInteger):
- if(verbose_opt[0])
- std::cout << "set field " << addname_opt[iname] << " to default " << string2type<int>(addvalue_opt[iname]) << std::endl;
- poDstFeature->SetField(addname_opt[iname].c_str(),string2type<int>(addvalue_opt[iname]));
- break;
- case(OFTString):
- if(verbose_opt[0])
- std::cout << "set field " << addname_opt[iname] << " to default " << addvalue_opt[iname] << std::endl;
- poDstFeature->SetField(addname_opt[iname].c_str(),addvalue_opt[iname].c_str());
- break;
- default:
- std::cerr << "Error: field type not supported" << std::endl;
- exit(1);
- break;
+ //set default values for new fields
+ if(verbose_opt[0])
+ std::cout << "set default values for new fields in feature " << ifeature << std::endl;
+ for(int iname=0;iname<addname_opt.size();++iname){
+ switch(fieldType[iname]){
+ case(OFTReal):
+ if(verbose_opt[0])
+ std::cout << "set field " << addname_opt[iname] << " to default " << string2type<float>(addvalue_opt[iname]) << std::endl;
+ poDstFeature->SetField(addname_opt[iname].c_str(),string2type<float>(addvalue_opt[iname]));
+ break;
+ case(OFTInteger):
+ if(verbose_opt[0])
+ std::cout << "set field " << addname_opt[iname] << " to default " << string2type<int>(addvalue_opt[iname]) << std::endl;
+ poDstFeature->SetField(addname_opt[iname].c_str(),string2type<int>(addvalue_opt[iname]));
+ break;
+ case(OFTString):
+ if(verbose_opt[0])
+ std::cout << "set field " << addname_opt[iname] << " to default " << addvalue_opt[iname] << std::endl;
+ poDstFeature->SetField(addname_opt[iname].c_str(),addvalue_opt[iname].c_str());
+ break;
+ default:
+ std::cerr << "Error: field type not supported" << std::endl;
+ exit(1);
+ break;
+ }
}
- }
- CPLErrorReset();
- if(verbose_opt[0])
- std::cout << "create feature" << std::endl;
- if(ogrWriter.createFeature( poDstFeature ) != OGRERR_NONE){
- const char* fmt;
- string errorString="Unable to translate feature %d from layer %s.\n";
- fmt=errorString.c_str();
- CPLError( CE_Failure, CPLE_AppDefined,
- fmt,
- poFeature->GetFID(), ogrWriter.getLayerName().c_str() );
+ CPLErrorReset();
+ if(verbose_opt[0])
+ std::cout << "create feature" << std::endl;
+ if(ogrWriter.createFeature( poDstFeature,ilayer ) != OGRERR_NONE){
+ const char* fmt;
+ string errorString="Unable to translate feature %d from layer %s.\n";
+ fmt=errorString.c_str();
+ CPLError( CE_Failure, CPLE_AppDefined,
+ fmt,
+ poFeature->GetFID(), ogrWriter.getLayerName().c_str() );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ OGRFeature::DestroyFeature( poFeature );
OGRFeature::DestroyFeature( poDstFeature );
+ progress=static_cast<float>(ifeature+1)/nfeature;
+ pfnProgress(progress,pszMessage,pProgressArg);
}
- OGRFeature::DestroyFeature( poFeature );
- OGRFeature::DestroyFeature( poDstFeature );
- progress=static_cast<float>(ifeature+1)/nfeature;
- pfnProgress(progress,pszMessage,pProgressArg);
}
if(verbose_opt[0])
std::cout << "replaced " << ifeature << " features" << std::endl;
diff --git a/src/apps/pkenhance.cc b/src/apps/pkenhance.cc
index 1fb9539..a518f79 100644
--- a/src/apps/pkenhance.cc
+++ b/src/apps/pkenhance.cc
@@ -136,8 +136,8 @@ int main(int argc,char **argv) {
pfnProgress(progress,pszMessage,pProgressArg);
for(int iband=0;iband<nband;++iband){
//calculate histograms
- int nbinRef=nbin_opt[0];
- int nbinInput=nbin_opt[0];
+ unsigned int nbinRef=nbin_opt[0];
+ unsigned int nbinInput=nbin_opt[0];
std::vector<unsigned long int> histRef(nbinRef);
std::vector<unsigned long int> histInput(nbinInput);
double minValueRef=0;
@@ -155,10 +155,10 @@ int main(int argc,char **argv) {
unsigned long int nsampleRef=referenceImg.getHistogram(histRef,minValueRef,maxValueRef,nbinRef,iband);
unsigned long int nsampleInput=inputImg.getHistogram(histInput,minValueInput,maxValueInput,nbinInput,iband);
//create cumulative historgrams
- for(int bin=0;bin<nbinRef;++bin){
+ for(unsigned int bin=0;bin<nbinRef;++bin){
histRef[bin]+=100.0*static_cast<double>(histRef[bin])/static_cast<double>(nsampleRef);
}
- for(int bin=0;bin<nbinInput;++bin)
+ for(unsigned int bin=0;bin<nbinInput;++bin)
histInput[bin]+=100.0*static_cast<double>(histInput[bin])/static_cast<double>(nsampleInput);
//match histograms
vector<double> lineBuffer(inputImg.nrOfCol());
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 5d6820a..ed911e7 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -42,14 +42,14 @@ using namespace std;
int main(int argc, char *argv[])
{
Optionpk<string> image_opt("i", "image", "Input image file");
- Optionpk<string> sample_opt("s", "sample", "Input sample file (shape) or class file (e.g. Corine CLC) if class option is set");
+ Optionpk<string> sample_opt("s", "sample", "Input sample vector file or class file (e.g. Corine CLC) if class option is set");
Optionpk<string> mask_opt("m", "mask", "Mask image file");
Optionpk<int> msknodata_opt("msknodata", "msknodata", "Mask value where image is invalid. If a single mask is used, more nodata values can be set. If more masks are used, use one value for each mask.", 1);
Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Use -1 to process every class in sample image, or leave empty to extract all valid data pixels from sample file");
Optionpk<string> output_opt("o", "output", "Output sample file (image file)");
Optionpk<string> ogrformat_opt("f", "f", "Output sample file format","ESRI Shapefile");
Optionpk<string> test_opt("test", "test", "Test sample file (use this option in combination with threshold<100 to create a training (output) and test set");
- Optionpk<string> bufferOutput_opt("bu", "bu", "Buffer output shape file");
+ // Optionpk<string> bufferOutput_opt("bu", "bu", "Buffer output shape file");
Optionpk<short> geo_opt("g", "geo", "geo coordinates", 1);
Optionpk<short> down_opt("down", "down", "down sampling factor. Can be used to create grid points", 1);
Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold for vector sample files. If using raster land cover maps as a sample file, 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);
@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
Optionpk<short> disc_opt("circ", "circular", "circular disc kernel boundary", 0);
Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or Integer)", "Real");
Optionpk<string> ltype_opt("lt", "ltype", "Label type: In16 or String", "Integer");
- Optionpk<string> fieldname_opt("bn", "bname", "Field name of output shape file", "B");
+ Optionpk<string> fieldname_opt("bn", "bname", "Attribute field name of extracted raster band", "B");
Optionpk<string> label_opt("cn", "cname", "name of the class label in the output vector file", "label");
Optionpk<bool> polygon_opt("polygon", "polygon", "create OGRPolygon as geometry instead of OGRPoint. Only if sample option is also of polygon type.", false);
Optionpk<int> band_opt("b", "band", "band index to crop. Use -1 to use all bands)", -1);
@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
output_opt.retrieveOption(argc,argv);
ogrformat_opt.retrieveOption(argc,argv);
test_opt.retrieveOption(argc,argv);
- bufferOutput_opt.retrieveOption(argc,argv);
+ // bufferOutput_opt.retrieveOption(argc,argv);
geo_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
@@ -773,8 +773,7 @@ int main(int argc, char *argv[])
std::cout << "number of layers: " << nlayer << endl;
for(int ilayer=0;ilayer<nlayer;++ilayer){
- OGRLayer *readLayer;
- readLayer = sampleReaderOgr.getDataSource()->GetLayer(ilayer);
+ OGRLayer *readLayer=sampleReaderOgr.getLayer(ilayer);
cout << "processing layer " << readLayer->GetName() << endl;
readLayer->ResetReading();
OGRLayer *writeLayer;
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index e24c5f6..0a0ce6c 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -195,12 +195,11 @@ int main(int argc,char **argv) {
cout << errorstring << endl;
exit(4);
}
- if(input.isGeoRef()){
- output.setProjection(input.getProjection());
- double gt[6];
- input.getGeoTransform(gt);
- output.setGeoTransform(gt);
- }
+ output.setProjection(input.getProjection());
+ double gt[6];
+ input.getGeoTransform(gt);
+ output.setGeoTransform(gt);
+
if(colorTable_opt.size()){
if(colorTable_opt[0]!="none"){
if(verbose_opt[0])
diff --git a/src/apps/pkgetmask.cc b/src/apps/pkgetmask.cc
index d7ecda2..8c0a861 100644
--- a/src/apps/pkgetmask.cc
+++ b/src/apps/pkgetmask.cc
@@ -27,8 +27,8 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<string> input_opt("i", "input", "Input image file");
Optionpk<short> band_opt("b", "band", "band(s) used for mask", 0);
- Optionpk<double> min_opt("min", "min", "Values smaller than min threshold(s) are masked as invalid. Use one threshold for each band", 0);
- Optionpk<double> max_opt("max", "max", "Values greater than max threshold(s) are masked as invalid. Use one threshold for each band", 0);
+ Optionpk<double> min_opt("min", "min", "Values smaller than min threshold(s) are masked as invalid. Use one threshold for each band");
+ Optionpk<double> max_opt("max", "max", "Values greater than max threshold(s) are masked as invalid. Use one threshold for each band");
Optionpk<string> operator_opt("p", "operator", "Operator: [AND,OR].", "OR");
Optionpk<unsigned short> data_opt("data", "data", "value(s) for valid pixels: between min and max", 1);
Optionpk<unsigned short> nodata_opt("nodata", "nodata", "value(s) for invalid pixels: not between min and max", 0);
@@ -93,23 +93,34 @@ int main(int argc,char **argv) {
assert(band_opt.size()>=0);
assert(band_opt.size()<=imgReader.nrOfBand());
- while(band_opt.size()>min_opt.size())
- min_opt.push_back(min_opt[0]);
- while(band_opt.size()>max_opt.size())
- max_opt.push_back(max_opt[0]);
- while(min_opt.size()>data_opt.size())
- data_opt.push_back(data_opt[0]);
- assert(min_opt.size()==max_opt.size());
- if(verbose_opt[0]){
- cout << "min,max values: ";
- for(int imin=0;imin<min_opt.size();++imin){
- cout << min_opt[imin] << "," << max_opt[imin];
- if(imin<min_opt.size()-1)
- cout << " ";
- else
- cout << endl;
- }
+ if(min_opt.size()&&max_opt.size()){
+ if(min_opt.size()!=max_opt.size())
+ cerr << "Error: number of min and max options must correspond if both min and max options are provide" << endl;
+ assert(min_opt.size()==max_opt.size());
+ }
+ if(min_opt.size()){
+ while(band_opt.size()>min_opt.size())
+ min_opt.push_back(min_opt[0]);
+ while(min_opt.size()>data_opt.size())
+ data_opt.push_back(data_opt[0]);
+ }
+ if(max_opt.size()){
+ while(band_opt.size()>max_opt.size())
+ max_opt.push_back(max_opt[0]);
+ while(max_opt.size()>data_opt.size())
+ data_opt.push_back(data_opt[0]);
}
+ // assert(min_opt.size()==max_opt.size());
+ // if(verbose_opt[0]){
+ // cout << "min,max values: ";
+ // for(int imin=0;imin<min_opt.size();++imin){
+ // cout << min_opt[imin] << "," << max_opt[imin];
+ // if(imin<min_opt.size()-1)
+ // cout << " ";
+ // else
+ // cout << endl;
+ // }
+ // }
vector< vector<float> > lineBuffer(band_opt.size());
for(int iband=0;iband<band_opt.size();++iband)
@@ -137,12 +148,12 @@ int main(int argc,char **argv) {
}
else if (imgReader.getColorTable()!=NULL)//copy colorTable from input image
imgWriter.setColorTable(imgReader.getColorTable());
- if(imgReader.isGeoRef()){
- imgWriter.setProjection(imgReader.getProjection());
- double gt[6];
- imgReader.getGeoTransform(gt);
- imgWriter.setGeoTransform(gt);//ulx,uly,imgReader.getDeltaX(),imgReader.getDeltaY(),0,0);
- }
+
+ imgWriter.setProjection(imgReader.getProjection());
+ double gt[6];
+ imgReader.getGeoTransform(gt);
+ imgWriter.setGeoTransform(gt);//ulx,uly,imgReader.getDeltaX(),imgReader.getDeltaY(),0,0);
+
if(nodata_opt.size())
imgWriter.GDALSetNoDataValue(nodata_opt[0]);
@@ -153,15 +164,42 @@ int main(int argc,char **argv) {
for(int icol=0;icol<imgReader.nrOfCol();++icol){
bool valid=(operator_opt[0]=="OR")?false:true;
unsigned short validValue=data_opt[0];
- for(int ivalid=0;ivalid<min_opt.size();++ivalid){
- bool validBand=false;
- // for(int iband=0;iband<band_opt.size();++iband){
- unsigned short theBand=(band_opt.size()==min_opt.size())? ivalid:0;
- if(lineBuffer[theBand][icol]>=min_opt[ivalid]&&lineBuffer[theBand][icol]<=max_opt[ivalid]){
- validValue=data_opt[ivalid];
- validBand=true;
- }
- valid=(operator_opt[0]=="OR")?valid||validBand : valid&&validBand;
+ if(min_opt.size()&&max_opt.size()){
+ assert(max_opt.size()==min_opt.size());
+ for(int ivalid=0;ivalid<min_opt.size();++ivalid){
+ bool validBand=false;
+ // for(int iband=0;iband<band_opt.size();++iband){
+ unsigned short theBand=(band_opt.size()==min_opt.size())? ivalid:0;
+ if(lineBuffer[theBand][icol]>=min_opt[ivalid]&&lineBuffer[theBand][icol]<=max_opt[ivalid]){
+ validValue=data_opt[ivalid];
+ validBand=true;
+ }
+ valid=(operator_opt[0]=="OR")?valid||validBand : valid&&validBand;
+ }
+ }
+ else if(min_opt.size()){
+ for(int ivalid=0;ivalid<min_opt.size();++ivalid){
+ bool validBand=false;
+ // for(int iband=0;iband<band_opt.size();++iband){
+ unsigned short theBand=(band_opt.size()==min_opt.size())? ivalid:0;
+ if(lineBuffer[theBand][icol]>=min_opt[ivalid]){
+ validValue=data_opt[ivalid];
+ validBand=true;
+ }
+ valid=(operator_opt[0]=="OR")?valid||validBand : valid&&validBand;
+ }
+ }
+ else if(max_opt.size()){
+ for(int ivalid=0;ivalid<max_opt.size();++ivalid){
+ bool validBand=false;
+ // for(int iband=0;iband<band_opt.size();++iband){
+ unsigned short theBand=(band_opt.size()==max_opt.size())? ivalid:0;
+ if(lineBuffer[theBand][icol]<=max_opt[ivalid]){
+ validValue=data_opt[ivalid];
+ validBand=true;
+ }
+ valid=(operator_opt[0]=="OR")?valid||validBand : valid&&validBand;
+ }
}
if(valid)
writeBuffer[icol]=validValue;
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 60f7d90..a0ffa5a 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box");
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box");
Optionpk<bool> hist_opt("hist", "hist", "Calculates histogram. Use --rel for a relative histogram output. ", false,0);
- Optionpk<short> nbin_opt("nbin", "nbin", "Number of bins used in histogram. Use 0 for all input values as integers", 0,0);
+ Optionpk<unsigned int> nbin_opt("nbin", "nbin", "Number of bins used in histogram. Use 0 for all input values as integers");
Optionpk<bool> type_opt("ot", "otype", "Returns data type", false,0);
Optionpk<bool> description_opt("d", "description", "Returns image description", false,0);
Optionpk<bool> metadata_opt("meta", "meta", "Shows meta data ", false,0);
@@ -289,7 +289,7 @@ int main(int argc, char *argv[])
hist_opt[0]=true;
if(hist_opt[0]){
assert(band_opt[0]<imgReader.nrOfBand());
- int nbin=nbin_opt[0];
+ unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
std::vector<unsigned long int> output;
minValue=0;
maxValue=0;
@@ -304,35 +304,24 @@ int main(int argc, char *argv[])
std::cout.precision(10);
for(int bin=0;bin<nbin;++bin){
// nsample+=output[bin];
- if(nbin>1){
- std::cout << (maxValue-minValue)*bin/(nbin-1)+minValue << " ";
- if(relative_opt[0])
- std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
- else
- std::cout << static_cast<double>(output[bin]) << std::endl;
- }
- else{
- std::cout << (maxValue-minValue)*bin/(2-1)+minValue << " ";
- if(relative_opt[0])
- std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
- else
- std::cout << static_cast<double>(output[bin]) << std::endl;
- }
-
- }
- }
- else{
- int minCol,minRow;
- if(src_min_opt.size()){
- assert(band_opt[0]<imgReader.nrOfBand());
- std::cout << "--min " << imgReader.getMin(minCol, minRow,band_opt[0]);
- }
- if(src_max_opt.size()){
- assert(band_opt[0]<imgReader.nrOfBand());
- assert(band_opt[0]<imgReader.nrOfBand());
- std::cout << "--max " << imgReader.getMax(minCol, minRow,band_opt[0]);
+ std::cout << minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin << " ";
+ if(relative_opt[0])
+ std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
+ else
+ std::cout << static_cast<double>(output[bin]) << std::endl;
}
}
+ // else{
+ // int minCol,minRow;
+ // if(src_min_opt.size()){
+ // assert(band_opt[0]<imgReader.nrOfBand());
+ // std::cout << "--min " << imgReader.getMin(minCol, minRow,band_opt[0]);
+ // }
+ // if(src_max_opt.size()){
+ // assert(band_opt[0]<imgReader.nrOfBand());
+ // std::cout << "--max " << imgReader.getMax(minCol, minRow,band_opt[0]);
+ // }
+ // }
if(projection_opt[0]){
if(imgReader.isGeoRef())
std::cout << " -a_srs " << imgReader.getProjection() << " ";
diff --git a/src/apps/pkndvi.cc b/src/apps/pkndvi.cc
index 3cf338f..4ac2b1e 100644
--- a/src/apps/pkndvi.cc
+++ b/src/apps/pkndvi.cc
@@ -150,12 +150,12 @@ int main(int argc, char *argv[])
if(description_opt.size())
outputWriter.setImageDescription(description_opt[0]);
//if input image is georeferenced, copy projection info to output image
- if(inputReader[0].isGeoRef()){
- outputWriter.setProjection(inputReader[0].getProjection());
- double ulx,uly,lrx,lry;
- inputReader[0].getBoundingBox(ulx,uly,lrx,lry);
- outputWriter.copyGeoTransform(inputReader[0]);
- }
+
+ outputWriter.setProjection(inputReader[0].getProjection());
+ double ulx,uly,lrx,lry;
+ inputReader[0].getBoundingBox(ulx,uly,lrx,lry);
+ outputWriter.copyGeoTransform(inputReader[0]);
+
if(colorTable_opt.size()){
if(colorTable_opt[0]!="none")
outputWriter.setColorTable(colorTable_opt[0]);
diff --git a/src/apps/pkreclass.cc b/src/apps/pkreclass.cc
index 21fefd8..e800bfc 100644
--- a/src/apps/pkreclass.cc
+++ b/src/apps/pkreclass.cc
@@ -227,14 +227,8 @@ int main(int argc, char *argv[])
if(inputReader.isGeoRef()){
for(int imask=0;imask<mask_opt.size();++imask)
assert(maskReader[imask].isGeoRef());
- outputWriter.setProjection(inputReader.getProjection());
- }
- else{
- for(int imask=0;imask<mask_opt.size();++imask){
- assert(maskReader[imask].nrOfCol()==inputReader.nrOfCol());
- assert(maskReader[imask].nrOfRow()==inputReader.nrOfRow());
- }
}
+ outputWriter.setProjection(inputReader.getProjection());
double ulx,uly,lrx,lry;
inputReader.getBoundingBox(ulx,uly,lrx,lry);
outputWriter.copyGeoTransform(inputReader);
@@ -277,14 +271,8 @@ int main(int argc, char *argv[])
bool masked=false;
if(mask_opt.size()>1){//multiple masks
for(int imask=0;imask<mask_opt.size();++imask){
- if(maskReader[imask].isGeoRef()){
- inputReader.image2geo(icol,irow,x,y);
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- }
- else{
- colMask=icol;
- rowMask=irow;
- }
+ inputReader.image2geo(icol,irow,x,y);
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
assert(rowMask>=0&&rowMask<maskReader[imask].nrOfRow());
try{
@@ -310,14 +298,8 @@ int main(int argc, char *argv[])
}
}
else if(mask_opt.size()){//potentially more invalid values for single mask
- if(maskReader[0].isGeoRef()){
- inputReader.image2geo(icol,irow,x,y);
- maskReader[0].geo2image(x,y,colMask,rowMask);
- }
- else{
- colMask=icol;
- rowMask=irow;
- }
+ inputReader.image2geo(icol,irow,x,y);
+ maskReader[0].geo2image(x,y,colMask,rowMask);
if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
assert(rowMask>=0&&rowMask<maskReader[0].nrOfRow());
try{
diff --git a/src/apps/pksetmask.cc b/src/apps/pksetmask.cc
index 841d7c7..6d21e60 100644
--- a/src/apps/pksetmask.cc
+++ b/src/apps/pksetmask.cc
@@ -130,12 +130,6 @@ int main(int argc, char *argv[])
for(int imask=0;imask<mask_opt.size();++imask)
assert(maskReader[imask].isGeoRef());
}
- else{
- for(int imask=0;imask<mask_opt.size();++imask){
- assert(maskReader[imask].nrOfCol()==inputReader.nrOfCol());
- assert(maskReader[imask].nrOfRow()==inputReader.nrOfRow());
- }
- }
assert(nodata_opt.size()==msknodata_opt.size());
assert(operator_opt.size()==msknodata_opt.size()||operator_opt.size()==1);
if(verbose_opt[0]){
@@ -185,16 +179,10 @@ int main(int argc, char *argv[])
for(icol=0;icol<inputReader.nrOfCol();++icol){
if(mask_opt.size()>1){//multiple masks
for(int imask=0;imask<mask_opt.size();++imask){
- if(maskReader[imask].isGeoRef()){
- inputReader.image2geo(icol,irow,x,y);
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- colMask=static_cast<int>(colMask);
- rowMask=static_cast<int>(rowMask);
- }
- else{
- colMask=icol;
- rowMask=irow;
- }
+ inputReader.image2geo(icol,irow,x,y);
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
+ colMask=static_cast<int>(colMask);
+ rowMask=static_cast<int>(rowMask);
bool masked=false;
if(rowMask>=0&&rowMask<maskReader[imask].nrOfRow()&&colMask>=0&&colMask<maskReader[imask].nrOfCol()){
if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask[imask])){
@@ -252,16 +240,10 @@ int main(int argc, char *argv[])
}
}
else{//potentially more invalid values for single mask
- if(maskReader[0].isGeoRef()){
- inputReader.image2geo(icol,irow,x,y);
- maskReader[0].geo2image(x,y,colMask,rowMask);
- colMask=static_cast<int>(colMask);
- rowMask=static_cast<int>(rowMask);
- }
- else{
- colMask=icol;
- rowMask=irow;
- }
+ inputReader.image2geo(icol,irow,x,y);
+ maskReader[0].geo2image(x,y,colMask,rowMask);
+ colMask=static_cast<int>(colMask);
+ rowMask=static_cast<int>(rowMask);
bool masked=false;
if(rowMask>=0&&rowMask<maskReader[0].nrOfRow()&&colMask>=0&&colMask<maskReader[0].nrOfCol()){
if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask[0])){
diff --git a/src/apps/pkstatascii.cc b/src/apps/pkstatascii.cc
index 7e9a233..71ba3c2 100644
--- a/src/apps/pkstatascii.cc
+++ b/src/apps/pkstatascii.cc
@@ -50,13 +50,13 @@ int main(int argc, char *argv[])
Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false);
Optionpk<bool> min_opt("min","min","calculate minimum value",false);
Optionpk<bool> max_opt("max","max","calculate maximum value",false);
- Optionpk<double> src_min_opt("src_min","src_min","start reading source from this minimum value",0);
- Optionpk<double> src_max_opt("src_max","src_max","stop reading source from this maximum value",0);
+ Optionpk<double> src_min_opt("src_min","src_min","start reading source from this minimum value");
+ Optionpk<double> src_max_opt("src_max","src_max","stop reading source from this maximum value");
Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false);
Optionpk<bool> histogram2d_opt("hist2d","hist2d","calculate 2-dimensional histogram based on two columns",false);
- Optionpk<short> nbin_opt("nbin","nbin","number of bins to calculate histogram",0);
+ Optionpk<short> nbin_opt("nbin","nbin","number of bins to calculate histogram");
Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false);
- Optionpk<double> kde_opt("kde","kde","bandwith of kernel density when producing histogram, use 0 for practical estimation based on Silverman's rule of thumb");
+ Optionpk<double> kde_opt("kde","kde","bandwith of kernel density when producing histogram, use 0 for practical estimation based on Silverman's rule of thumb. Leave empty if no kernel density is required");
Optionpk<bool> correlation_opt("cor","correlation","calculate Pearson produc-moment correlation coefficient between two columns (defined by -c <col1> -c <col2>",false);
Optionpk<bool> rmse_opt("rmse","rmse","calculate root mean square error between two columns (defined by -c <col1> -c <col2>",false);
Optionpk<bool> reg_opt("reg","regression","calculate linear regression error between two columns (defined by -c <col1> -c <col2>",false);
@@ -129,14 +129,16 @@ int main(int argc, char *argv[])
asciiReader.setMaxRow(range_opt[1]);
asciiReader.readData(dataVector,col_opt);
assert(dataVector.size());
- double minValue=src_min_opt[0];
- double maxValue=src_max_opt[0];
+ double minValue=0;
+ double maxValue=0;
+
+ unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
if(histogram_opt[0]){
- if(nbin_opt[0]<1){
+ if(nbin<1){
std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
if(maxValue<=minValue)
stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minValue,maxValue);
- nbin_opt[0]=maxValue-minValue+1;
+ nbin=maxValue-minValue+1;
}
}
double minX=src_min_opt[0];
@@ -145,7 +147,7 @@ int main(int argc, char *argv[])
double maxY=(src_max_opt.size()==2)? src_max_opt[1] : src_max_opt[0];
if(histogram2d_opt[0]){
assert(col_opt.size()==2);
- if(nbin_opt[0]<1){
+ if(nbin<1){
std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
if(maxValue<=minValue){
stat.minmax(dataVector[0],dataVector[0].begin(),dataVector[0].end(),minX,maxX);
@@ -155,7 +157,7 @@ int main(int argc, char *argv[])
maxValue=(maxX>maxY)? maxX:maxY;
if(verbose_opt[0])
std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
- nbin_opt[0]=maxValue-minValue+1;
+ nbin=maxValue-minValue+1;
}
}
for(int icol=0;icol<col_opt.size();++icol){
@@ -198,14 +200,14 @@ int main(int argc, char *argv[])
else
sigma=1.06*sqrt(stat.var(dataVector[icol]))*pow(dataVector[icol].size(),-0.2);
}
- assert(nbin_opt[0]);
+ assert(nbin);
if(verbose_opt[0]){
if(sigma>0)
std::cout << "calculating kernel density estimate with sigma " << sigma << " for col " << icol << std::endl;
else
std::cout << "calculating histogram for col " << icol << std::endl;
}
- stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin_opt[0],minValue,maxValue,sigma);
+ stat.distribution(dataVector[icol],dataVector[icol].begin(),dataVector[icol].end(),statVector[icol],nbin,minValue,maxValue,sigma);
if(verbose_opt[0])
std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
}
@@ -227,7 +229,8 @@ int main(int argc, char *argv[])
}
if(histogram_opt[0]){
for(int irow=0;irow<statVector.begin()->size();++irow){
- std::cout << (maxValue-minValue)*irow/(nbin_opt[0]-1)+minValue << " ";
+
+ std::cout << minValue+static_cast<double>(maxValue-minValue)*(irow+0.5)/nbin << " ";
for(int icol=0;icol<col_opt.size();++icol){
if(relative_opt[0])
std::cout << 100.0*static_cast<double>(statVector[icol][irow])/static_cast<double>(dataVector[icol].size());
@@ -240,7 +243,7 @@ int main(int argc, char *argv[])
}
}
if(histogram2d_opt[0]){
- assert(nbin_opt[0]);
+ assert(nbin);
assert(dataVector.size()==2);
assert(dataVector[0].size()==dataVector[1].size());
double sigma=0;
@@ -251,22 +254,22 @@ int main(int argc, char *argv[])
else
sigma=1.06*sqrt(sqrt(stat.var(dataVector[0]))*sqrt(stat.var(dataVector[0])))*pow(dataVector[0].size(),-0.2);
}
- assert(nbin_opt[0]);
+ assert(nbin);
if(verbose_opt[0]){
if(sigma>0)
std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for cols " << col_opt[0] << " and " << col_opt[1] << std::endl;
else
std::cout << "calculating 2d histogram for cols " << col_opt[0] << " and " << col_opt[1] << std::endl;
- std::cout << "nbin: " << nbin_opt[0] << std::endl;
+ std::cout << "nbin: " << nbin << std::endl;
}
std::vector< std::vector<double> > histVector;
- stat.distribution2d(dataVector[0],dataVector[1],histVector,nbin_opt[0],minX,maxX,minY,maxY,sigma);
- for(int binX=0;binX<nbin_opt[0];++binX){
+ stat.distribution2d(dataVector[0],dataVector[1],histVector,nbin,minX,maxX,minY,maxY,sigma);
+ for(int binX=0;binX<nbin;++binX){
std::cout << std::endl;
- for(int binY=0;binY<nbin_opt[0];++binY){
+ for(int binY=0;binY<nbin;++binY){
double value=0;
value=static_cast<double>(histVector[binX][binY])/dataVector[0].size();
- std::cout << (maxX-minX)*binX/(nbin_opt[0]-1)+minX << " " << (maxY-minY)*binY/(nbin_opt[0]-1)+minY << " " << value << std::endl;
+ std::cout << minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin << " " << minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin << " " << value << std::endl;
}
}
}
diff --git a/src/apps/pkstatogr.cc b/src/apps/pkstatogr.cc
index a088185..139e721 100644
--- a/src/apps/pkstatogr.cc
+++ b/src/apps/pkstatogr.cc
@@ -32,11 +32,13 @@ int main(int argc, char *argv[])
Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false);
Optionpk<bool> min_opt("min","min","calculate minimum value",0);
Optionpk<bool> max_opt("max","max","calculate maximum value",0);
- Optionpk<double> src_min_opt("min","min","set minimum value",0);
- Optionpk<double> src_max_opt("max","max","set maximum value",0);
+ Optionpk<double> src_min_opt("src_min","src_min","set minimum value for histogram");
+ Optionpk<double> src_max_opt("src_max","src_max","set maximum value for histogram");
+ Optionpk<double> nodata_opt("nodata","nodata","set nodata value(s)");
Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false);
- Optionpk<short> nbin_opt("nbin", "nbin", "number of bins", 0);
+ Optionpk<unsigned int> nbin_opt("nbin", "nbin", "number of bins");
Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false);
+ Optionpk<double> kde_opt("kde","kde","bandwith of kernel density when producing histogram, use 0 for practical estimation based on Silverman's rule of thumb. Leave empty if no kernel density is required");
Optionpk<bool> mean_opt("mean","mean","calculate mean value",false);
Optionpk<bool> median_opt("median","median","calculate median value",false);
Optionpk<bool> stdev_opt("stdev","stdev","calculate standard deviation",false);
@@ -52,9 +54,11 @@ int main(int argc, char *argv[])
max_opt.retrieveOption(argc,argv);
src_min_opt.retrieveOption(argc,argv);
src_max_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
histogram_opt.retrieveOption(argc,argv);
nbin_opt.retrieveOption(argc,argv);
relative_opt.retrieveOption(argc,argv);
+ kde_opt.retrieveOption(argc,argv);
mean_opt.retrieveOption(argc,argv);
median_opt.retrieveOption(argc,argv);
stdev_opt.retrieveOption(argc,argv);
@@ -83,23 +87,34 @@ int main(int argc, char *argv[])
statfactory::StatFactory stat;
//todo: implement ALL
+ stat.setNoDataValues(nodata_opt);
for(int ifield=0;ifield<fieldname_opt.size();++ifield){
if(verbose_opt[0])
std::cout << "field: " << ifield << std::endl;
theData.clear();
inputReader.readData(theData,OFTReal,fieldname_opt[ifield],0,verbose_opt[0]);
std::vector<double> binData;
- double minValue=src_min_opt[0];
- double maxValue=src_max_opt[0];
+ double minValue=0;
+ double maxValue=0;
+ stat.minmax(theData,theData.begin(),theData.end(),minValue,maxValue);
+ if(src_min_opt.size())
+ minValue=src_min_opt[0];
+ if(src_max_opt.size())
+ maxValue=src_max_opt[0];
+ unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
+
if(histogram_opt[0]){
- if(nbin_opt[0]<1){
- if(maxValue<=minValue)
- stat.minmax(theData,theData.begin(),theData.end(),minValue,maxValue);
- nbin_opt[0]=maxValue-minValue+1;
+ double sigma=0;
+ if(kde_opt.size()){
+ if(kde_opt[0]>0)
+ sigma=kde_opt[0];
+ else
+ sigma=1.06*sqrt(stat.var(theData))*pow(theData.size(),-0.2);
}
- assert(nbin_opt[0]);
+ if(nbin<1)
+ nbin=(maxValue-minValue+1);
try{
- stat.distribution(theData,theData.begin(),theData.end(),binData,nbin_opt[0],minValue,maxValue);
+ stat.distribution(theData,theData.begin(),theData.end(),binData,nbin,minValue,maxValue,sigma);
}
catch(std::string theError){
std::cerr << "Warning: all identical values in data" << std::endl;
@@ -116,25 +131,27 @@ int main(int argc, char *argv[])
std::cout << " --mean " << theMean;
if(stdev_opt[0])
std::cout << " --stdev " << sqrt(theVar);
- if(minmax_opt[0]){
- std::cout << " -min " << stat.min(theData);
- std::cout << " -max " << stat.max(theData);
+ if(minmax_opt[0]||min_opt[0]||max_opt[0]){
+ if(minmax_opt[0])
+ std::cout << " --min " << minValue << " --max " << maxValue << " ";
+ else{
+ if(min_opt[0])
+ std::cout << " --min " << minValue << " ";
+ if(max_opt[0])
+ std::cout << " --max " << maxValue << " ";
+ }
}
- if(min_opt[0])
- std::cout << " -min " << stat.min(theData);
- if(max_opt[0])
- std::cout << " -max " << stat.max(theData);
if(median_opt[0])
std::cout << " -median " << stat.median(theData);
if(size_opt[0])
std::cout << " -size " << theData.size();
std::cout << std::endl;
if(histogram_opt[0]){
- for(int bin=0;bin<nbin_opt[0];++bin){
+ for(int ibin=0;ibin<nbin;++ibin){
if(relative_opt[0])
- std::cout << (maxValue-minValue)*bin/(nbin_opt[0]-1)+minValue << " " << 100.0*static_cast<double>(binData[bin])/theData.size() << std::endl;
+ std::cout << minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin << " " << 100.0*static_cast<double>(binData[ibin])/theData.size() << std::endl;
else
- std::cout << (maxValue-minValue)*bin/(nbin_opt[0]-1)+minValue << " " << binData[bin] << std::endl;
+ std::cout << minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin << " " << binData[ibin] << std::endl;
}
}
}
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index 9d0a9a9..0c8acc1 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -370,8 +370,8 @@ double ImgReaderGdal::getMin(int& x, int& y, int band) const{
for(int irow=0;irow<nrOfRow();++irow){
readData(lineBuffer,GDT_Float64,irow,band);
for(int icol=0;icol<nrOfCol();++icol){
- bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
- if(valid){
+ // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
+ if(!isNoData(lineBuffer[icol])){
if(!init){
y=irow;
x=icol;
@@ -399,8 +399,9 @@ double ImgReaderGdal::getMax(int& x, int& y, int band) const{
for(int irow=0;irow<nrOfRow();++irow){
readData(lineBuffer,GDT_Float64,irow,band);
for(int icol=0;icol<nrOfCol();++icol){
- bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
- if(valid){
+ // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
+ // if(valid){
+ if(!isNoData(lineBuffer[icol])){
if(!init){
y=irow;
x=icol;
@@ -442,8 +443,9 @@ void ImgReaderGdal::getMinMax(int startCol, int endCol, int startRow, int endRow
for(int irow=startCol;irow<endRow+1;++irow){
readData(lineBuffer,GDT_Float64,startCol,endCol,irow,band);
for(int icol=0;icol<lineBuffer.size();++icol){
- bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
- if(valid){
+ // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
+ // if(valid){
+ if(!isNoData(lineBuffer[icol])){
if(!init){
minValue=lineBuffer[icol];
maxValue=lineBuffer[icol];
@@ -482,8 +484,9 @@ void ImgReaderGdal::getMinMax(double& minValue, double& maxValue, int band, bool
for(int irow=0;irow<nrOfRow();++irow){
readData(lineBuffer,GDT_Float64,irow,band);
for(int icol=0;icol<nrOfCol();++icol){
- bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
- if(valid){
+ // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
+ // if(valid){
+ if(!isNoData(lineBuffer[icol])){
if(!init){
minValue=lineBuffer[icol];
maxValue=lineBuffer[icol];
@@ -503,7 +506,7 @@ void ImgReaderGdal::getMinMax(double& minValue, double& maxValue, int band, bool
}
}
-unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max, int& nbin, int theBand) const{
+unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max, unsigned int& nbin, int theBand) const{
double minValue=0;
double maxValue=0;
getMinMax(minValue,maxValue,theBand);
@@ -517,7 +520,7 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi
if(maxValue>minValue){
if(nbin==0)
nbin=maxValue-minValue+1;
- scale=static_cast<double>(nbin-1)/(maxValue-minValue);
+ scale=static_cast<double>(nbin)/(maxValue-minValue);
}
else
nbin=1;
@@ -541,7 +544,6 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi
else{//scale to [0:nbin]
int theBin=static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue));
assert(theBin>=0);
- assert(theBin!=nbin);
assert(theBin<nbin);
++histvector[theBin];
// else if(lineBuffer[icol]==maxValue)
@@ -606,8 +608,9 @@ void ImgReaderGdal::getRefPix(double& refX, double &refY, int band) const
for(int irow=0;irow<nrOfRow();++irow){
readData(lineBuffer,GDT_Float64,irow,band);
for(int icol=0;icol<nrOfCol();++icol){
- bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
- if(valid){
+ // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
+ // if(valid){
+ if(!isNoData(lineBuffer[icol])){
validCol+=icol+1;
++nvalidCol;
validRow+=irow+1;
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index 21c1da0..e32e02b 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -78,7 +78,7 @@ public:
/* }; */
}
int getNoDataValues(std::vector<double>& noDataValues) const;
- bool isNoData(double value) const{return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
+ bool isNoData(double value) const{if(m_noDataValues.empty()) return false;else return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
int pushNoDataValue(double noDataValue);
CPLErr GDALSetNoDataValue(double noDataValue, int band=0) {getRasterBand(band)->SetNoDataValue(noDataValue);};
bool covers(double x, double y) const;
@@ -97,7 +97,7 @@ public:
void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double& minValue, double& maxValue) const;
void getMinMax(double& minValue, double& maxValue, int band=0, bool exhaustiveSearch=false) const;
double getMin(int& col, int& row, int band=0) const;
- unsigned long int getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max,int& nbin, int theBand=0) const;
+ unsigned long int getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max,unsigned int& nbin, int theBand=0) const;
double getMax(int& col, int& row, int band=0) const;
void getRefPix(double& refX, double &refY, int band=0) const;
void getRange(std::vector<short>& range, int Band=0) const;
diff --git a/src/imageclasses/ImgReaderOgr.cc b/src/imageclasses/ImgReaderOgr.cc
index 7adf657..500dac9 100644
--- a/src/imageclasses/ImgReaderOgr.cc
+++ b/src/imageclasses/ImgReaderOgr.cc
@@ -229,48 +229,48 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
int nband=0;
if(verbose)
std::cout << "reading shape file " << m_filename << std::endl;
- try{
- //only retain bands in fields
- getFields(fields);
- std::vector<std::string>::iterator fit=fields.begin();
- if(verbose>1)
- std::cout << "reading fields: ";
- while(fit!=fields.end()){
- if(verbose)
- std::cout << *fit << " ";
+ for(int ilayer=0;ilayer<getLayerCount();++ilayer){
+ try{
+ //only retain bands in fields
+ getFields(fields,ilayer);
+ std::vector<std::string>::iterator fit=fields.begin();
+ if(verbose>1)
+ std::cout << "reading fields: ";
+ while(fit!=fields.end()){
+ if(verbose)
+ std::cout << *fit << " ";
// size_t pos=(*fit).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ");
- if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
- if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
- int theBand=atoi((*fit).substr(1).c_str());
- if(bands.size()){
- bool validBand=false;
- for(int iband=0;iband<bands.size();++iband){
- if(theBand==bands[iband])
- validBand=true;
+ if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
+ if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
+ int theBand=atoi((*fit).substr(1).c_str());
+ if(bands.size()){
+ bool validBand=false;
+ for(int iband=0;iband<bands.size();++iband){
+ if(theBand==bands[iband])
+ validBand=true;
+ }
+ if(validBand)
+ ++fit;
+ else
+ fields.erase(fit);
}
- if(validBand)
- ++fit;
else
- fields.erase(fit);
+ ++fit;
}
- else
+ else if((*fit)=="B" || (*fit)=="b" || (*fit)=="Band")//B is only band
++fit;
}
- else if((*fit)=="B" || (*fit)=="b" || (*fit)=="Band")//B is only band
- ++fit;
+ else
+ fields.erase(fit);
}
- else
- fields.erase(fit);
- }
- if(verbose)
- std::cout << std::endl;
- if(verbose){
- std::cout << "fields:";
+ if(verbose)
+ std::cout << std::endl;
+ if(verbose){
+ std::cout << "fields:";
for(std::vector<std::string>::iterator fit=fields.begin();fit!=fields.end();++fit)
std::cout << " " << *fit;
std::cout << std::endl;
- }
- for(int ilayer=0;ilayer<m_datasource->GetLayerCount();++ilayer){
+ }
if(!nband){
if(verbose)
std::cout << "reading data" << std::endl;
@@ -284,11 +284,11 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
if(verbose)
std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
}
- }
- catch(std::string e){
- std::ostringstream estr;
- estr << e << " " << m_filename;
- throw(estr.str());
+ catch(std::string e){
+ std::ostringstream estr;
+ estr << e << " " << m_filename;
+ throw(estr.str());
+ }
}
if(verbose)
std::cout << "total number of samples read " << totalSamples << std::endl;
@@ -308,39 +308,39 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
int nband=0;
if(verbose)
std::cout << "reading shape file " << m_filename << std::endl;
- try{
- //only retain bands in fields
- getFields(fields);
- std::vector<std::string>::iterator fit=fields.begin();
- if(verbose)
- std::cout << "reading fields: ";
- while(fit!=fields.end()){
+ for(int ilayer=0;ilayer<getLayerCount();++ilayer){
+ try{
+ //only retain bands in fields
+ getFields(fields,ilayer);
+ std::vector<std::string>::iterator fit=fields.begin();
if(verbose)
- std::cout << *fit << " ";
- // size_t pos=(*fit).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ");
- if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
- if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
- int iband=atoi((*fit).substr(1).c_str());
- if((start||end)&&(iband<start||iband>end))
- fields.erase(fit);
- else
+ std::cout << "reading fields: ";
+ while(fit!=fields.end()){
+ if(verbose)
+ std::cout << *fit << " ";
+ // size_t pos=(*fit).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ");
+ if((*fit).substr(0,1)=="B"||(*fit).substr(0,1)=="b"){
+ if((*fit).substr(1).find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_ ")!=std::string::npos){
+ int iband=atoi((*fit).substr(1).c_str());
+ if((start||end)&&(iband<start||iband>end))
+ fields.erase(fit);
+ else
+ ++fit;
+ }
+ else if(*fit=="B" || *fit=="b"|| *fit=="Band")
++fit;
}
- else if(*fit=="B" || *fit=="b"|| *fit=="Band")
- ++fit;
+ else
+ fields.erase(fit);
+ }
+ if(verbose)
+ std::cout << std::endl;
+ if(verbose){
+ std::cout << "fields:";
+ for(std::vector<std::string>::iterator fit=fields.begin();fit!=fields.end();++fit)
+ std::cout << " " << *fit;
+ std::cout << std::endl;
}
- else
- fields.erase(fit);
- }
- if(verbose)
- std::cout << std::endl;
- if(verbose){
- std::cout << "fields:";
- for(std::vector<std::string>::iterator fit=fields.begin();fit!=fields.end();++fit)
- std::cout << " " << *fit;
- std::cout << std::endl;
- }
- for(int ilayer=0;ilayer<m_datasource->GetLayerCount();++ilayer){
if(!nband){
if(verbose)
std::cout << "reading data" << std::endl;
@@ -354,14 +354,14 @@ unsigned int ImgReaderOgr::readDataImageShape(std::map<std::string,Vector2d<floa
if(verbose)
std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
}
+ catch(std::string e){
+ std::ostringstream estr;
+ estr << e << " " << m_filename;
+ throw(estr.str());
+ }
+ if(verbose)
+ std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
}
- catch(std::string e){
- std::ostringstream estr;
- estr << e << " " << m_filename;
- throw(estr.str());
- }
- if(verbose)
- std::cout << ": " << nsample << " samples read with " << nband << " bands" << std::endl;
if(verbose)
std::cout << "total number of samples read " << totalSamples << std::endl;
return totalSamples;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list