conversion.cc

Go to the documentation of this file.

/*
 * Copyright (c) 1999-2003 Bert Kampes
 * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
 *
 * This file is part of Doris, the Delft o-o radar interferometric software.
 *
 * Doris program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Doris is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Publications that contain results produced by the Doris software should
 * contain an acknowledgment. (For example: The interferometric processing
 * was performed using the freely available Doris software package developed
 * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
 * University of Technology, or include a reference to: Bert Kampes and
 * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
 * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
 *
 */
/****************************************************************
 * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/conversion.cc,v $
 * $Revision: 3.14 $
 * $Date: 2005/04/06 15:29:03 $
 * $Author: kampes $
 *
 * Some utilities like deg2rad
 * we need also UTM2WGS etc here to include other DEMs (class?)
 ****************************************************************/


#include "matrixbk.hh"

#include "constants.hh"                 // structs, SOL
#include "readinput.hh"                 // input structs
#include "orbitbk.hh"                    // my orbit class
#include "slcimage.hh"                   // my slc image class
#include "exceptions.hh"                 // my exceptions class

#include "conversion.hh"                // header file, inlines
#include "ioroutines.hh"                // error etc.
#include "utilities.hh"                 // splint etc.
#include "refsystems.hh"                // RADIUS

#include <cmath>                        // sin, sqrt, rint, etc.


/****************************************************************
 *    pol2xyz                                                   *
 *                                                              *
 * Converts     lat,lon,hei(above surface) (sphere)             * 
 *         to   x,y,z (cartesian)                               *
 * uses constants in constants.h (pi, radius)                   *
 *                                                              *
 * lat is positive for North, lon is negative for West.         *
 *                                                              *
 *    Bert Kampes, 21-Dec-1998                                  *
 ****************************************************************/
void pol2xyz(
        cn &xyz,
        real8 phi,
        real8 lon,
        real8 height) //default =0)
  {
  TRACE_FUNCTION("pol2xyz (BK 21-Dec-1998)");
  real8 Rph = RADIUS + height;
  xyz.x = Rph *cos(phi)*cos(lon);
  xyz.y = Rph *cos(phi)*sin(lon);
  xyz.z = Rph *sin(phi);
  }  // END pol2xyz


/****************************************************************
 *    xyz2pol                                                   *
 *                                                              *
 * Converts x,y,z (cartesian)                                   *
 *           to lat,lon,hei(above surface) (sphere)             * 
 * uses constants in constants.h (pi, radius)                   *
 *                                                              *
 * lat is positive for North, lon is negative for West.         *
 *                                                              *
 *    Bert Kampes, 21-Dec-1998                                  *
 ****************************************************************/
void xyz2pol(
        const cn &xyz,
        real8    &phi,
        real8    &lambda,
        real8    &height)
  {
  TRACE_FUNCTION("xyz2pol (BK 21-Dec-1998)");
  real8 r = sqrt(sqr(xyz.x)+sqr(xyz.y));
  phi     = atan2(xyz.z,r);
  lambda  = atan2(xyz.y,xyz.x);
  height  = (r/cos(phi)) - RADIUS;
  }  // END xyz2pol



/****************************************************************
 *    xyz2ell                                                   *
 *                                                              *
 * Converts geocentric cartesian coord. in the XXXX             *
 *  reference frame to geodetic ones, using                     *
 *  method of bowring see globale en locale geodetische systemen*
 * input:                                                       *
 *  - ellipsinfo, xyz, (phi,lam,hei)                            *
 * output:                                                      *
 *  - void (updated lam<-pi,pi> ?, phi<-pi,pi>,hei)             *
 *                                                              *
 *                                                              *
 *    Bert Kampes, 05-Jan-1999                                  *
 ****************************************************************/
void xyz2ell(
        const input_ell &ell,
        const cn        &xyz,
        real8           &phi,
        real8           &lambda,
        real8           &height)
  {
  TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999)");
  real8 r    = sqrt(sqr(xyz.x)+sqr(xyz.y));
  real8 nu   = atan2((xyz.z*ell.a),(r*ell.b));
  real8 sin3 = pow(sin(nu),3);
  real8 cos3 = pow(cos(nu),3);
 
  phi        = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
  lambda     = atan2(xyz.y,xyz.x);
  real8 N    = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
  height     = (r/cos(phi)) - N;
  } // END xyz2ell



