Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

conversion.cc

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /****************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/conversion.cc,v $
00032  * $Revision: 3.14 $
00033  * $Date: 2005/04/06 15:29:03 $
00034  * $Author: kampes $
00035  *
00036  * Some utilities like deg2rad
00037  * we need also UTM2WGS etc here to include other DEMs (class?)
00038  ****************************************************************/
00039 
00040 
00041 #include "matrixbk.hh"
00042 
00043 #include "constants.hh"                 // structs, SOL
00044 #include "readinput.hh"                 // input structs
00045 #include "orbitbk.hh"                    // my orbit class
00046 #include "slcimage.hh"                   // my slc image class
00047 #include "exceptions.hh"                 // my exceptions class
00048 
00049 #include "conversion.hh"                // header file, inlines
00050 #include "ioroutines.hh"                // error etc.
00051 #include "utilities.hh"                 // splint etc.
00052 #include "refsystems.hh"                // RADIUS
00053 
00054 #include <cmath>                        // sin, sqrt, rint, etc.
00055 
00056 
00057 /****************************************************************
00058  *    pol2xyz                                                   *
00059  *                                                              *
00060  * Converts     lat,lon,hei(above surface) (sphere)             * 
00061  *         to   x,y,z (cartesian)                               *
00062  * uses constants in constants.h (pi, radius)                   *
00063  *                                                              *
00064  * lat is positive for North, lon is negative for West.         *
00065  *                                                              *
00066  *    Bert Kampes, 21-Dec-1998                                  *
00067  ****************************************************************/
00068 void pol2xyz(
00069         cn &xyz,
00070         real8 phi,
00071         real8 lon,
00072         real8 height) //default =0)
00073   {
00074   TRACE_FUNCTION("pol2xyz (BK 21-Dec-1998)");
00075   real8 Rph = RADIUS + height;
00076   xyz.x = Rph *cos(phi)*cos(lon);
00077   xyz.y = Rph *cos(phi)*sin(lon);
00078   xyz.z = Rph *sin(phi);
00079   }  // END pol2xyz
00080 
00081 
00082 /****************************************************************
00083  *    xyz2pol                                                   *
00084  *                                                              *
00085  * Converts x,y,z (cartesian)                                   *
00086  *           to lat,lon,hei(above surface) (sphere)             * 
00087  * uses constants in constants.h (pi, radius)                   *
00088  *                                                              *
00089  * lat is positive for North, lon is negative for West.         *
00090  *                                                              *
00091  *    Bert Kampes, 21-Dec-1998                                  *
00092  ****************************************************************/
00093 void xyz2pol(
00094         const cn &xyz,
00095         real8    &phi,
00096         real8    &lambda,
00097         real8    &height)
00098   {
00099   TRACE_FUNCTION("xyz2pol (BK 21-Dec-1998)");
00100   real8 r = sqrt(sqr(xyz.x)+sqr(xyz.y));
00101   phi     = atan2(xyz.z,r);
00102   lambda  = atan2(xyz.y,xyz.x);
00103   height  = (r/cos(phi)) - RADIUS;
00104   }  // END xyz2pol
00105 
00106 
00107 
00108 /****************************************************************
00109  *    xyz2ell                                                   *
00110  *                                                              *
00111  * Converts geocentric cartesian coord. in the XXXX             *
00112  *  reference frame to geodetic ones, using                     *
00113  *  method of bowring see globale en locale geodetische systemen*
00114  * input:                                                       *
00115  *  - ellipsinfo, xyz, (phi,lam,hei)                            *
00116  * output:                                                      *
00117  *  - void (updated lam<-pi,pi> ?, phi<-pi,pi>,hei)             *
00118  *                                                              *
00119  *                                                              *
00120  *    Bert Kampes, 05-Jan-1999                                  *
00121  ****************************************************************/
00122 void xyz2ell(
00123         const input_ell &ell,
00124         const cn        &xyz,
00125         real8           &phi,
00126         real8           &lambda,
00127         real8           &height)
00128   {
00129   TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999)");
00130   real8 r    = sqrt(sqr(xyz.x)+sqr(xyz.y));
00131   real8 nu   = atan2((xyz.z*ell.a),(r*ell.b));
00132   real8 sin3 = pow(sin(nu),3);
00133   real8 cos3 = pow(cos(nu),3);
00134  
00135   phi        = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
00136   lambda     = atan2(xyz.y,xyz.x);
00137   real8 N    = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
00138   height     = (r/cos(phi)) - N;
00139   } // END xyz2ell
00140 
00141 
00142 
00143 /****************************************************************
00144  *    xyz2ell                                                   *
00145  *                                                              *
00146  * Converts geocentric cartesian coord. in the XXXX             *
00147  *  reference frame to geodetic ones, using                     *
00148  *  method of bowring see globale en locale geodetische systemen*
00149  * input:                                                       *
00150  *  - ellipsinfo, xyz, (phi,lam)                                *
00151  * output:                                                      *
00152  *  - void (updated lam<-pi,pi> ?, phi<-pi,pi>)                 *
00153  *                                                              *
00154  *                                                              *
00155  *    Bert Kampes, 05-Jan-1999                                  *
00156  ****************************************************************/
00157 void xyz2ell(
00158         const input_ell &ell,
00159         const cn        &xyz,
00160         real8           &phi,
00161         real8           &lambda)
00162   {
00163   TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999");
00164   real8 r    = sqrt(sqr(xyz.x)+sqr(xyz.y));
00165   real8 nu   = atan2((xyz.z*ell.a),(r*ell.b));
00166   real8 sin3 = pow(sin(nu),3);
00167   real8 cos3 = pow(cos(nu),3);
00168  
00169   phi        = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
00170   lambda     = atan2(xyz.y,xyz.x);
00171   // real8 N    = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
00172   // height     = (r/cos(phi)) - N;
00173   } // END xyz2ell
00174 
00175 
00176 
00177 /****************************************************************
00178  *    ell2xyz                                                   *
00179  *                                                              *
00180  * Converts wgs84 ellipsoid cn to geocentric cartesian coord.   *
00181  * input:                                                       *
00182  *  - ellipsinfo, phi,lam,hei (x,y,z)                           *
00183  * output:                                                      *
00184  *  - void (updated xyz)                                        *
00185  *                                                              *
00186  *                                                              *
00187  *    Bert Kampes, 05-Jan-1999                                  *
00188  ****************************************************************/
00189 void ell2xyz(
00190         const input_ell &ell,
00191         cn    &xyz,
00192         real8  phi,
00193         real8  lambda,
00194         real8  height)
00195   {
00196 /*
00197   // --- remove some slow io, since called very often for step ---
00198   // --- COMPREFDEM #%// Bert Kampes, 24-Sep-2004 ---
00199   #ifdef __DEBUG
00200   TRACE_FUNCTION("ell2xyz (BK 05-Jan-1999)");
00201   #endif
00202   real8 N   = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
00203   real8 Nph = N + height;
00204   #ifdef __DEBUG
00205   DEBUG << "N = " << N;
00206   DEBUG.print();
00207   #endif
00208 
00209   // ______Return these______
00210   xyz.x =  Nph * cos(phi)  * cos(lambda); 
00211   xyz.y =  Nph * cos(phi)  * sin(lambda); 
00212   xyz.z = (Nph - ell.e2*N) * sin(phi);
00213 */
00214 
00215   const real8 sinphi = sin(phi);// tmp variable
00216   const real8 N      = ell.a / sqrt(1.0-ell.e2*sqr(sinphi));
00217   const real8 Nph    = N + height;
00218   const real8 Nph_cosphi = Nph*cos(phi);// tmp variable
00219   // ______Return updated xyz______
00220   xyz.x =  Nph_cosphi      * cos(lambda);
00221   xyz.y =  Nph_cosphi      * sin(lambda);
00222   xyz.z = (Nph - ell.e2*N) * sinphi;
00223   } // END ell2xyz
00224 
00225 
00226 
00227 
00228 /****************************************************************
00229  *    tiepoint                                                  *
00230  *                                                              *
00231  * For a given tiepoint, compute the pixel coordinates and      *
00232  * timing etc.  This can then be used in matlab to correct the  *
00233  * unwrapped phase by hand.                                     *
00234  * This routine is called only if tiepoint is given.            *
00235  *                                                              *
00236  * Input:                                                       *
00237  *  lat/lon/hei of a tiepoint                                   *
00238  * Output:                                                      *
00239  *  INFO the details                                            *
00240  *                                                              *
00241  #%// BK 14-May-2004                                            *
00242  ****************************************************************/
00243 void tiepoint(
00244         const input_gen     &generalinput,
00245         const slcimage      &master,
00246         const slcimage      &slave,
00247               orbit         &masterorbit,
00248               orbit         &slaveorbit,
00249         const input_ell     &ellips)
00250   {
00251   TRACE_FUNCTION("tiepoint (Bert Kampes 14-May-2004)")
00252   if (abs(generalinput.tiepoint.x) < 0.001 &&
00253       abs(generalinput.tiepoint.y) < 0.001 &&
00254       abs(generalinput.tiepoint.z) < 0.001)
00255     {
00256     INFO.print("No tiepoint given.");
00257     return;
00258     }
00259   if (masterorbit.npoints()==0)
00260     {
00261     INFO.print("Exiting tiepoint, no orbit data master available.");
00262     return;
00263     }
00264   DEBUG.precision(30); // permanent; INFO << setp(13)<<q;// once
00265   DEBUG.width(14);
00266   INFO.precision(13); // permanent; INFO << setp(13)<<q;// once
00267   INFO.width(14);
00268   const int32 MAXITER   = 10;
00269   const real8 CRITERPOS = 1e-6;
00270   const real8 CRITERTIM = 1e-10;
00271   const real8 m_minpi4cdivlam = (-4.0*PI*SOL)/master.wavelength;
00272   const real8 s_minpi4cdivlam = (-4.0*PI*SOL)/slave.wavelength;
00273   DEBUG << "Master wavelength = " << master.wavelength;
00274   DEBUG.print();
00275   DEBUG << "Slave wavelength  = " << slave.wavelength;
00276   DEBUG.print();
00277 
00278 
00279   ;// ______ Convert to X,Y,Z ______
00280   INFO.print("Computing radar coordinates and unwrapped phase of tie point");
00281   real8 tiepointphi    = deg2rad(generalinput.tiepoint.x);// lat in [dec.degrees]
00282   real8 tiepointlambda = deg2rad(generalinput.tiepoint.y);// lon in [dec.degrees]
00283   real8 tiepointheight = generalinput.tiepoint.z;//          height in [m]
00284   DEBUG << "TIEPOINT: lat/lon/hei [rad/rad/m]: " 
00285         << tiepointphi << " " << tiepointlambda << " " << tiepointheight;
00286   DEBUG.print();
00287   cn tiepointpos;// compute phi/lambda/h-->x,y,z
00288   ell2xyz(ellips,tiepointpos,tiepointphi,tiepointlambda,tiepointheight);
00289   INFO  << "TIEPOINT lat/lon/hei --> x,y,z (WGS84): "
00290         << tiepointphi << ", " << tiepointlambda << ", " << tiepointheight
00291         << " --> "
00292         << tiepointpos.x << " " << tiepointpos.y << " " << tiepointpos.z;
00293   INFO.print();
00294 
00295   // ______ Convert to line/pix in master ______
00296   real8 tiepointline;// returned
00297   real8 tiepointpix;// returned
00298   int32 n_iter = xyz2lp(tiepointline, tiepointpix,
00299                         master, masterorbit, tiepointpos,
00300                         MAXITER, CRITERTIM);
00301   INFO << "TIEPOINT line/pix in master: " << tiepointline << " " << tiepointpix;
00302   INFO.print();
00303 
00304   // ______ Compute master azimuth/range time of this pixel______
00305   const real8 m_aztime = master.line2ta(tiepointline);
00306   const real8 m_trange = master.pix2tr(tiepointpix);
00307   DEBUG << "TIEPOINT master azimuth time: " << m_aztime;
00308   DEBUG.print();
00309   DEBUG << "TIEPOINT master range time:   " << m_trange;
00310   DEBUG.print();
00311 
00312   // ______ Ellipsoid reference phase for this point master(l,p) ______
00313   // ______ Compute xyz for slave satellite from P ______
00314   cn P_ref;// point at H=0 (ellips)
00315   lp2xyz(tiepointline,tiepointpix,
00316          ellips,master,masterorbit,P_ref,
00317          MAXITER,CRITERPOS);
00318   INFO  << "Pixel of TIEPOINT on ellipsoid: x,y,z (WGS84): "
00319         << P_ref.x << " " << P_ref.y << " " << P_ref.z;
00320   INFO.print();
00321 
00322   // ______ Check if slave is there ______
00323   if (slaveorbit.npoints()==0)
00324     {
00325     INFO.print("tiepoint: no orbit data slave available cannot compute more.");
00326     DEBUG.reset();
00327     INFO.reset();
00328     }
00329   else
00330     {
00331     // ______ Compute xyz for slave satellite from P ______
00332     real8 s_aztime;                             // returned
00333     real8 s_trange;                             // returned
00334     xyz2t(s_aztime,s_trange,
00335           slave, slaveorbit, tiepointpos,MAXITER,CRITERTIM);
00336     DEBUG << "TIEPOINT slave azimuth time: " << s_aztime;
00337     DEBUG.print();
00338     DEBUG << "TIEPOINT slave range time:   " << s_trange;
00339     DEBUG.print();
00340     INFO << "TIEPOINT line/pix in slave: " 
00341          << slave.ta2line(s_aztime) << " "
00342          << slave.tr2pix (s_trange);
00343     INFO.print();
00344 
00345     // ___ Phase of this point ___
00346     real8 phase = m_minpi4cdivlam*m_trange - s_minpi4cdivlam*s_trange;// [rad]
00347     INFO << "TIEPOINT unwrapped phase (incl. reference phase): "
00348          << phase;
00349     INFO.print();
00350   
00351     // _____ Times for slave image ______
00352     real8 s_tazi_ref;                               // returned
00353     real8 s_trange_ref;                             // returned
00354     xyz2t(s_tazi_ref,s_trange_ref,slave,
00355           slaveorbit, P_ref,
00356           MAXITER,CRITERTIM);
00357 
00358     // ___ report reference phase corrected unwrapped phase of tiepoint ___
00359     // ___ real8 m_trange_ref = m_trange;// by definition
00360     // ___ t_corrected = t_m - t_m_ref - (t_s - t_s_ref) = t_s_ref-t_s;
00361     real8 phase_ref_corrected = s_minpi4cdivlam*(s_trange_ref-s_trange);// [rad]
00362     INFO << "TIEPOINT unwrapped phase (ellipsoid ref.phase corrected) [rad]: "
00363          << phase_ref_corrected;
00364     INFO.print();
00365     INFO.print("This phase can be used to correct the unwrapped interferogram");
00366     INFO.print("make sure that you compute the correct pixel position");
00367     INFO.print("the line/pix here are in original master coordinates.");
00368 
00369 
00370     // ______ Check if zero-Doppler iteration is correct by going back ______
00371     // ------ Go from P_ref to M and to S, then back from M to P_refcheckM
00372     // ------ and from S to P_refcheckS; Bert Kampes, 06-Apr-2005
00373     real8 t_az, t_rg;// tmp variables
00374     // --- from P_ref to orbits ---
00375     // ______ master: from P_ref to M ____
00376     n_iter     = xyz2t(t_az,t_rg,master,masterorbit,P_ref,MAXITER,CRITERTIM);// return t
00377     const cn M = masterorbit.getxyz(t_az);
00378     // ______ master: from M to P_ref_checkM ____
00379     cn P_ref_checkM;
00380     n_iter     = lp2xyz(master.ta2line(t_az), master.tr2pix(t_rg),
00381                     ellips, master, masterorbit, 
00382                     P_ref_checkM, MAXITER, CRITERPOS);
00383     // ______ slave: from P_ref to S ____
00384     n_iter     = xyz2t(t_az,t_rg,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);// return t
00385     const cn S = slaveorbit.getxyz(t_az);
00386     // ______ slave: from S to P_ref_checkS ____
00387     cn P_ref_checkS;
00388     n_iter     = lp2xyz(slave.ta2line(t_az), slave.tr2pix(t_rg),
00389                     ellips, slave, slaveorbit, 
00390                     P_ref_checkS, MAXITER, CRITERPOS);
00391     // ______ finally back again from P_ref to slave and master ______
00392     cn M_check;
00393     n_iter     = xyz2orb(M_check,master,masterorbit,P_ref,MAXITER,CRITERTIM);
00394     cn S_check;
00395     n_iter     = xyz2orb(S_check,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);
00396     // ______ Report: ______
00397     INFO.print("Checking zero-Doppler iteration by doing forward and inverse transform:");
00398     DEBUG << "M(x,y,z) = " << M.x     << "; " << M.y     << "; "<< M.z;
00399     DEBUG.print();
00400     DEBUG << "M'(x,y,z)= " << M_check.x     << "; " << M_check.y     << "; "<< M_check.z;
00401     DEBUG.print();
00402     DEBUG << "P(x,y,z) = " << P_ref.x << "; " << P_ref.y << "; "<< P_ref.z;
00403     DEBUG.print();
00404     DEBUG << "P'(x,y,z)= " << P_ref_checkM.x << "; " << P_ref_checkM.y << "; "<< P_ref_checkM.z;
00405     DEBUG.print();
00406     DEBUG << "P''(x,y,z)=" << P_ref_checkS.x << "; " << P_ref_checkS.y << "; "<< P_ref_checkS.z;
00407     DEBUG.print();
00408     DEBUG << "S(x,y,z) = " << S.x     << "; " << S.y     << "; "<< S.z;
00409     DEBUG.print();
00410     DEBUG << "S'(x,y,z)= " << S_check.x     << "; " << S_check.y     << "; "<< S_check.z;
00411     DEBUG.print();
00412     if (abs(M.dist(M_check))>0.001)// assume 1 mm or so?
00413       WARNING.print("check failed for M_check");
00414     if (abs(S.dist(S_check))>0.001)// assume 1 mm or so?
00415       WARNING.print("check failed for S_check");
00416     if (abs(P_ref.dist(P_ref_checkM))>0.001)// assume 1 mm or so?
00417       WARNING.print("check failed for P_ref_checkM");
00418     if (abs(P_ref.dist(P_ref_checkS))>0.001)// assume 1 mm or so?
00419       WARNING.print("check failed for P_ref_checkS");
00420     INFO.print("If no warnings seen then zero-Doppler is within 1 mm accurate.");
00421     }// slave orbit present
00422 
00423   // ______ Tidy up ______
00424   PROGRESS.print("Finished tiepoint.");
00425   DEBUG.reset();
00426   INFO.reset();
00427   } // END tiepoint
00428 

Generated on Fri Apr 22 15:57:49 2005 for Doris by doxygen 1.3.6