[pktools] 183/375: pkcrop corrected bounding box in case of user defined dx dy, support crop for non-georeferenced i mages

Bas Couwenberg sebastic at xs4all.nl
Wed Dec 3 21:54:12 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 9ede8a0c10cc14232720b29cbabbe3d362a3da51
Author: Pieter Kempeneers <kempenep at gmail.com>
Date:   Sun Jan 26 12:52:49 2014 +0100

    pkcrop corrected bounding box in case of user defined dx dy, support crop for non-georeferenced i
    mages
---
 ChangeLog                         |  4 ++++
 src/apps/pkcrop.cc                | 25 +++++++++++++++----------
 src/apps/pkinfo.cc                |  2 +-
 src/imageclasses/ImgReaderGdal.cc | 12 ++++++------
 src/imageclasses/ImgReaderGdal.h  |  4 ++--
 src/imageclasses/ImgWriterGdal.cc | 28 ++++++++++++++--------------
 src/imageclasses/ImgWriterGdal.h  |  4 ++--
 7 files changed, 44 insertions(+), 35 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 72fd89f..9709d29 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -246,8 +246,12 @@ version 2.5
 	new utility (moved from pklas2ig)
 
 version 2.5.1
+ - pkinfo
+	-stats instead of -stat to be conform to gdalinfo
  - pkextract
 	support sum rule
  - pkfilter
 	support of density filter
 	debug down option (forgot to adapt xres and yres)
+ - pkcrop
+	correct bounding box when dx_opt and dy_opt are set
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index 899c415..d27483b 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -18,6 +18,7 @@ You should have received a copy of the GNU General Public License
 along with pktools.  If not, see <http://www.gnu.org/licenses/>.
  ***********************************************************************/
 #include <assert.h>
+#include <cstdlib>
 #include <string>
 #include <list>
 #include <iostream>
@@ -144,8 +145,11 @@ int main(int argc, char *argv[])
   if(dy_opt.size())
     dy=dy_opt[0];
 
+  bool isGeoRef=false;
   for(int iimg=0;iimg<input_opt.size();++iimg){
     imgReader.open(input_opt[iimg]);
+    if(!isGeoRef)
+      isGeoRef=imgReader.isGeoRef();
     if(dx_opt.empty()){
       if(!iimg||imgReader.getDeltaX()<dx)
         dx=imgReader.getDeltaX();
@@ -200,15 +204,15 @@ int main(int argc, char *argv[])
   }
   else if(cx_opt.size()&&cy_opt.size()&&nx_opt.size()&&ny_opt.size()){
     ulx_opt[0]=cx_opt[0]-nx_opt[0]/2.0;
-    uly_opt[0]=cy_opt[0]+ny_opt[0]/2.0;
+    uly_opt[0]=(isGeoRef) ? cy_opt[0]+ny_opt[0]/2.0 : cy_opt[0]-ny_opt[0]/2.0;
     lrx_opt[0]=cx_opt[0]+nx_opt[0]/2.0;
-    lry_opt[0]=cy_opt[0]-ny_opt[0]/2.0;
+    lry_opt[0]=(isGeoRef) ? cy_opt[0]-ny_opt[0]/2.0 : cy_opt[0]+ny_opt[0]/2.0;
   }
   else if(cx_opt.size()&&cy_opt.size()&&ns_opt.size()&&nl_opt.size()){
     ulx_opt[0]=cx_opt[0]-ns_opt[0]*dx/2.0;
-    uly_opt[0]=cy_opt[0]+nl_opt[0]*dy/2.0;
+    uly_opt[0]=(isGeoRef) ? cy_opt[0]+nl_opt[0]*dy/2.0 : cy_opt[0]-nl_opt[0]*dy/2.0;
     lrx_opt[0]=cx_opt[0]+ns_opt[0]*dx/2.0;
-    lry_opt[0]=cy_opt[0]-nl_opt[0]*dy/2.0;
+    lry_opt[0]=(isGeoRef) ? cy_opt[0]-nl_opt[0]*dy/2.0 : cy_opt[0]+nl_opt[0]*dy/2.0;
   }
   if(ulx_opt[0]<cropulx)
     cropulx=ulx_opt[0];
@@ -288,8 +292,9 @@ int main(int argc, char *argv[])
       }
       imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
       imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