/****************************************************************
 *    xyz2ell                                                   *
 *                                                              *
 * Converts geocentric cartesian coord. in the XXXX             *
 *  reference frame to geodetic ones, using                     *
 *  method of bowring see globale en locale geodetische systemen*
 * input:                                                       *
 *  - ellipsinfo, xyz, (phi,lam)                                *
 * output:                                                      *
 *  - void (updated lam<-pi,pi> ?, phi<-pi,pi>)                 *
 *                                                              *
 *                                                              *
 *    Bert Kampes, 05-Jan-1999                                  *
 ****************************************************************/
void xyz2ell(
        const input_ell &ell,
        const cn        &xyz,
        real8           &phi,
        real8           &lambda)
  {
  TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999");
  real8 r    = sqrt(sqr(xyz.x)+sqr(xyz.y));
  real8 nu   = atan2((xyz.z*ell.a),(r*ell.b));
  real8 sin3 = pow(sin(nu),3);
  real8 cos3 = pow(cos(nu),3);
 
  phi        = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
  lambda     = atan2(xyz.y,xyz.x);
  // real8 N    = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
  // height     = (r/cos(phi)) - N;
  } // END xyz2ell



/****************************************************************
 *    ell2xyz                                                   *
 *                                                              *
 * Converts wgs84 ellipsoid cn to geocentric cartesian coord.   *
 * input:                                                       *
 *  - ellipsinfo, phi,lam,hei (x,y,z)                           *
 * output:                                                      *
 *  - void (updated xyz)                                        *
 *                                                              *
 *                                                              *
 *    Bert Kampes, 05-Jan-1999                                  *
 ****************************************************************/
void ell2xyz(
        const input_ell &ell,
        cn    &xyz,
        real8  phi,
        real8  lambda,
        real8  height)
  {
/*
  // --- remove some slow io, since called very often for step ---
  // --- COMPREFDEM #%// Bert Kampes, 24-Sep-2004 ---
  #ifdef __DEBUG
  TRACE_FUNCTION("ell2xyz (BK 05-Jan-1999)");
  #endif
  real8 N   = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
  real8 Nph = N + height;
  #ifdef __DEBUG
  DEBUG << "N = " << N;
  DEBUG.print();
  #endif

  // ______Return these______
  xyz.x =  Nph * cos(phi)  * cos(lambda); 
  xyz.y =  Nph * cos(phi)  * sin(lambda); 
  xyz.z = (Nph - ell.e2*N) * sin(phi);
*/

  const real8 sinphi = sin(phi);// tmp variable
  const real8 N      = ell.a / sqrt(1.0-ell.e2*sqr(sinphi));
  const real8 Nph    = N + height;
  const real8 Nph_cosphi = Nph*cos(phi);// tmp variable
  // ______Return updated xyz______
  xyz.x =  Nph_cosphi      * cos(lambda);
  xyz.y =  Nph_cosphi      * sin(lambda);
  xyz.z = (Nph - ell.e2*N) * sinphi;
  } // END ell2xyz




/****************************************************************
 *    tiepoint                                                  *
 *                                                              *
 * For a given tiepoint, compute the pixel coordinates and      *
 * timing etc.  This can then be used in matlab to correct the  *
 * unwrapped phase by hand.                                     *
 * This routine is called only if tiepoint is given.            *
 *                                                              *
 * Input:                                                       *
 *  lat/lon/hei of a tiepoint                                   *
 * Output:                                                      *
 *  INFO the details                                            *
 *                                                              *
 #%// BK 14-May-2004                                            *
 ****************************************************************/
