McIDAS Programmer's Manual
Version 2015
[Search Manual] [Table of Contents] [Go to Previous] [Go to Next]
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.
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. |
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.
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.
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 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 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.
This section contains some navigation-specific terms that may be unfamiliar to you. They are defined below.
The McIDAS-X library provides Application Program Interfaces for geographic calculations as well as grid, image and frame navigation. These APIs are described below.
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 |
|---|---|
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 |
|---|---|
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.
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 |
|---|---|
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) |
Use the appropriate image navigation API functions from the table below if your application must perform one of these tasks:
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.
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.
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.
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 |
|---|---|
nvxdmsp, nvxgoes, nvxgraf, nvxgvar, nvxlamb, nvxmerc, nvxmsat, nvxps, nvxradr, nvxrect, nvxsin, nvxtiro |
|
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.
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.
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.
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.
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.
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.
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.
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:
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.
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.
All instances of 'TANC' navigation must be specified with these parameters.
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:
or a McIDAS-X area format image file:
call araput(aranum,4*DIRSIZ,4*NVWDS,nvblk) |
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.
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.
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.
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]