Two subroutines are provided for Earth location/navigation for NESDIS Operational Products that are mapped onto a polar stereographic map base. These subroutines provide the capability of converting the line, sample location of a data point on the polar stereographic base to latitude and longitude (IJTOLL) or the reverse (LLTOIJ). With proper specification of the input parameters, these subroutines may be used in conjunction with the following:
Subroutines IJTOLL and LLTOIJ are very similar in structure to each other. They are described together where possible to avoid repetition. IJTOLL converts the i and j coordinates of a point to latitude and longitude, whereas LLTOIJ converts the latitude and longitude of a point to the i and j coordinates. Both require the origin of the map to be in the upper left corner. Several parameters must be input in order to tailor IJTOLL and LLTOIJ to a specific user's needs. Both subroutines have identical variables in their argument lists. These variables are defined in Table A-1.
|PRMLON||The prime longitude of the Polar Stereographic map base. This is defined as the meridian which is perpendicular to the base of the map and which goes from the pole to the bottom of the map. The map may be oriented as desired. The prime longitude can be expressed as either positive 0 to 360W or positive 0 to 180W and negative 0 to 180E.|
|SCALE||The number of desired grid points between the pole and the equator.|
|CENTI||The i coordinate representing the center of the array.|
|CENTJ||The j coordinate representing the center of the array|
|RI||The i coordinate of the point to be converted.|
|RJ||The j coordinate of the point to be converted.|
|PNTLAT||The latitude of the point written as positive 0 to 90N or negative 0 to 90S.|
|PNTLON||The longitude of the point written as either positive 0 to 360W or positive 0 to 180W and negative 0 to 180E.|
|Note: IJTOLL has RI and RJ for inputs with the corresponding PNTLAT and PNTLON as outputs, and LLTOIJ uses PNTLAT and PNTLON as inputs to calculate the resultant RI and RJ as outputs.|
Subroutine IJTOLL computes the angle between the horizontal and a line drawn from the pole to the point. From this angle and the prime longitude, the longitude of the point is obtained. IJTOLL also computes the latitude angle using various ratios and some simple trigonometry.
The equations governing the relation between i and j coordinates and latitude and longitude follow.
The distance from the point to the pole (in units of grid points) was derived from the Pythagorean theorem and is defined as:
The angle (in degrees) governing the latitude can easily be derived using the diagram of the Earth's cross-section in Figure A-1, and is defined as:
The angle (in degrees) governing the longitude is also easily derived from the diagram of a polar stereographic map and the i, j coordinates array in Figure A-2, and is defined as:
Thus, from Figure A-1, the latitude is:
From Figure A-2, the longitude is defined as:
Subroutines IJTOLL and LLTOIJ can be used for an array containing one or two hemispheres. For the case of two hemispheres, IJTOLL and LLTOIJ assume that they are arranged as shown in Figure A-3. For purposes of illustration only, the array shown is dimensioned as two 512 x 512s or one 512 x 1024.
These equations make up the algorithm for converting from i, j coordinates to latitude and longitude of a point. There is one special case where these equations blow up, which happens when the point is at a pole. This makes POLDIS=0.0 and is handled by setting the latitude to ±90 and the longitude to the prime longitude. There are a series of tests in IJTOLL which compensate for the longitude according to which quadrant the point resides. These tests produce the correct longitude for both the Northern and Southern Hemispheres.
If an i, j value is input which falls outside the Earth's disk but inside the array itself, the algorithm checks whether POLDIS exceeds SCALE. If it does not, then it sets the latitude and longitude to -999.9.
The algorithms which LLTOIJ are based upon can easily be derived by working the IJTOLL equations backwards. LLTOIJ computes the distance of the projected point from the pole using the latitude angle and some simple trigonometry. In addition, LLTOIJ computes the angular difference between the prime longitude and the longitude of the point. These two parameters are then used to compute the i and j coordinates of the point.
In Figure A-3, the upper 512 x 512 array contains the Northern Hemisphere while the lower 512 x 512 array contains the Southern Hemisphere. Subroutine IJTOLL will automatically return the latitude and longitude (or LLTOIJ will return the i and j coordinates) for the Southern Hemisphere although the CENTI and CENTJ input parameters specifically refer to the Northern Hemisphere.
The arrangement shown in Figure A-3 is not mandatory in order for IJTOLL and LLTOIJ to function properly. For example, the data for the Northern Hemisphere could be contained in one 512 x 512 array and the data for the Southern Hemisphere could be contained in an entirely separate 512 x 512 array. For that case, IJTOLL and LLTOIJ could be used for the Northern Hemisphere with the proper input arguments. IJTOLL and LLTOIJ could also be used for the Southern Hemisphere using the same input arguments except the prime longitude which would have to be changed. (It is assumed that the orientation of the hemispheres with the arrays is the same as in Figure A-3.)
IJTOLL and LLTOIJ represent the i and j coordinates in such a way that the i values represent the columns and the j values represent the rows of the array. Some of the principal points of the Earth's disk have been labeled with their proper i and j coordinates in parentheses in Figure A-3.
For the example in Figure A-3, the input arguments would be set as follows:
PRMLON = +80.0
SCALE = 256.0
CENTI = 256.0
CENTJ = 256.0
The input parameters can be adjusted so that IJTOLL and LLTOIJ can be applied to any case. For example, the array for the Northern Hemisphere is a 512 x 512 array. But unlike the case shown in Figure A-3, the Earth's disk is smaller than the boundary of the array. There are 8 blank-filled rows at the top and 4 at the bottom. Plus, there are 4 blank-filled columns on the left and 8 on the right of the Earth's disk. In this case, SCALE = 256.0.
CENTI = 250.0 + 4.0 = 254.0 and the CENTJ = 250.0 + 8.0 = 258.0. See Figure A-4 for an illustration of this example.
In LLTOIJ, there is a problem with obtaining the proper i and j coordinates for any point along the equator in the Southern Hemisphere. This is caused by the fact that the computer represents -0.0 as 0.0; therefore, making it impossible to represent the Southern Hemisphere equator. To obtain the approximate i and j coordinates for the equator in the Southern Hemisphere, input a number between 0.0 and -1.0 for the latitude. Or alternately input the latitude as 0.0, and add the value of 2 x CENTJ to the RJ value obtained. (Adding 2 x CENTJ to RJ would give the RJ the correct value for the Southern Hemisphere.)
A listing of the FORTRAN code for subroutines IJTOLL and LLTOIJ follows.
SUBROUTINE IJTOLL (PRMLON, SCALE, CENTI, CENTJ, RI, RJ, PNTLAT, PNTLON)
C THIS SUBROUTINE CONVERTS I AND J COORDINATES OF A POINT ON C A POLAR STEREOGRAPHIC MAP TO LATITUDE AND LONGITUDE. C C INPUT C PRMLON PRIME LONGITUDE OF MAP BASE IN DEGREES. CAN BE C EITHER POSITIVE 0-360W OR POSITIVE 0-180W AND C NEGATIVE 0-180E. C SCALE DISTANCE FROM POLE TO EQUATOR IN GRID POINTS. C CENTI, CENTER OF THE ARRAY IN I AND J COORDINATES. C CENTJ C RI, RJ THE I AND J COORDINATE OF THE POINT TO BE CONVERTED. C OUTPUT C PNTLAT THE LATITUDE OF THE POINT IN DEGREES (POSITIVE C 0-90N AND NEGATIVE 0-90S). C PNTLON THE LONGITUDE OF THE POINT IN DEGREES WITH SAME C SIGN CONVENTION AS PRMLON. SOUTH=2.0*CENTJ RY=RJ RADDEG=0.01745329 C CHECK WHETHER POINT IS IN SOUTHERN HEMISPHERE AND ADJUST C ACCORDINGLY. IF (RJ.GT.SOUTH) RJ=2.0*SOUTH-RJ C COMPUTE THE DISTANCE BETWEEN THE POLE AND THE POINT. IF C POLDIS IS GREATER THAN SCALE, SET PNTLAT AND PNTLON TO C -999.9 TO INDICATE AREA OUTSIDE OF EARTH'S DISK. IF POLDIS C EQUALS SCALE SET PNTLAT AND PNTLON EQUAL TO A POLE. POLDIS=SQRT(ABS(CENTI-RI)**2.0+ABS(CENTJ-RJ)**2.0) IF (POLDIS.GT.SCALE) GO TO 50 IF (POLDIS.EQ.0) GO TO 30 C COMPUTE PNTLAT AND PNTLON FROM ANG1 AND ANG2, THEN ADJUST C PNTLAT ACCORDING TO RELATION WITH CENTER OF ARRAY. ANG1=ATAN(POLDIS/SCALE) ANG2=ACOS((CENTJ-RJ)/POLDIS) IF (RJ.GT.CENTJ) ANG2=ACOS((RJ-CENTJ)/POLDIS) PNTLON=PRMLON-ANG2/RADDEG IF (RI.GT.CENTI.AND.RJ.LE.CENTJ) PNTLON=PRMLON+ANG2/RADDEG IF (RI.LT.CENTI.AND.RJ.GE.CENTJ) PNTLON=PRMLON+ANG2/RADDEG PNTLAT=90.0-2.0*ANG1/RADDEG IF (RI.LE.CENTI.AND.RJ.LT.CENTJ) PNTLON=PNTLON+180.0 IF (RI.GT.CENTI.AND.RJ.LE.CENTJ) PNTLON=PNTLON+180.0 IF (RY.GT.SOUTH) PNTLAT=-PNTLAT RJ=RY RETURN 30 PNTLAT=90.0 PNTLON=PRMLON IF (RY.GT.SOUTH) PNTLAT=-PNTLAT RJ=RY RETURN 50 PNTLAT=-999.9 PNTLON=-999.9 RJ=RY RETURN END
SUBROUTINE LLTOIJ (PRMLON, SCALE, PNTLAT, PNTLON, CENTI, CENTJ, RI, RJ)
C THIS SUBROUTINE CONVERTS LATITUDE AND LONGITUDE OF A POINT C ON A POLAR STEREOGRAPHIC MAP TO I AND J COORDINATES. C C INPUT C PRMLON PRIME LONGITUDE OF MAP BASE IN DEGREES. CAN BE C EITHER POSITIVE 0-360W OR POSITIVE 0-180W AND C NEGATIVE 0-180E. C SCALE DISTANCE FROM POLE TO EQUATOR IN GRID POINTS. C CENTI, CENTER OF THE ARRAY IN I AND J COORDINATES. C CENTJ C PNTLAT THE LATITUDE OF THE POINT TO BE CONVERTED C IN DEGREES (POSITIVE 0-90N AND NEGATIVE 0-90S). C PNTLON THE LONGITUDE OF THE POINT TO BE CONVERTED C IN DEGREES WITH SAME SIGN CONVENTION AS PRMLON. C C OUTPUT C RI THE I COORDINATE OF THE POINT. C RJ THE J COORDINATE OF THE POINT. C SOUTH=2.0*CENTJ RADDEG=0.01745329 C COMPUTE THE DISTANCE OF THE PROJECTED POINT FROM THE C POLE. ANG1=(90.0-ABS(PNTLAT))*RADDEG/2.0 POLDIS=SCALE*TAN(ANG1) C COMPUTE THE DIFFERENCE IN RADIANS BETWEEN THE PRIME C LONGITUDE AND THE LONGITUDE OF THE POINT. ANG2=(PRMLON-PNTLON)*RADDEG C COMPUTE THE I AND J COORDINATES OF THE POINT (RI AND RJ). RI=CENTI+POLDIS*SIN(ANG2) IF(PRMLON.LT.0.0) RI=CENTI-POLDIS*SIN(ANG2) IF(PNTLAT.LT.0.0) GO TO 10 RJ=CENTJ+POLDIS*COS(ANG2) RETURN 10 RJ=CENTJ-POLDIS*COS(ANG2) IF(PNTLAT.LT.0.0) RJ=RJ+SOUTH RETURN END Source: http://www.ncdc.noaa.gov