void tiepoint(
        const input_gen     &generalinput,
        const slcimage      &master,
        const slcimage      &slave,
              orbit         &masterorbit,
              orbit         &slaveorbit,
        const input_ell     &ellips)
  {
  TRACE_FUNCTION("tiepoint (Bert Kampes 14-May-2004)")
  if (abs(generalinput.tiepoint.x) < 0.001 &&
      abs(generalinput.tiepoint.y) < 0.001 &&
      abs(generalinput.tiepoint.z) < 0.001)
    {
    INFO.print("No tiepoint given.");
    return;
    }
  if (masterorbit.npoints()==0)
    {
    INFO.print("Exiting tiepoint, no orbit data master available.");
    return;
    }
  DEBUG.precision(30); // permanent; INFO << setp(13)<<q;// once
  DEBUG.width(14);
  INFO.precision(13); // permanent; INFO << setp(13)<<q;// once
  INFO.width(14);
  const int32 MAXITER   = 10;
  const real8 CRITERPOS = 1e-6;
  const real8 CRITERTIM = 1e-10;
  const real8 m_minpi4cdivlam = (-4.0*PI*SOL)/master.wavelength;
  const real8 s_minpi4cdivlam = (-4.0*PI*SOL)/slave.wavelength;
  DEBUG << "Master wavelength = " << master.wavelength;
  DEBUG.print();
  DEBUG << "Slave wavelength  = " << slave.wavelength;
  DEBUG.print();


  ;// ______ Convert to X,Y,Z ______
  INFO.print("Computing radar coordinates and unwrapped phase of tie point");
  real8 tiepointphi    = deg2rad(generalinput.tiepoint.x);// lat in [dec.degrees]
  real8 tiepointlambda = deg2rad(generalinput.tiepoint.y);// lon in [dec.degrees]
  real8 tiepointheight = generalinput.tiepoint.z;//          height in [m]
  DEBUG << "TIEPOINT: lat/lon/hei [rad/rad/m]: " 
        << tiepointphi << " " << tiepointlambda << " " << tiepointheight;
  DEBUG.print();
  cn tiepointpos;// compute phi/lambda/h-->x,y,z
  ell2xyz(ellips,tiepointpos,tiepointphi,tiepointlambda,tiepointheight);
  INFO  << "TIEPOINT lat/lon/hei --> x,y,z (WGS84): "
        << tiepointphi << ", " << tiepointlambda << ", " << tiepointheight
        << " --> "
        << tiepointpos.x << " " << tiepointpos.y << " " << tiepointpos.z;
  INFO.print();

  // ______ Convert to line/pix in master ______
  real8 tiepointline;// returned
  real8 tiepointpix;// returned
  int32 n_iter = xyz2lp(tiepointline, tiepointpix,
                        master, masterorbit, tiepointpos,
                        MAXITER, CRITERTIM);
  INFO << "TIEPOINT line/pix in master: " << tiepointline << " " << tiepointpix;
  INFO.print();

  // ______ Compute master azimuth/range time of this pixel______
  const real8 m_aztime = master.line2ta(tiepointline);
  const real8 m_trange = master.pix2tr(tiepointpix);
  DEBUG << "TIEPOINT master azimuth time: " << m_aztime;
  DEBUG.print();
  DEBUG << "TIEPOINT master range time:   " << m_trange;
  DEBUG.print();

  // ______ Ellipsoid reference phase for this point master(l,p) ______
  // ______ Compute xyz for slave satellite from P ______
  cn P_ref;// point at H=0 (ellips)
  lp2xyz(tiepointline,tiepointpix,
         ellips,master,masterorbit,P_ref,
         MAXITER,CRITERPOS);
  INFO  << "Pixel of TIEPOINT on ellipsoid: x,y,z (WGS84): "
        << P_ref.x << " " << P_ref.y << " " << P_ref.z;
  INFO.print();

  // ______ Check if slave is there ______
  if (slaveorbit.npoints()==0)
    {
    INFO.print("tiepoint: no orbit data slave available cannot compute more.");
    DEBUG.reset();
    INFO.reset();
    }
  else
    {
    // ______ Compute xyz for slave satellite from P ______
    real8 s_aztime;                             // returned
    real8 s_trange;                             // returned
    xyz2t(s_aztime,s_trange,
          slave, slaveorbit, tiepointpos,MAXITER,CRITERTIM);
    DEBUG << "TIEPOINT slave azimuth time: " << s_aztime;
    DEBUG.print();
    DEBUG << "TIEPOINT slave range time:   " << s_trange;
    DEBUG.print();
    INFO << "TIEPOINT line/pix in slave: " 
         << slave.ta2line(s_aztime) << " "
         << slave.tr2pix (s_trange);
    INFO.print();

    // ___ Phase of this point ___
    real8 phase = m_minpi4cdivlam*m_trange - s_minpi4cdivlam*s_trange;// [rad]
    INFO << "TIEPOINT unwrapped phase (incl. reference phase): "
         << phase;
    INFO.print();
  
    // _____ Times for slave image ______
    real8 s_tazi_ref;                               // returned
    real8 s_trange_ref;                             // returned
    xyz2t(s_tazi_ref,s_trange_ref,slave,
          slaveorbit, P_ref,
          MAXITER,CRITERTIM);

    // ___ report reference phase corrected unwrapped phase of tiepoint ___
    // ___ real8 m_trange_ref = m_trange;// by definition
    // ___ t_corrected = t_m - t_m_ref - (t_s - t_s_ref) = t_s_ref-t_s;
    real8 phase_ref_corrected = s_minpi4cdivlam*(s_trange_ref-s_trange);// [rad]
    INFO << "TIEPOINT unwrapped phase (ellipsoid ref.phase corrected) [rad]: "
         << phase_ref_corrected;
    INFO.print();
    INFO.print("This phase can be used to correct the unwrapped interferogram");
    INFO.print("make sure that you compute the correct pixel position");
    INFO.print("the line/pix here are in original master coordinates.");


    // ______ Check if zero-Doppler iteration is correct by going back ______
    // ------ Go from P_ref to M and to S, then back from M to P_refcheckM
    // ------ and from S to P_refcheckS; Bert Kampes, 06-Apr-2005
    real8 t_az, t_rg;// tmp variables
    // --- from P_ref to orbits ---
    // ______ master: from P_ref to M ____
    n_iter     = xyz2t(t_az,t_rg,master,masterorbit,P_ref,MAXITER,CRITERTIM);// return t
    const cn M = masterorbit.getxyz(t_az);
    // ______ master: from M to P_ref_checkM ____
    cn P_ref_checkM;
    n_iter     = lp2xyz(master.ta2line(t_az), master.tr2pix(t_rg),
                    ellips, master, masterorbit, 
                    P_ref_checkM, MAXITER, CRITERPOS);
    // ______ slave: from P_ref to S ____
    n_iter     = xyz2t(t_az,t_rg,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);// return t
    const cn S = slaveorbit.getxyz(t_az);
    // ______ slave: from S to P_ref_checkS ____
    cn P_ref_checkS;
    n_iter     = lp2xyz(slave.ta2line(t_az), slave.tr2pix(t_rg),
                    ellips, slave, slaveorbit, 
                    P_ref_checkS, MAXITER, CRITERPOS);
    // ______ finally back again from P_ref to slave and master ______
    cn M_check;
    n_iter     = xyz2orb(M_check,master,masterorbit,P_ref,MAXITER,CRITERTIM);
    cn S_check;
    n_iter     = xyz2orb(S_check,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);
    // ______ Report: ______
    INFO.print("Checking zero-Doppler iteration by doing forward and inverse transform:");
    DEBUG << "M(x,y,z) = " << M.x     << "; " << M.y     << "; "<< M.z;
    DEBUG.print();
    DEBUG << "M'(x,y,z)= " << M_check.x     << "; " << M_check.y     << "; "<< M_check.z;
    DEBUG.print();
    DEBUG << "P(x,y,z) = " << P_ref.x << "; " << P_ref.y << "; "<< P_ref.z;
    DEBUG.print();
    DEBUG << "P'(x,y,z)= " << P_ref_checkM.x << "; " << P_ref_checkM.y << "; "<< P_ref_checkM.z;
    DEBUG.print();
    DEBUG << "P''(x,y,z)=" << P_ref_checkS.x << "; " << P_ref_checkS.y << "; "<< P_ref_checkS.z;
    DEBUG.print();
    DEBUG << "S(x,y,z) = " << S.x     << "; " << S.y     << "; "<< S.z;
    DEBUG.print();
    DEBUG << "S'(x,y,z)= " << S_check.x     << "; " << S_check.y     << "; "<< S_check.z;
    DEBUG.print();
    if (abs(M.dist(M_check))>0.001)// assume 1 mm or so?
      WARNING.print("check failed for M_check");
    if (abs(S.dist(S_check))>0.001)// assume 1 mm or so?
      WARNING.print("check failed for S_check");
    if (abs(P_ref.dist(P_ref_checkM))>0.001)// assume 1 mm or so?
      WARNING.print("check failed for P_ref_checkM");
    if (abs(P_ref.dist(P_ref_checkS))>0.001)// assume 1 mm or so?
      WARNING.print("check failed for P_ref_checkS");
    INFO.print("If no warnings seen then zero-Doppler is within 1 mm accurate.");
    }// slave orbit present

  // ______ Tidy up ______
  PROGRESS.print("Finished tiepoint.");
  DEBUG.reset();
  INFO.reset();
  } // END tiepoint

---