Subject: Re:  Meteosat Navigation/Geolocation
From: Instituto de Aeronautica e Espaco Centro Tecnico Aeroespacial <iaeaca@eu.ansp.br>
Date: Wed, 10 Jul 1996 08:04:09 -0300

You can use this routine:



        SUBROUTINE CONVPD(IDIR,PIXEL,LINE,LONG,LAT,IRET)

c Navigation for full resolution METEOSAT images
c Converts image coordinates (line/column) into geographical coordinates
c (latitude/longitude) and vice-versa
c
c       IDIR    I*4   Input     +1 : conversion line/column --> lat/long
c                               -1 : conversion lat/long --> line/column
c       PIXEL   R*4   I or O    pixel column (1-2500)
c       LINE    R*4   I or O    pixel line (1-2500)
c       LONG    R*4   O or I    longitude (degree)
c       LAT     R*4   O or I    latitude (degree)
c       IRET    I*4   Output    =0 if OK
C--------------------------------------------------------------------------

        
        IMPLICIT REAL*8(A-H,O-Z)
        REAL*4 LINE,PIXEL,LONG,LAT
        INTEGER*4 IRET,IDIR
        LOGICAL*1 FIRST
        DATA PI,H / 3.141592653589793 D0, 42164.0 D0 /
        DATA RE / 6378.144 D0 /
        DATA FIRST /.TRUE./

C       initialisations - first call
        IF(FIRST) THEN
        FLAT =1.D0/298.257D0
        RP   =(1.D0-FLAT)*RE
        RAD  =PI/180.0D0
        DEG  =180.0D0/PI
        RE2  =RE*RE
        RP2  =RP*RP
        REP2 =RE2/RP2
        RPE2 =RP2/RE2
        RE2H2=RE2-H*H
        EPSI2=(RE2-RP2)/RE2
        QN   =4.0D-5*PI
        FIRST=.FALSE.
        ENDIF

C       image coord. --> geographic coord.
        IF(IDIR.EQ.1) THEN
        ALPHA=QN*(PIXEL-1250.5)
        BETA =QN*(LINE-1250.5)
        SINA =DSIN(ALPHA)
        SINB =DSIN(BETA)
        COSA =DCOS(ALPHA)
        COSB =DCOS(BETA)
        A =COSB*COSB+REP2*SINB*SINB
        B =COSA*COSB*H
        DISCR=B*B+A*RE2H2
C       check visibility
        IF(DISCR.LT.0.0D0) THEN
          IRET=1
          LONG=0.0
          LAT=0.0
        ELSE
          IRET=0
          RNORM=(B-DSQRT(DISCR))/A
          R1=H-RNORM*COSA*COSB
          R2= -RNORM*SINA*COSB
          R3=  RNORM*SINB
          RL=DSQRT(R1*R1+R2*R2+R3*R3)
          IF(RL.GT.0.0D0) THEN
            RXY=DSQRT(R1*R1+R2*R2)
C           avoid pb near pi/2
            IF(DABS(R1)/RL.GT.1.0D-16) THEN
              RLONG=DATAN(R2/R1)
            ELSE
              RLONG=PI/2.0D0
            ENDIF
            IF(R1.LT.0.0D0) RLONG=RLONG+DSIGN(PI,R2)
            IF(RXY/RL.GT.1.0D-16) THEN
              RLAT=DATAN(R3*REP2/RXY)
            ELSE
              RLAT=DSIGN((PI/2.0D0),R3)
              RLONG=0.0D0
            ENDIF
          ELSE
            RLAT=0.0D0
            RLONG=0.0D0
          ENDIF
          LAT =DEG*RLAT
          LONG=DEG*RLONG
        ENDIF

C       geographic coord. --> image coord.
        ELSE IF(IDIR.EQ.-1) THEN
        RLAT =RAD*LAT
        RLONG=RAD*LONG
        RLAT=DATAN(DTAN(RLAT)*RPE2)
        COSB=DCOS(RLAT)
        SINB=DSIN(RLAT)
        SINA=DSIN(RLONG)
        COSA=DCOS(RLONG)
        RL  =RP/DSQRT(1.0D0-EPSI2*COSB*COSB)
        R1  =H-RL*COSB*COSA
        R2  = -RL*COSB*SINA
        R3  =  RL*SINB
        RNORM=DSQRT(R1*R1+R2*R2+R3*R3)
        ALPHA=DATAN(R2/R1)
        SINA =DSIN(ALPHA)
        COSA =DCOS(ALPHA)
        SINB =R3/RNORM
        BETA =DASIN(SINB)
        COSB =DCOS(BETA)
        A    =COSB*COSB+REP2*SINB*SINB
        B    =COSA*COSB*H
        DISCR=B-A*RNORM
        IF(DISCR.LT.0.0D0) THEN
          IRET=1
          PIXEL=0.
          LINE =0.
        ELSE
          IRET=0
          PIXEL=ALPHA/QN+1250.5D0
          LINE =BETA /QN+1250.5D0
        ENDIF

C       wrong value for IDIR
        ELSE
        IRET=-1
        ENDIF
        RETURN
        END






XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Luiz Augusto Toledo Machado
CTA/IAE/ACA
Sao Jose dos Campos/SP
Brazil
e-mail: iaeaca@eu.ansp.br
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx