McIDAS Programmer's Manual
Version 2003

[Search Manual] [Table of Contents] [Go to Previous] [Go to Next]


McIDAS navigation

All geophysical data must have a location to be meaningful. McIDAS-X stores geophysical data in the form closest to its native spatial structure; for example, remotely sensed data is stored in the image coordinates of the scanner line and element, while NWP model output is stored as model grid points. Since earth coordinates (for example, latitude, longitude and height above mean sea level) are easiest to interpret, McIDAS-X uses a process called navigation to convert a dataset's native coordinate system to earth coordinates and back again.

Navigating datasets implies the ability to transform any dataset into the spatial coordinates of any other dataset. McIDAS-X' ability to integrate and display a variety of datasets is made possible by its navigation subsystems.

This section describes:

 

For more information about McIDAS-X coordinate systems, see Chapter 2, Learning the Basics. For additional information about map projections, see Map Projections -- A Working Manual, USGS Professional Paper 1395, U.S. Government Printing Office, Washington, 1987.

Basic concepts

To select and use the best McIDAS-X navigation subsystem for a given task, you need to understand the nature of navigation transforms. Navigation systems and transforms are described below, along with some navigation-specific terminology used throughout this section.

Navigation transforms

A navigation transform is a set of equations for converting a dataset's image or native coordinates to and from earth coordinates. Each dataset has its own navigation transform. However, navigation transforms can be grouped into classes or types. Within a type, the differences between transforms are quantitative only. Each type has a set of variables, or parameters, whose values uniquely specify an instance of that type. For example, the parameters for a tangent cone projection include the standard latitude and longitude, the scale factor, and the location of the pole relative to the projection surface's origin.

All navigated McIDAS-X data has an associated type and a full set of parameters defining an instance.

Navigation subsystems

McIDAS-X uses three navigation subsystems.

Grid navigation

Grid navigation converts earth coordinates to and from the row and column locations in the grid structure. This subsystem is separate from image and frame navigation, and supports fewer projections These projection types are implemented in a single set of functions. To add other types, you would have to modify the functions and recompile all grid applications. Since only a single instance can be in use at one time, you must repeatedly reinitialize the subsystem to use multiple projections within an application.

For example, to compare two grids, A and B, you initialize grid navigation for A, compute the earth location of a grid point a in A, reinitialize the subsystem using B's parameters, and transform the earth location of a to a grid in B. This process is then repeated for every grid point being compared. Most applications that process or use grids use grid navigation.

 

For more information about grid projection formats, see the Grid data section earlier in this chapter.

Image navigation

Image navigation converts image coordinates to and from earth coordinates. Each image navigation type is implemented in its own source file or module, and each navigable image dataset has an associated navigation block that identifies the type and contains a complete set of parameters needed to define an instance of that type.

Each module implements the same function names and calling sequence. Since McIDAS-X uses the type identifier to select the correct module, you can add new types without modifying existing modules. In addition, up to three instances of navigation, of the same or different types, can be active simultaneously in a single application through different slots. This requires a one- or two-step initialization process before use and is typically used by applications that work directly with image datasets.

Frame navigation

Frame navigation converts frame coordinates to and from earth coordinates. Frame navigation uses the image navigation subsystem but through a simpler API. It is sufficient for most applications that just add information to an existing display.

Terminology

This section contains some navigation-specific terms that may be unfamiliar to you. They are defined below.

Using the navigation APIs

The McIDAS-X library provides Application Program Interfaces for geographic calculations as well as grid, image and frame navigation. These APIs are described below.

Geographic calculation API

Many calculations just require a change from one form of earth or planetary coordinate to another; for example, changing from a terrestrial (earth-relative) to celestial (star-relative) coordinate system, or converting a position from a Cartesian to spherical representation. The table below lists the functions for performing these tasks.

Functions Description

m0distnc

calculates the great-circle distance on a sphere

geolat

performs geodetic/geocentric latitude conversions

nllxyz, nxyzll

perform spherical/Cartesian conversions

raerac, racrae

perform terrestrial/celestial longitude conversions

solarp

calculates solar position

cartll, llcart, llobl, llopt

perform generic spherical/Cartesian conversions

Grid navigation API

The table below shows the grid navigation API. You will use the grddef function with the header to initialize the subsystem, then call ijll and llij to transform grid coordinates to and from earth coordinates, respectively.

Function Description

grddef

initializes the grid transformation

ijll

converts grid coordinates to earth coordinates

llij

converts earth coordinates to grid coordinates

No API functions currently exist for creating grid headers for a projection. When creating McIDAS-X grid files, you must determine the projection to use and then insert the projection type and parameters into the grid header. The grid header format is described in the GRIDnnnn data structure in Chapter 6.

Frame navigation API

Use the two functions below to determine the earth coordinates of a pixel or the pixel location of a point in earth coordinates.

Function Description

illtv

converts earth coordinates to image and frame coordinates

itvll

converts frame coordinates to image and earth coordinates

The only prerequisite for these functions is that the frame has navigation associated with it. The McIDAS-X IMGDISP command creates frame navigation when it displays an image dataset. The McIDAS-X MAP, PTDISP and GRDDISP commands, which can define a map base before generating a display, also establish frame navigation.

The two examples below assume that frame navigation information exists. The first code fragment uses the illtv function to convert an earth location (zlat,zlon) to a frame location (itvlin,itvele), as is done in the McIDAS-X PC command.

zlat = 43.13
zlon = 89.35

if(illtv(iframe,gnum,zlat,zlon,ilin,iele,itvlin,itvele).ne.0)then
   call edest('unable to navigate point',0) 
   return
endif

The second code fragment uses itvll to determine the earth location of the cursor, as is done in the McIDAS-X E command.

     if( mcmoubtn( 0, left, right, itvlin, itvele ) ) then
C        // handle cursor-positioning error here
     end if
     frame = mcgetimageframenumber()
     rc    = itvll(frame,itvlin,itvele,ilin,iele,rlat,rlon,iscene)

Image navigation API

Use the appropriate image navigation API functions from the table below if your application must perform one of these tasks:

Function Description

nvprep

dynamically loads a navigation module for a specified slot, given a set of parameters

nvset

initializes navigation for a specified slot from a specified source

nvxeas

converts earth coordinates to image coordinates; due to dynamic linking, this function is called nv1eas, nv2eas or nv3eas from applications, depending on the navigation slot used

nvxini

initializes the navigation module; due to dynamic linking, this function is called nv1ini, nv2ini or nv3ini from applications, depending on the navigation slot used

nvxopt

performs special navigation services and is module-dependent; due to dynamic linking, this function is called nv1opt, nv2opt or nv3opt from applications, depending on the navigation slot used

nvxsae

converts image coordinates to earth coordinates; due to dynamic linking, this function is called nv1sae, nv2sae or nv3sae from applications, depending on the navigation slot used

To use image navigation services, the application must do the following:

The example below, from remap.pgm, illustrates the initialization and use of two instances (slots) of navigation at the same time.

The lines below assume that navigation blocks were already obtained from the source (snav) and destination (dnav) areas. To link the appropriate navigation module to a slot and initialize it, use nvprep. The lines below initialize slot 1 for the source and slot 2 for the destination. The nv1ini and nv2ini calls set the earth coordinate form to Cartesian for both slots.

c --- check for initialization
      if( init.eq.0 ) then
c ------ initialize the source navigation
         if( nvprep( 1, snav(1) ).ne.0 ) then
        navinit = -1
        return
         endif
c ------ initialize the destination navigation
         if( nvprep( 2, dnav(1) ).ne.0 ) then
        navinit = -1
        return
         endif
     init = 1
      endif
c ------ initialize the navigations
c ------ trans can be 'LL  ' for latitude longitude
c ------ or 'XYZ ' for cartesian 
      itrans = lit( trans(1:4) )
      call nv1ini(2, itrans)
     call nv2ini(2, itrans)

The lines below map image coordinate (dline,delement) in the destination image to image coordinate (sline,selement) in the source image by converting from destination image coordinates to earth coordinates, then from earth coordinates to source image coordinates.

c --- convert destination line/element to latitude/longitude
      status = nv2sae( dline, delement, 0., dlat, dlon, dz )
      if( status.ne.0 ) then
     crossnav = -1
     return
      endif
c --- convert latitude/longitude to source line/element
      status = nv1eas( dlat, dlon, dz, sline, selement, dz)
      if( status.ne.0 ) then
     crossnav = -2
     return
      endif

When the navigation block source is a local McIDAS-X area format image file or a frame whose number is known and you only intend to use slot 1, you can use nvset. It reads the navigation block out of an area or frame and calls nvprep.

nvset ('AREA', area_num)
  -or-
nvset ('FRAME', frame_num)

In addition to the commonly used conversions between earth and image coordinates available through nvxsae and nvxeas, other special services are available through the nvxopt function. Some of these services and the modules that support them are listed in the table below.

Special services Modules that support them

oblate radius

nvxmerc, nvxps, nvxrect

orbital period

nvxdmsp

parallax correction

nvxgoes, nvxgvar, nvxmsat

satellite location/subpoint

nvxdmsp, nvxgoes, nvxgraf, nvxgvar, nvxlamb, nvxmerc, nvxmsat, nvxps, nvxradr, nvxrect, nvxsin, nvxtiro

view angles (satellite, solar, etc.)

nvxgoes, nvxgvar

The sample code below shows the use of nvxopt. When using nvxopt, you need to be aware that the input/output can be different for the different satellite navigations, and you'll have to check the nvx sat.dlm file to check what is appropriate for input/output.

      real         xin(6)
      real         xout(6)

      ... 

C --- initialize slot 1 navigation (navblk is assumed to be
C --- filled with valid parameters)
      if( nvprep( 1, navblk ) .ne.0 ) then
C         (error handler)
      end if

C --- call nvxopt(), and extract the outputs
      if( nv1opt( lit('SPOS'), xin, xout ) .ne.0 ) then
C         (error handler)
      end if

      lat = xout(1)
      lon = xout(2)

Implementing an image navigation module

You can extend image navigation to include a new type and make it available to all McIDAS-X commands without changing any applications or existing navigation modules. This section describes the process for implementing a tangent-cone map projection. The complete source listing is in the next section titled Example navigation module code.

Requirements

The new module must adhere to strict guidelines governing the following:

Navigation block contents

The navigation block must consist of an array of no more than 640 Fortran integers that are sufficient to uniquely specify a particular instance of the new navigation type. The first word of the block must be four ASCII characters that uniquely identify the type, stored in an integer variable. The block contents are obvious once the algorithm is fully specified and understood.

Module name

The navigation module name must be nvxtype.dlm, where type is the same four characters specified in the first word of the block. The nvprep function uses the type from the block to link the application with the correct module.

Functions

The navigation module must implement the functions nvxini, nvxsae, nvxeas and nvxopt, with calling sequences following the documentation blocks in the tangent cone example. In particular, nvxini must support modes 1 (initialize) and 2 (set earth coordinate mode 'LL' (Latitude-longitude) or 'XYZ' (Cartesian)). The nvxopt function must provide, at minimum, an 'SPOS' (subpoint) service.

Parameters

Between calls, the navigation module must remember the navigation parameters uniquely defining an instance. This is typically accomplished by storing these parameters, and quantities derived from them, in a named common block that is shared among all four API functions in the module.

Invalid inputs

Satellites and map projections cannot physically represent all earth locations. The limits vary from type to type and even from instance to instance. Applications, therefore, cannot be expected to validate inputs to navigation modules. The new module must be able to recognize inputs it cannot handle and return an error status, typically a negative integer.

Algorithms

A navigation algorithm is a set of equations for mapping McIDAS-X image coordinates to earth coordinates and vice versa. For satellite navigation types, this typically involves a lengthy expression in vector notation. Map projections are simpler. Some manipulation and derivation of additional relationships is often necessary to extend a published relationship into a full algorithm.

