Re: What place is this???

From: Dave Kennison (kennison@XXXXXX)
Date: Wed Nov 01 1995 - 14:24:40 MST


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:45:34 MDT