-      // ncropcol=ceil((croplrx-cropulx)/dx);
-      // ncroprow=ceil((cropuly-croplry)/dy);
+      //test
+      ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
+      ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
     }
     else{
       double magicX=1,magicY=1;
@@ -309,8 +314,8 @@ int main(int argc, char *argv[])
       imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
       imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
 
-      ncropcol=ceil((croplrx-cropulx)/dx);
-      ncroprow=ceil((cropuly-croplry)/dy);
+      ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
+      ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
       uli=floor(uli);
       ulj=floor(ulj);
       lri=floor(lri);
@@ -373,7 +378,7 @@ int main(int argc, char *argv[])
       gt[2]=0;
       gt[3]=cropuly;
       gt[4]=0;
-      gt[5]=-dy;
+      gt[5]=(imgReader.isGeoRef())? -dy : dy;
       imgWriter.setGeoTransform(gt);
       if(projection_opt.size()){
 	if(verbose_opt[0])
@@ -483,7 +488,7 @@ int main(int argc, char *argv[])
               imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,readBand,theResample);
 	    // for(int ib=0;ib<ncropcol;++ib){
 	    for(int ib=0;ib<imgWriter.nrOfCol();++ib){
-	      assert(imgWriter.image2geo(ib,irow,x,y));
+	      imgWriter.image2geo(ib,irow,x,y);
 	      //lookup corresponding row for irow in this file
 	      imgReader.geo2image(x,y,readCol,readRow);
 	      if(readCol<0||readCol>=imgReader.nrOfCol()){
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index e2d9294..c1dfd41 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -42,7 +42,7 @@ int main(int argc, char *argv[])
   Optionpk<bool>  minmax_opt("mm", "minmax", "Shows min and max value of the image ", false,0);
   Optionpk<bool>  min_opt("min", "min", "Shows min value of the image ", false,0);
   Optionpk<bool>  max_opt("max", "max", "Shows max value of the image ", false,0);
-  Optionpk<bool>  stat_opt("stat", "stat", "Shows statistics (min,max, mean and stdDev of the image)", false,0);
+  Optionpk<bool>  stat_opt("stats", "stats", "Shows statistics (min,max, mean and stdDev of the image)", false,0);
   Optionpk<double>  src_min_opt("src_min", "src_min", "Sets minimum for histogram (does not calculate min value: use -mm instead)");
   Optionpk<double>  src_max_opt("src_max", "src_max", "Sets maximum for histogram (does not calculate min value: use -mm instead)");
   Optionpk<bool>  relative_opt("rel", "rel", "Calculates relative histogram in percentage", false,0);
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index 7a7cc76..6bd8093 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -23,7 +23,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include <iostream>
 
 ImgReaderGdal::ImgReaderGdal(void)
-  : m_gds(NULL), m_isGeoRef(false), m_ncol(0), m_nrow(0), m_nband(0)
+  : m_gds(NULL), m_ncol(0), m_nrow(0), m_nband(0)
 {}
 
 void ImgReaderGdal::open(const std::string& filename)//, double magicX, double magicY)
@@ -56,7 +56,7 @@ void ImgReaderGdal::setCodec()//double magicX, double magicY)
   m_ncol= m_gds->GetRasterXSize();
   m_nrow= m_gds->GetRasterYSize();
   m_nband= m_gds->GetRasterCount();
-  m_isGeoRef=( static_cast<std::string>(m_gds->GetProjectionRef())  != "" );
+  // m_isGeoRef=( static_cast<std::string>(m_gds->GetProjectionRef())  != "" );
   // m_magic_x=magicX;
   // m_magic_y=magicY;
   double adfGeoTransform[6];
@@ -242,7 +242,7 @@ bool ImgReaderGdal::getBoundingBox(double& ulx, double& uly, double& lrx, double
   uly=gt[3];
   lrx=gt[0]+nrOfCol()*gt[1]+nrOfRow()*gt[2];
   lry=gt[3]+nrOfCol()*gt[4]+nrOfRow()*gt[5];
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // ulx=m_ulx;
     // uly=m_uly;
     // lrx=ulx+nrOfCol()*m_delta_x;
@@ -272,7 +272,7 @@ bool ImgReaderGdal::getCentrePos(double& x, double& y) const
   //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
   x=gt[0]+(nrOfCol()/2.0)*gt[1]+(nrOfRow()/2.0)*gt[2];
   y=gt[3]+(nrOfCol()/2.0)*gt[4]+(nrOfRow()/2.0)*gt[5];
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // x=m_ulx+(nrOfCol()/2.0)*m_delta_x;
     // y=m_uly-(nrOfRow()/2.0)*m_delta_y;
     return true;
@@ -304,7 +304,7 @@ bool ImgReaderGdal::geo2image(double x, double y, double& i, double& j) const
     i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
     j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
   }
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // double ulx=m_ulx;
     // double uly=m_uly;
     // i=(x-ulx)/m_delta_x;
@@ -334,7 +334,7 @@ bool ImgReaderGdal::image2geo(double i, double j, double& x, double& y) const
 
   x=gt[0]+(0.5+i)*gt[1]+(0.5+j)*gt[2];
   y=gt[3]+(0.5+i)*gt[4]+(0.5+j)*gt[5];
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // x=m_ulx+(0.5+i)*m_delta_x;
     // y=m_uly-(0.5+j)*m_delta_y;
     return true;
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index e32e02b..adf903c 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -41,7 +41,7 @@ public:
   int nrOfCol(void) const { return m_ncol;};
   int nrOfRow(void) const { return m_nrow;};
   int nrOfBand(void) const { return m_nband;};
-  bool isGeoRef() const {return m_isGeoRef;};
+  bool isGeoRef() const {double gt[6];getGeoTransform(gt);if(gt[5]<0) return true;else return false;};
   std::string getProjection(void) const;
   std::string getProjectionRef(void) const;
   std::string getGeoTransform() const;
@@ -127,7 +127,7 @@ protected:
   /* double m_uly; */
   /* double m_delta_x; */
   /* double m_delta_y; */
-  bool m_isGeoRef;
+  /* bool m_isGeoRef; */
   std::vector<double> m_noDataValues;
   std::vector<double> m_scale;
   std::vector<double> m_offset;
diff --git a/src/imageclasses/ImgWriterGdal.cc b/src/imageclasses/ImgWriterGdal.cc
index 5569bb3..f5f9bc0 100644
--- a/src/imageclasses/ImgWriterGdal.cc
+++ b/src/imageclasses/ImgWriterGdal.cc
@@ -29,7 +29,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 
 //---------------------------------------------------------------------------
 ImgWriterGdal::ImgWriterGdal(void)
-  : m_gds(NULL), m_isGeoRef(false), m_ncol(0), m_nrow(0), m_nband(0)
+  : m_gds(NULL), m_ncol(0), m_nrow(0), m_nband(0)
 {}
 
 // ImgWriterGdal::ImgWriterGdal(void)
@@ -46,7 +46,7 @@ ImgWriterGdal::~ImgWriterGdal(void)
 //---------------------------------------------------------------------------
 void ImgWriterGdal::open(const std::string& filename, const ImgReaderGdal& imgSrc, const std::vector<std::string>& options)
 {
-  m_isGeoRef=imgSrc.isGeoRef();
+  // m_isGeoRef=imgSrc.isGeoRef();
   m_filename=filename;
   m_ncol=imgSrc.nrOfCol();
   m_nrow=imgSrc.nrOfRow();
@@ -76,7 +76,7 @@ void ImgWriterGdal::open(const std::string& filename, const ImgReaderGdal& imgSr
 
 void ImgWriterGdal::open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::vector<std::string>& options)
 {
-  m_isGeoRef=false;
+  // m_isGeoRef=false;
   m_filename = filename;
   m_ncol = ncol;
   m_nrow = nrow;
@@ -124,12 +124,12 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
   // interleaveList << "INTERLEAVE=" << m_interleave;
   // papszOptions = CSLAddString(papszOptions,(interleaveList.str()).c_str());
   m_gds=poDriver->Create(m_filename.c_str(),m_ncol,m_nrow,m_nband,m_type,papszOptions);
-  if(imgSrc.isGeoRef()){
+  // if(imgSrc.isGeoRef()){
     setProjection(imgSrc.getProjection());
     double gt[6];
     imgSrc.getGeoTransform(gt);
     setGeoTransform(gt);
-  }
+  // }
   m_gds->SetMetadata(imgSrc.getMetadata() ); 
   m_gds->SetMetadataItem( "TIFFTAG_DOCUMENTNAME", m_filename.c_str());
   std::string versionString="pktools ";
@@ -263,7 +263,7 @@ std::string ImgWriterGdal::getProjection(void) const
 
 //---------------------------------------------------------------------------
 void ImgWriterGdal::setGeoTransform(double* gt){
-  m_isGeoRef=true;
+  // m_isGeoRef=true;
   m_gt[0]=gt[0];
   m_gt[1]=gt[1];
   m_gt[2]=gt[2];
@@ -304,8 +304,8 @@ void ImgWriterGdal::copyGeoTransform(const ImgReaderGdal& imgSrc)
 
 std::string ImgWriterGdal::setProjectionProj4(const std::string& projection)
 {
-  if(!m_isGeoRef)
-    m_isGeoRef=true;
+  // if(!m_isGeoRef)
+  //   m_isGeoRef=true;
 
   OGRSpatialReference theRef;
   theRef.SetFromUserInput(projection.c_str());
@@ -339,8 +339,8 @@ std::string ImgWriterGdal::setProjectionProj4(const std::string& projection)
 
 void ImgWriterGdal::setProjection(const std::string& projection)
 {
-  if(!m_isGeoRef)
-    m_isGeoRef=true;
+  // if(!m_isGeoRef)
+  //   m_isGeoRef=true;
   OGRSpatialReference oSRS;
   char *pszSRS_WKT = NULL;
   assert(m_gds);
@@ -384,7 +384,7 @@ bool ImgWriterGdal::getBoundingBox(double& ulx, double& uly, double& lrx, double
   uly=gt[3];
   lrx=gt[0]+nrOfCol()*gt[1]+nrOfRow()*gt[2];
   lry=gt[3]+nrOfCol()*gt[4]+nrOfRow()*gt[5];
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // ulx=m_ulx;
     // uly=m_uly;
     // lrx=ulx+nrOfCol()*m_delta_x;
@@ -414,7 +414,7 @@ bool ImgWriterGdal::getCentrePos(double& x, double& y) const
   //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
   x=gt[0]+(nrOfCol()/2.0)*gt[1]+(nrOfRow()/2.0)*gt[2];
   y=gt[3]+(nrOfCol()/2.0)*gt[4]+(nrOfRow()/2.0)*gt[5];
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // x=m_ulx+nrOfCol()/2.0*m_delta_x;
     // y=m_uly-nrOfRow()/2.0*m_delta_y;
     return true;
@@ -442,7 +442,7 @@ bool ImgWriterGdal::geo2image(double x, double y, double& i, double& j) const
     i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
     j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
   }
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // double ulx=m_ulx;
     // double uly=m_uly;
     // i=(x-ulx)/m_delta_x;
@@ -473,7 +473,7 @@ bool ImgWriterGdal::image2geo(double i, double j, double& x, double& y) const
   x=gt[0]+(0.5+i)*gt[1]+(0.5+j)*gt[2];
   y=gt[3]+(0.5+i)*gt[4]+(0.5+j)*gt[5];
 
-  if(m_isGeoRef){
+  if(isGeoRef()){
     // x=m_ulx+(0.5+i)*m_delta_x;
     // y=m_uly-(0.5+j)*m_delta_y;
     return true;
diff --git a/src/imageclasses/ImgWriterGdal.h b/src/imageclasses/ImgWriterGdal.h
index 10fc7a5..98da98b 100644
--- a/src/imageclasses/ImgWriterGdal.h
+++ b/src/imageclasses/ImgWriterGdal.h
@@ -60,7 +60,7 @@ public:
   bool covers(double ulx, double  uly, double lrx, double lry) const;
   bool geo2image(double x, double y, double& i, double& j) const;
   bool image2geo(double i, double j, double& x, double& y) const;
-  bool isGeoRef() const {return m_isGeoRef;};
+  bool isGeoRef() const {double gt[6];getGeoTransform(gt);if(gt[5]<0) return true;else return false;};
   // void getMagicPixel(double& x, double& y) const {x=m_magic_x;y=m_magic_y;};
   double getDeltaX(void) const {double gt[6];getGeoTransform(gt);return gt[1];};
   double getDeltaY(void) const {double gt[6];getGeoTransform(gt);return -gt[5];};
@@ -92,7 +92,7 @@ protected:
   /* double m_uly; */
   /* double m_delta_x; */
   /* double m_delta_y; */
-  bool m_isGeoRef;
+  /* bool m_isGeoRef; */
   // std::string m_interleave;
   // std::string m_compression;
   std::vector<std::string> m_options;

-- 
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