The example presented here will show the development of a northern hemisphere tangent cone projection algorithm. The geometry of the tangent cone is shown in Figure 5-6 below. The cone is tangent to the planet at one parallel of latitude and the planetary axis pierces the apex of the cone. A line between a point on the Earth's surface and the pole opposite the cone apex maps that point onto the cone. The imaginary cone is then cut along a meridian opposite an arbitrary standard longitude and flattened. Each point on the Earth's surface in colatitude ψ and longitude λ corresponds to a point (R,θ) in polar coordinates on the flattened cone.

Figure 5-6. Tangent cone projection geometry.


 


A published form of the equations defining the location (R,θ) on the flattened cone in terms of an earth location (ψ,λ) is:

(1)

(2)

Also provided is an equation for the ratio of the distance between two closely spaced points on the cone to the distance between the corresponding points on the Earth's surface.
(3)
As you can see from (1) and (2), the planetary radius, a the standard colatitude ψ0, and the standard longitude λ0 define an instance of the tangent cone. In equation 3, m is the map scale expressed as a unitless ratio of projection plane distance to map distance.

This is not yet a complete algorithm ready for McIDAS-X implementation. The inverse relationship capable of converting earth location to projection coordinate (R,θ) must be derived, the earth location (λ,ψ) related to its equivalent in McIDAS-X (LAT, LON), and the projection coordinate (R,θ) related to McIDAS-X image coordinates. The inverse is:

 

(4)

 

and
(5)
The relationships between McIDAS-X planetary coordinates and (λ,ψ) are:
(6)
and
(7)
Those between projection and image coordinates are:

(8)

and
(9)
 

Adding (6) - (9) to make a complete McIDAS-X navigation algorithm adds additional parameters to make the complete set ψ0, λ0, m, L0, and E0, where the latter specify the scale in km per pixel and the location of the pole of the projection surface in image coordinates.

The last step before implementation is to examine the completed algorithm for limits and singularities; these will serve as the basis for input validation. The limits on the inputs are summarized below.

Parameters  

0 < ψ0 < π/2

-π < λ0 ≤ π

standard colatitude is confined to the northern hemisphere

standard longitude must be a legal value

Earth coordinates  
0 ≤ ψ < π

-π < λ ≤ π

all colatitudes are navigable except the south pole; note that map scale M(ψ) is undefined at the north pole,

however all longitudes are navigable

Image coordinates
R ≥ 0

π cos(ψ0) < θ ≤ π cos( ψ0)

no valid earth location maps to -R

area in the split region of the flattened cone isn't navigable; exclude it as input

Implementation

Implementing the navigation module consists of three steps:

  1. Specifying the navigation block format
  2. Creating or modifying the routines and applications that will create instances of the new type
  3. creating the navigation module for the new type

Each of these steps is discussed below.

Defining the navigation block

The contents of the navigation block follow directly from the parameters provided in the previous Algorithms section. The five parameters are:

A sixth parameter, planetary radius, could also be added but is not included in the examples. The parameters can be represented in different ways, such as standard colatitude in place of latitude. However, since the applications will create navigation blocks, you should use units and quantities convenient for McIDAS-X, as shown in the lines of code below. This code is taken from the documentation block at the beginning of nvxini and is shown in lines 0017 through 0022 in the complete source listing that follows.

*$   param( 1) = 'TANC'
*$   param( 2) = image line of pole   *10000
*$   param( 3) = image element of pole*10000
*$   param( 4) = km per pixel         *10000
*$   param( 6) = standard latitude    *10000
*$   param( 7) = standard longitude   *10000

All instances of 'TANC' navigation must be specified with these parameters.

Creating navigation blocks in applications

Applications that create navigation blocks must provide values for a complete set of parameters. Sometimes they must be derived from user inputs using the navigation transformation itself. Below is a code fragment from an application (not the navigation module itself) that creates a tangent-cone navigation block. It assumes that the center latitude and longitude (clat and clon), the standard latitude and longitude (slat and slon), and the scale of standard latitude (sscale), in km per pixel, are already provided. First, the earth locations are converted from McIDAS-X to projection form.

    psi   = PI/2.d0 - D2R*clat
    psi_0 = PI/2.d0 - D2R*slat
    lam   =         - D2R*clon
    lam_0 =         - D2R*slon

Next, the radius of the specified image coordinate center point (relative to the projection center or the pole) is computed.

    radius = A * tan(psi_0) *
$       (tan(psi/2.d0)/tan(psi_0/2.d0))**cos(psi_0)
    theta  = cos(psi_0)*(lam-lam_0)

Then the pole's image coordinate is computed using simple trigonometry.

    call status = mcfsize(mcgetimageframenumber, nlins, neles)
    lin_c  = dble(nlins) / 2.d0
    ele_c  = dble(neles) / 2.d0


    lin_0  = lin_c - radius*cos(theta)/sscale
    ele_0  = ele_c - radius*sin(theta)/sscale

Once the parameters are computed, the navigation block can be filled and written either to the image dataset:

C   // Create and insert the navigation block

    nvblk( 1) = lit('TANC') ! Tangent Cone nav type
    nvblk( 2) = nint(10000*lin_0)   ! image line of pole
    nvblk( 3) = nint(10000*ele_0)   ! image element of pole
    nvblk( 4) = nint(10000*sscale)  ! scale km/pixel at slat
    nvblk( 5) = nint(10000*slat)    ! standard latitude
    nvblk( 6) = nint(10000*slon)    ! standard longitude

    status = mcaput(dataset, nsort, sort, adir, nvblk, calblk)

or a McIDAS-X area format image file:

    call araput(aranum,4*DIRSIZ,4*NVWDS,nvblk)

Creating the navigation module

If the algorithm and navigation block are complete and unambiguous, implementing a map projection as a McIDAS-X navigation module is generally straightforward. Your major design decision is what to store in the common block. At minimum, you should include the navigation parameters or quantities derived from them without loss of information. This is necessary for the module to remember the characteristics of the instance.

