Paul Michael writes:
> I've a "backwards" question. Is it possible to use the EZMAP and/or
> EZMAPA routines to obtain Area Identifiers associated with a Latitude-
> Longitude pair? (Thew Area Identifiers being the ones liste on pp 271-286
> of the Aug. 1987 Users Guide).
>
> Actualy I could do with less - could I tell whether a lat-long point is
> over water or land?
The program attached below illustrates the use of a function called IWHERE
having arguments RLAT, RLON, and IAMA. For a given latitude RLAT and
longitude RLON, the value of IWHERE(RLAT,RLON,IAMA) is 1 if that point is
over water, greater than 1 if that point is over land, 0 or less if neither
(which should not be possible). IAMA is an area map array generated by a
call to the subroutine IWINIT. This program could easily be beefed up to
return the area identifiers of which you speak.
Program follows:
PROGRAM TESTIT
C
C Define the size of the area map to be used.
C
PARAMETER (LAMA=125000)
C
C Declare the area map and a couple of arrays in which to put the
C coordinates of a box to be filled.
C
DIMENSION IAMA(LAMA),XCBX(4),YCBX(4)
C
C Open GKS.
C
CALL OPNGKS
C
C Define a color to use for ocean and another color to use for land.
C
CALL GSCR (1,2,0.,1.,1.)
CALL GSCR (1,3,1.,0.,1.)
C
C Call IWINIT to generate the area map IWHERE will use.
C
CALL IWINIT (IAMA,LAMA)
C
C Call SET to define user coordinates to be longitude and latitude.
C
CALL SET (.02,.98,.26,.74,-180.,180.,-90.,90.,1)
C
C For each box in a 360x180 grid of boxes, use IWHERE to find out if
C the center point of the box is over neither ocean nor land (ITMP < 1),
C over ocean (ITMP = 1) or over some land area (ITMP > 1), and color
C the box appropriately.
C
DO 102 J=-90,89
RLAT=.5+REAL(J)
DO 101 I=-180,179
RLON=.5+REAL(I)
ITMP=IWHERE(RLAT,RLON,IAMA)
IF (ITMP.NE.0) THEN
IF (ITMP.EQ.1) THEN
CALL GSFACI (2)
ELSE
CALL GSFACI (3)
END IF
XCBX(1)=RLON-.5
YCBX(1)=RLAT-.5
XCBX(2)=RLON+.5
YCBX(2)=RLAT-.5
XCBX(3)=RLON+.5
YCBX(3)=RLAT+.5
XCBX(4)=RLON-.5
YCBX(4)=RLAT+.5
CALL GFA (4,XCBX,YCBX)
END IF
101 CONTINUE
102 CONTINUE
C
C Advance the frame.
C
CALL FRAME
C
C Close GKS.
C
CALL CLSGKS
C
C Done.
C
STOP
C
END
SUBROUTINE IWINIT (IAMA,LAMA)
C
C Generate an area map for the function IWHERE to use.
C
DIMENSION IAMA(LAMA)
C
C Initialize EZMAP. Note that I'm assuming EZMAP is in its default
C state at this point. By default, we get a cylindrical equidistant
C projection and nothing but continental outlines. Vertical stripping
C is turned off so that only group 1 edges go into the area map.
C
CALL MPSETI ('VS - VERTICAL STRIPPING FLAG',0)
CALL MAPINT
C
C Override the SET call done by EZMAP, forcing the map to fill the
C entire plotter frame.
C
CALL GETSET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
CALL SET ( 0., 1., 0., 1.,XWDL,XWDR,YWDB,YWDT,LNLG)
C
C Initialize the area-map array.
C
CALL ARINAM (IAMA,LAMA)
C
C Uncomment the following lines to see what's going into the area map.
C
C CALL MAPLOT
C CALL FRAME
C
C Put continental outlines into the area map.
C
CALL MAPBLA (IAMA)
C
C Pre-process the area map.
C
CALL ARPRAM (IAMA,0,0,0)
C
C Uncomment the following statement to print the amount of space that
C was actually used in the area map.
C
C PRINT * , 'SPACE USED IN AREA MAP IS ',
C + IAMA(1)-(IAMA(6)-IAMA(5)-1)
C
C Done.
C
RETURN
C
END
FUNCTION IWHERE (RLAT,RLON,IAMA)
C
C For a given latitude RLAT and longitude RLON, the value of IWHERE is
C 1 if that point is over water, greater than 1 if that point is over
C land, 0 or less if neither (which should not be possible). IAMA is
C an area map array generated by a call to the subroutine IWINIT.
C
DIMENSION IAMA(*)
C
C Translate the given latitude and longitude into X and Y positions
C commensurate with what's in the area map and then translate those
C into the current user system, whatever it may be. This ensures
C that ARGTAI will look at the proper point in the area map, no matter
C what SET call has been done.
C
XPOS=CFUX(MOD((RLON+540.)/360.,1.))
YPOS=CFUY((RLAT+90.)/180.)
C
C Call the areas routine ARGTAI to get area-identifier information
C about the point (XPOS,YPOS). We should get back just one number
C in the "array" IAAI.
C
CALL ARGTAI (IAMA,XPOS,YPOS,IAAI,IAGI,1,NOIR,1)
C
C Interpret the number. IAAI < 1 implies "not over map". For other
C values, return the recommended color index suggested by the routine
C MAPACI, which will always be 1 for ocean areas and greater than 1
C for land areas.
C
IF (IAAI.LE.0) THEN
IWHERE=0
ELSE
IWHERE=MAPACI(IAAI)
END IF
C
C Done.
C
RETURN
C
END
This archive was generated by hypermail 2b29 : Wed Jun 28 2000 - 09:40:26 MDT