Because navigation modules are called many times, you should precompute additional quantities in the initialization routine nvxini and store them in common, especially if doing so eliminates trigonometric function calls in the other functions. Thus, the variables Coscl, Tancl, Tancl2, and Mxtheta are included in the common block below. This code is located in lines 0099 through 0120 of the source listing that follows.

C     // Common block contents: precomputed intermediate values

      real      Coscl       ! cosine(Colat0)
      real      Tancl       ! tangent(Colat0)
      real      Tancl2      ! tangent(Colat0/2)
      real      Mxtheta     ! limit of angle from std. lon
C               ! on projection surface

C     // Common block contents: constants

      real      D2R     ! degrees to radians factor
      real      Pi
      real      Badreal     ! returned when navigation
C                       ! cannot be done
      real      Erad        ! Earth radius
      logical       Init        ! initialized flag
      logical       Latlon      ! .TRUE. for lat/lon I/O


      COMMON/TANC/ Lin0, Ele0, Scale, Lon0, Colat0,
     &  Coscl, Tancl, Tancl2, Mxtheta,
     &  D2R, Pi, Badreal, Erad, Init, Latlon


Also included in the block are physical constants and two flags; Init is set when nvxini runs successfully, and Latlon is set if the earth coordinate mode is latitude and longitude rather than Cartesian.

The overall structure of nvxini is a large if-block to handle the two options. Option 1 (lines 0144 - 0220 in the source listing) initializes the instance by range checking the navigation parameters and precomputing the quantities to retain in the common block. Option 2 (lines 0222 - 0233) sets the earth coordinate mode based on word param(1).

The forward (image to Earth) navigation routine's one subtlety is input validation. Points within the split of the projection plane, which is shown in Figure 5-6, are not navigable.

Recognizing such points is possible only after the conversion from McIDAS-X image coordinates to projection coordinates is complete. The code below is from lines 0370 through 0384.

C  // Compute radius and bearing from pole

   dx   = Scale*(lin-Lin0)
   dy   = Scale*(ele-Ele0)

   radius   = sqrt(dx*dx+dy*dy)
   theta    = atan2(dy,dx)


C  // Apply range checking on theta to determine if point is navigable

   if ( theta.le.-Mxtheta .or. theta.gt.Mxtheta ) then
      NVXSAE = -1
      return
   end if

When range checking is done, the projection is inverted and the colatitude and longitude (in radians, east positive) is converted to latitude and longitude (in degrees, west positive); see lines 0389 - 0401. A call to nllxyz to generate Cartesian coordinate output is made if the Latlon flag was cleared by an earlier call to nvxini (line 0404). Inverse (earth-to-image) navigation in nvxeas follows a similar pattern, except the earth-coordinate option must be handled at the beginning (line 0544).

The nvxopt function performs special services. Since the contents of the argument vector depend on the option selected, the routine consists of a large if-block (lines 0716 - 0742) with one branch per recognized option and the necessary input validation done within each branch. As in nvxini, the code must be able to recognize options that it doesn't know and return an appropriate error status. All core McIDAS-X navigation modules have an SPOS option that returns the subpoint at a given time (for satellites) or the latitude and longitude of the center (maps). The example that follows implements only an SCAL to return the map scale factor at any latitude.

1. Saucier, W. J. 1983: Principles of meteorological analysis. Dover Publications, Inc., New York. 438 pp.

Example navigation module code

0001    C   THIS IS SSEC PROPRIETARY SOFTWARE - ITS USE IS RESTRICTED.
0002    
0003    C *** McIDAS Revision History ***
0004    C *** McIDAS Revision History ***
0005    
0006    *$ Name:
0007    *$  nvxini - Initialize navigation for tangent cone projection
0008    *$
0009    *$ Interface:
0010    *$  integer function
0011    *$  nvxint(integer option, integer param(*))
0012    *$
0013    *$ Input:
0014    *$  option  - 1 to set or change projection parameters
0015    *$  option  - 2 set output option
0016    *$  param   - For option 1:
0017    *$        param( 1) = 'TANC'
0018    *$        param( 2) = image line    of pole*10000
0019    *$        param( 3) = image element of pole*10000
0020    *$        param( 4) = km per pixel         *10000
0021    *$        param( 6) = standard latitude    *10000
0022    *$        param( 7) = standard longitude   *10000
0023    *$          for option 2:
0024    *$          param( 1) = 'LL' or 'XYZ'
0025    *$
0026    *$ Input and Output:
0027    *$  none
0028    *$
0029    *$ Output:
0030    *$  none
0031    *$
0032    *$ Return values:
0033    *$   0  - success
0034    *$  -3  - invalid or inconsistent navigation parameters
0035    *$  -4  - invalid navigation parameter type
0036    *$  -5  - invalid nvxini() option
0037    *$
0038    *$ Remarks:
0039    *$  Latitudes and longitudes are in degrees, West positive.
0040    *$  Projection parameters must be in the following ranges:
0041    *$      0. <  standard latitude  <  90.
0042    *$      -180. <= standard longitude < 180.
0043    *$      0. <  scale
0044    *$  Accuracy may suffer near the standard latitude limits.
0045    *$
0046    *$  The projection algorithm is adapted from that in
0047    *$  Saucier, W. J. 1989: Principles of meteorological analysis.
0048    *$  Dover Publications, Inc. 433 pp.
0049    *$
0050    *$ Categories: 
0051    *$  navigation  
0052    
0053    C     // CODING CONVENTION note: function declarations and common
0054    C     // block declarations are all capitalized to be recognizable
0055    C     // to script 'convdlm;' this is necessary for a correct build
0056    C     // in MCIDAS-X. For the same reason, avoid referring to
0057    C     // function or common block names in uppercase elsewhere
0058    
0059          INTEGER FUNCTION NVXINI(option,param)
0060    
0061    
0062          implicit      NONE
0063    
0064    
0065    C     // Interface variables (formal arguments)
0066    
0067          integer           option      ! initialization option
0068          integer           param(*)        ! navigation parameters or
0069    C                   ! output coordinate type
0070    
0071    C     // Local variable definitions
0072    
0073          character*4           navtyp      ! codicil type
0074          character*4           outcoord        ! output coordinate type
0075          real          lat0        ! standard latitude
0076          character*80          cbuf        ! text output buffer
0077    
0078    
0079    
0080    C     /////////////////////////////////////////////////////////
0081    C     Common block variables and declaration.
0082    
0083    C     ALL CODE BETWEEN THE '//////' SEPARATORS MUST BE 
0084    C     DUPLICATED EXACTLY IN EACH NAVIGATION ROUTINE
0085    
0086    C     (A more maintenance-safe version would use ENTRY points
0087    C     rather than separate functions for the navigation API
0088    C     but entry points cannot be processed by 'convdlm.')
0089    
0090    C     // Common block contents: projection parameters
0091    
0092          real          Lin0        ! image line of pole
0093          real          Ele0        ! image element of pole
0094          real          Scale       ! km per unit image coordinate
0095    C                   ! (pixel)
0096          real          Lon0        ! standard longitude
0097          real          Colat0      ! standard colatitude
0098    
0099    C     // Common block contents: precomputed intermediate values
0100    
0101          real          Coscl       ! cosine(Colat0)
0102          real          Tancl       ! tangent(Colat0)
0103          real          Tancl2      ! tangent(Colat0/2)
0104          real          Mxtheta     ! limit of angle from std. lon
0105    C                   ! on projection surface
0106    
0107    C     // Common block contents: constants
0108    
0109          real          D2R     ! degrees to radians factor
0110          real          Pi
0111          real          Badreal     ! returned when navigation
0112    C                   ! cannot be done
0113          real          Erad        ! Earth radius
0114          logical           Init        ! initialized flag
0115          logical           Latlon      ! .TRUE. for lat/lon I/O
0116    
0117    
0118          COMMON/TANC/ Lin0, Ele0, Scale, Lon0, Colat0,
0119         &  Coscl, Tancl, Tancl2, Mxtheta,
0120         &  D2R, Pi, Badreal, Erad, Init, Latlon
0121          
0122    C     End of common block variables and declaration.
0123    C     /////////////////////////////////////////////////////////
0124    
0125    
0126    
0127    C     // Begin initialization process by setting values of constants.
0128    
0129          Erad      = 6370.         ! This value of Erad is ok 
0130    C                   ! for low precision nav
0131    C                   ! where a spherical Earth is
0132    C                   ! adequate (Saucier, p. 32)
0133          Pi        = acos(-1.)
0134          D2R       = Pi / 180.
0135          Badreal       = -1.E10            ! obvious unreasonable value
0136    C                   ! for nav transform result
0137    
0138    C     // Process initialization options. Only one option, initialize
0139    C     // navigation parameters, is supported in this demo version,
0140    C     // but a 'hook' is left for an additional option to set the
0141    C     // output coordinate to something other than lat/lon
0142    
0143          if( option.eq.1 ) then
0144              
0145             call DDEST('nvxini(tanc) option=1',0)
0146    
0147             call movwc(param(1),navtyp)
0148             if( navtyp.eq.'TANC') then
0149    
0150    C           // Unpack tangent cone projection parameters
0151    
0152                Lin0            = param(2) / 10000.
0153                Ele0            = param(3) / 10000.
0154                Scale           = param(4) / 10000.
0155                lat0            = param(5) / 10000.
0156                Lon0            = param(6) / 10000.
0157    
0158                write(cbuf,'('' nvxini: lat0, Lon0 '',2F12.4)')
0159         *        lat0, Lon0
0160                call DDEST(cbuf,0)
0161    
0162    C           // apply range checking
0163    
0164                if(Scale.le.0. ) then
0165                   call DDEST('nvxini(tanc) scale is negative',0)
0166                   Init         = .FALSE.
0167                   NVXINI           = -3
0168                   return
0169                end if
0170    
0171                if(lat0.le.0. .or. lat0.ge.90. ) then
0172                   call DDEST('nvxini(tanc) standard lat out of range',0)
0173                   Init         = .FALSE.
0174                   NVXINI   = -3
0175                   return
0176                end if
0177    
0178                if(Lon0.le.-180. .or. Lon0.gt.180. ) then
0179                   call DDEST('nvxini(tanc) standard lon out of range',0)
0180                   Init         = .FALSE.
0181                   NVXINI   = -3
0182                   return
0183                end if
0184    
0185    C           // convert degrees to radians and latitude to colatitude.
0186    C           // Account for McIDAS longitude convention
0187    
0188                Lon0            = -Lon0 * D2R
0189                Colat0          =  Pi/2. - D2R*lat0
0190    
0191                write(cbuf,'('' nvxini:  Colat0, Lon_0 '',2F12.4)')
0192         *        Colat0, Lon0
0193                call DDEST(cbuf,0)
0194    
0195    C           // Compute intermediate quantities
0196    
0197                Coscl           = cos(Colat0)
0198                Tancl           = tan(Colat0)
0199                Tancl2          = tan(Colat0/2.)
0200                Mxtheta         = Pi*Coscl
0201    
0202                write(cbuf,'('' nvxini: Coscl,   Tancl'', 2F7.4)')
0203         *        Coscl, Tancl
0204                call DDEST(cbuf,0)
0205                write(cbuf,'('' nvxini: Tancl2, Mxtheta '', 2F7.4)')
0206         *        tancl2, Mxtheta
0207                call DDEST(cbuf,0)
0208    
0209            Latlon      = .TRUE.
0210    
0211    
0212             else           ! option=1 but type not 'TANC'
0213    
0214                call DDEST('nvxini(tanc) parameter type bad',0)
0215                Init            = .FALSE.
0216                NVXINI          = -4
0217                return
0218    
0219             end if
0220    
0221          else if ( option .eq. 2) then
0222    
0223         call movwc(param(1),outcoord)
0224         if( outcoord.eq.'LL' ) then
0225             Latlon     = .TRUE.
0226         else if( outcoord.eq.'XYZ') then
0227             Latlon     = .FALSE.
0228         else
0229                 call DDEST('option=2 coord '//outcoord//' not supported',0)
0230                 Init           = .FALSE.
0231                 NVXINI         = -5
0232         end if
0233    
0234          else          ! option not 1 or 2
0235    
0236             call DDEST('nvxini(tanc) unrecognized output option ',option)
0237             NVXINI     =  -4
0238             return
0239    
0240          end if
0241    
0242          NVXINI        = 0
0243          Init      = .TRUE.
0244    
0245          return
0246          end
0247        
0248    
0249    
0250    
0251    *$ Name:
0252    *$  nvxsae - Compute earth coordinates from image coordinates
0253    *$
0254    *$ Interface:
0255    *$  integer function
0256    *$  nvxsae( real lin, real ele, real dummy,
0257    *$              real e1,  real  e2, real    e3 )
0258    *$
0259    *$ Input:
0260    *$  lin - image line
0261    *$  ele - image element
0262    *$  dummy   - (unused)
0263    *$
0264    *$ Input and Output:
0265    *$  none
0266    *$
0267    *$ Output:
0268    *$  e1  - latitude  or x
0269    *$  e2  - longitude or y
0270    *$  e3  - height    or z
0271    *$
0272    *$ Return values:
0273    *$   0  - success
0274    *$  -1  - input data physically valid but not navigable
0275    *$        given the specified projection
0276    *$  -6  - module not initialized
0277    *$
0278    *$ Remarks:
0279    *$  The navigation module must first be initialized with
0280    *$  a call to nvxini(). The output form (lat,lon) or (x,y,z)
0281    *$  depends on the last call to nvxini() with option 2.
0282    *$
0283    *$ Categories: 
0284    *$  navigation  
0285    
0286          INTEGER FUNCTION NVXSAE( lin, ele, dummy, e1, e2, e3 )
0287    
0288          implicit      NONE
0289    
0290    C     // Interface variables (formal arguments)
0291    
0292          real          lin     ! image line to navigate
0293          real          ele     ! image element to navigate
0294          real          dummy       ! (unused argument)
0295          real          e1      ! Earth coordinate 1
0296          real          e2      ! Earth coordinate 2
0297          real          e3      ! Earth coordinate 3
0298    
0299    C     // Local variables
0300    
0301          real          lat     ! latitude  (McIDAS convention)
0302          real          lon     ! longitude (McIDAS convention)
0303          real          hgt     ! height
0304          real          dx      ! zonal displacement from pole
0305    C                   ! on projection surface
0306          real          dy      ! meridional displacement from pole
0307          real          radius      ! distance from pole on projection
0308          real          theta       ! angle from standard longitude on
0309    C                   ! projection surface
0310          real          colat       ! colatitude of navigated point
0311    
0312    
0313    C     ///////////////////////////////////////////////////////
0314    C     Common block variables and declaration.
0315    
0316    C     ALL CODE BETWEEN THE '//////' SEPARATORS MUST BE 
0317    C     DUPLICATED EXACTLY IN EACH NAVIGATION ROUTINE
0318    
0319    C     (A more maintenance-safe version would use ENTRY points
0320    C     rather than separate functions for the navigation API
0321    C     but entry points cannot be processed by 'convdlm.')
0322    
0323    C     // Common block contents: projection parameters
0324    
0325          real          Lin0        ! image line of pole
0326          real          Ele0        ! image element of pole
0327          real          Scale       ! km per unit image coordinate
0328    C                   ! (pixel)
0329          real          Lon0        ! standard longitude
0330          real          Colat0      ! standard colatitude
0331    
0332    C     // Common block contents: pre-computed intermediate values
0333    
0334          real          Coscl       ! cosine(Colat0)
0335          real          Tancl       ! tangent(Colat0)
0336          real          Tancl2      ! tangent(Colat0/2)
0337          real          Mxtheta     ! limit of angle from std. lon
0338    C                   ! on projection surface
0339    
0340    C     // Common block contents: constants
0341    
0342          real          D2R     ! degrees to radians factor
0343          real          Pi
0344          real          Badreal     ! returned when navigation
0345    C                   ! cannot be done
0346          real          Erad        ! Earth radius
0347          logical           Init        ! initialized flag
0348          logical           Latlon      ! .TRUE. for lat/lon I/O
0349    
0350          COMMON/TANC/ Lin0, Ele0, Scale, Lon0, Colat0,
0351         &  Coscl, Tancl, Tancl2, Mxtheta,
0352         &  D2R, Pi, Badreal, Erad, Init, Latlon
0353          
0354    C     End of common block variables and declaration.
0355    C     ///////////////////////////////////////////////////
0356    
0357          e1        = Badreal
0358          e2        = Badreal
0359          e3        = Badreal
0360    
0361    
0362    C     // verify initialized module
0363    
0364          if(.not.Init) then
0365              NVXSAE = -6
0366              return
0367          end if
0368    
0369    
0370    C     // Compute radius and bearing from pole
0371    
0372          dx        = Scale*(lin-Lin0)
0373          dy        = Scale*(ele-Ele0)
0374    
0375          radius        = sqrt(dx*dx+dy*dy)
0376          theta     = atan2(dy,dx)
0377    
0378    
0379    C     // Apply range checking on theta to determine if point is navigable
0380    
0381          if ( theta.le.-Mxtheta .or. theta.gt.Mxtheta ) then
0382             NVXSAE = -1
0383             return
0384          end if
0385    
0386    C     // Forward navigation: compute longitude and colatitude
0387    C     // from radius and theta
0388    
0389          lon       = Lon0 + theta/Coscl
0390          if(lon.le.-Pi) lon = lon + 2.d0*Pi
0391          if(lon.gt. Pi) lon = lon - 2.d0*Pi
0392    
0393          colat     = 2.*atan( Tancl2 * (radius/(Erad*Tancl))**(1./Coscl))
0394    
0395    C     // Rescale to McIDAS convention (degrees, West positive).
0396    C     // Apply conversion to Cartesian coordinates if 'XYZ' set
0397    C     // as output form. Set return code for success.
0398    
0399          lon       = -lon/D2R
0400          lat       = 90. - colat/D2R
0401          hgt       = 0.
0402    
0403          if(.not.Latlon) then
0404              call nllxyz(lat,lon,e1,e2,e3)
0405          else
0406          e1 = lat
0407          e2 = lon
0408          e3 = 0.
0409          end if
0410    
0411          NVXSAE        = 0
0412    
0413          return
0414          end
0415    *$ Name:
0416    *$  nvxeas - Compute image coordinates from earth coordinates 
0417    *$
0418    *$ Interface:
0419    *$  integer function
0420    *$  nvxeas( real e1, real e2, real e3,
0421    *$          real lin, real ele, real dummy)
0422    *$
0423    *$ Input:
0424    *$  e1  - latitude  or x
0425    *$  e2  - longitude or y
0426    *$  e3  - height    or z
0427    *$
0428    *$ Input and Output:
0429    *$  none
0430    *$
0431    *$ Output:
0432    *$  lin - image line
0433    *$  ele - image element
0434    *$  dummy   - (unused)
0435    *$
0436    *$ Return values:
0437    *$   0  - success
0438    *$  -1  - input data physically valid but not navigable
0439    *$        given the specified projection
0440    *$  -2  - input data exceed physical limits
0441    *$  -6  - module not initialized
0442    *$
0443    *$ Remarks:
0444    *$  The navigation module must first be initialized with
0445    *$  a call to nvxini(). The input form (lat,lon) or (x,y,z)
0446    *$  depends on the last call to nvxini() with option 2.
0447    *$  Input longitude may be in the range -360 to +360;
0448    *$  values outside this range will not be de-navigated.
0449    *$  Height (hgt) is ignored.
0450    *$
0451    *$ Categories: 
0452    *$  navigation  
0453    
0454          INTEGER FUNCTION NVXEAS( e1, e2, e3, lin, ele, dummy)
0455    
0456          implicit      NONE
0457    
0458    C     // Interface variables (formal arguments)
0459    
0460          real          e1      ! Earth coordinate 1
0461          real          e2      ! Earth coordinate 2
0462          real          e3      ! Earth coordinate 3
0463          real          lin     ! image line to navigate
0464          real          ele     ! image element to navigate
0465          real          dummy       ! (unused argument)
0466    
0467    C     // Local variables
0468    
0469          real          lat     ! latitude  (McIDAS convention)
0470          real          lon     ! longitude (McIDAS convention)
0471          real          hgt     ! height
0472          real          in_lon      ! input longitude (radians,
0473    C                   ! East positive)
0474          real          colat       ! colatitude
0475          real          radius      ! distance from pole on projection
0476          real          theta       ! angle from standard longitude on
0477    C                   ! projection surface
0478    
0479    C     ////////////////////////////////////////////////////////
0480    C     Common block variables and declaration.
0481    
0482    C     ALL CODE BETWEEN THE '//////' SEPARATORS MUST BE 
0483    C     DUPLICATED EXACTLY IN EACH NAVIGATION ROUTINE
0484    
0485    C     (A more maintenance-safe version would use ENTRY points
0486    C     rather than separate functions for the navigation API
0487    C     but entry points cannot be processed by 'convdlm.')
0488    
0489    C     // Common block contents: projection parameters
0490    
0491          real          Lin0        ! image line of pole
0492          real          Ele0        ! image element of pole
0493          real          Scale       ! km per unit image coordinate
0494    C                   ! (pixel)
0495          real          Lon0        ! standard longitude
0496          real          Colat0      ! standard colatitude
0497    
0498    C     // Common block contents: precomputed intermediate values
0499    
0500          real          Coscl       ! cosine(Colat0)
0501          real          Tancl       ! tangent(Colat0)
0502          real          Tancl2      ! tangent(Colat0/2)
0503          real          Mxtheta     ! limit of angle from std. lon
0504    C                   ! on projection surface
0505    
0506    C     // Common block contents: constants
0507    
0508          real          D2R     ! degrees to radians factor
0509          real          Pi
0510          real          Badreal     ! returned when navigation
0511    C                   ! cannot be done
0512          real          Erad        ! Earth radius
0513          logical           Init        ! initialized flag
0514          logical           Latlon      ! .TRUE. for lat/lon I/O
0515    
0516    
0517          COMMON/TANC/ Lin0, Ele0, Scale, Lon0, Colat0,
0518         &  Coscl, Tancl, Tancl2, Mxtheta,
0519         &  D2R, Pi, Badreal, Erad, Init, Latlon
0520          
0521    C     End of common block variables and declaration.
0522    C     /////////////////////////////////////////////////////
0523    
0524          lin       = Badreal
0525          ele       = Badreal
0526          dummy     = Badreal
0527    
0528    C     // verify that module is initialized
0529    
0530          if(.not.init) then
0531              NVXEAS = -6
0532              return
0533          end if
0534    
0535    C     // Preprocess input values. If mode is 'XYZ' first convert
0536    C     // from Cartesian to lat/lon. If mode is 'LL' just transcribe
0537    C     // from arguments.
0538    
0539          if(Latlon) then
0540          lat = e1
0541          lon = e2
0542          hgt = e3
0543          else
0544          call nxyzll( e1, e2, e3, lat, lon)
0545          hgt = 0.
0546          end if
0547    
0548    C     // check that input values are physically possible and
0549    C     // then convert to radians and East positive
0550    
0551          if ( lat.lt.-90. .or. lat.gt.90. ) then
0552             NVXEAS = -2
0553             return
0554          end if
0555          if( lon.le.-360..or.lon.gt.360.) then
0556              NVXEAS = -2
0557              return
0558          end if
0559    
0560          if( lat.eq.-90. .or. lat.eq.90. ) then
0561              NVXEAS = -1
0562              return
0563          end if
0564    
0565          colat     = Pi/2. - D2R*lat
0566          in_lon        = -D2R*lon
0567    
0568    C     // map longitude into range -Pi to Pi
0569    
0570          if(in_lon.le.-Pi) in_lon = in_lon + 2.*Pi
0571          if(in_lon.gt. Pi) in_lon = in_lon - 2.*Pi
0572    
0573    
0574    C     // Now trap South Pole. Though a physically possible latitude,
0575    C     // tan(colat/2) -> infinity there so it is not navigable
0576    
0577          if ( colat.eq.Pi ) then
0578            NVXEAS      = -1
0579            return
0580          end if
0581    
0582    
0583    C     // Compute radius and theta of point on projection surface.
0584    C     // Theta is tricky; you have to compute offset relative
0585    C     // to standard longitude, force that into -pi to +pi range,
0586    C     // and THEN scale by cos(Colat0)
0587    
0588          radius        = Erad * Tancl *( tan(colat/2.)/Tancl2 ) ** Coscl
0589          theta     = in_lon-Lon0
0590    
0591          if(theta.le.-Pi) theta = theta + 2.*Pi
0592          if(theta.gt. Pi) theta = theta - 2.*Pi
0593    
0594          theta     = Coscl * theta
0595    
0596    
0597    C     // Compute line and element
0598    
0599          lin       = Lin0 + radius*cos(theta)/Scale
0600          ele       = Ele0 + radius*sin(theta)/Scale
0601          dummy     = 0.
0602    
0603          NVXEAS        = 0
0604    
0605          return
0606          end
0607    
0608    
0609    *$ Name:
0610    *$  nvxopt - Perform supplemental navigation operations
0611    *$
0612    *$ Interface:
0613    *$  integer function
0614    *$  nvxopt(integer option, real xin(*),
0615    *$    real xout(*) )
0616    *$ Input:
0617    *$  option  - 'SCAL' compute projection scale
0618    *$  xin(1)  - latitude
0619    *$
0620    *$ Input and Output:
0621    *$  none
0622    *$
0623    *$ Output:
0624    *$  xout(1) - km per pixel at given latitude
0625    *$ Return values:
0626    *$   0  - success
0627    *$  -1  - input latitude physically valid, but projection
0628    *$        undefined or scale infinite there
0629    *$  -2  - input latitude exceeds physical limits
0630    *$  -5  - unrecognized option
0631    *$  -6  - module not initialized
0632    *$
0633    *$ Remarks:
0634    *$  The navigation module must first be initialized by
0635    *$  a call to nvxini().  Latitude is in degrees, north positive,
0636    *$  and must lie between -90. and +90.
0637    *$
0638    *$ Categories: 
0639    *$  navigation  
0640    
0641          INTEGER FUNCTION NVXOPT( option, xin, xout)
0642    
0643          implicit      NONE
0644    
0645    C     // Interface variables (formal arguments)
0646    
0647          integer           option      ! special service name (character
0648    C                   ! stored as integer)
0649          real          xin(*)      ! input vector
0650          real          xout(*)     ! output vector
0651    
0652    C     // Local variables
0653    
0654          character*4           copt        ! special service (character form)
0655          real          colat       ! input colatitude
0656    
0657    C     ////////////////////////////////////////////////////
0658    C     Common block variables and declaration.
0659    
0660    C     ALL CODE BETWEEN THE '//////' SEPARATORS MUST BE 
0661    C     DUPLICATED EXACTLY IN EACH NAVIGATION ROUTINE
0662    
0663    C     (A more maintenance-safe version would use ENTRY points
0664    C     rather than separate functions for the navigation API
0665    C     but entry points cannot be processed by 'convdlm.')
0666    
0667    C     // Common block contents: projection parameters
0668    
0669          real          Lin0        ! image line of pole
0670          real          Ele0        ! image element of pole
0671          real          Scale       ! km per unit image coordinate
0672    C                   ! (pixel)
0673          real          Lon0        ! standard longitude
0674          real          Colat0      ! standard colatitude
0675    
0676    C     // Common block contents: precomputed intermediate values
0677    
0678          real          Coscl       ! cosine(Colat0)
0679          real          Tancl       ! tangent(Colat0)
0680          real          Tancl2      ! tangent(Colat0/2)
0681          real          Mxtheta     ! limit of angle from std. lon
0682    C                   ! on projection surface
0683    
0684    C     // Common block contents: constants
0685    
0686          real          D2R     ! degrees to radians factor
0687          real          Pi
0688          real          Badreal     ! returned when navigation
0689    C                   ! cannot be done
0690          real          Erad        ! Earth radius
0691          logical           Init        ! initialized flag
0692          logical           Latlon      ! .TRUE. for lat/lon I/O
0693    
0694    
0695          COMMON/TANC/ Lin0, Ele0, Scale, Lon0, Colat0,
0696         &  Coscl, Tancl, Tancl2, Mxtheta,
0697         &  D2R, Pi, Badreal, Erad, Init, Latlon
0698          
0699    C     End of common block variables and declaration.
0700    C     ////////////////////////////////////////////////////
0701    
0702    
0703          xout(1)       = Badreal
0704    
0705    C     // verify initialized module
0706    
0707          if(.not.init) then
0708              NVXOPT        = -6
0709              return
0710          end if
0711    
0712    C     // Extract and interpret the option
0713    
0714          call movwc(option,copt)
0715    
0716          if(copt.eq.'SCAL') then
0717    
0718    C        //  Compute colatitude and make sure it is
0719    C        //  physically possible and navigable
0720    
0721         if ( xin(1).gt.90. .or. xin(1).lt.-90. ) then
0722             NVXOPT     = -2
0723             return
0724         else if ( xin(1).eq.90. .or. xin(1).eq.-90. ) then
0725             NVXOPT     = -1
0726             return
0727         end if
0728    
0729         colat  = Pi/2. - D2R*xin(1)
0730    
0731    C        // Now compute actual scale for this colatitude
0732    
0733         xout(1)    = scale
0734         *           *(sin(Colat0)*(tan(colat/2.)/Tancl2)**Coscl)/sin(colat)
0735    
0736    C     else if(copt.eq.'????') 
0737    C        // Add code for additional options here
0738    
0739          else
0740         NVXOPT     = -5
0741         return
0742          end if
0743    
0744          NVXOPT = 0
0745
0746          return
0747          end

[Search Manual] [Table of Contents] [Go to Previous] [Go to Next]