The question was as follows:
> how can I draw an ARROW from point A (alon, alat) to point B (blon,blat)
> on a map, any suggestion is welcome, thanks.
>
> Yanjun
Appended below is a test program that may do what you want. I have made
some assumptions that may or may not be correct for your purposes, but I
think what I have done is the most difficult thing you might want to do:
it draws the arrow along the shortest great circle route from point A to
point B, using short segments (the size of the segments being determined
by a PARAMETER statement in subroutine DRSGCR). The width and length of
the arrowhead are set in a DATA statement in the main program.
This was constructed from bits and pieces of code that I had lying about
and may not be the simplest solution, but I think it does work. Tell me
if it doesn't do what you want and I will try to modify it accordingly.
Thanks,
Dave Kennison
Program follows:
----------------
PROGRAM TESTIT
C
C This program draws an arrow from the point A to the point B on the
C surface of the globe.
C
C NOTE: See the routine DRSGCR for parameters controlling the size of
C the little segments used to draw the arrow.
C
C Define some points A and B to try out.
C
DIMENSION ALAT(4),ALON(4),BLAT(4),BLON(4)
C
DATA ALAT(1),ALON(1),BLAT(1),BLON(1) / -15., -20., 25., 25. /
DATA ALAT(2),ALON(2),BLAT(2),BLON(2) / -25., -10., 15., -20. /
DATA ALAT(3),ALON(3),BLAT(3),BLON(3) / -10., -10., -15., 25. /
DATA ALAT(4),ALON(4),BLAT(4),BLON(4) / -15., -15., -25., 25. /
C
C Define the length and width of the arrowhead.
C
DATA ALTH,AWTH / 2. , 1. / ! arrowhead length and width
C
C Open GKS.
C
CALL OPNGKS
C
C Define a color index corresponding to the color red.
C
CALL GSCR (1,2,1.,0.,0.)
C
C Initialize EZMAP.
C
CALL MAPROJ ('OR',0.,0.,0.)
CALL MAPSET ('AN',40.,40.,40.,40.)
C
C Draw a simple map.
C
CALL MAPDRW
C
C Change the current polyline color to red and double the line width.
C
CALL PLOTIT (0.,0.,2)
CALL GSPLCI (2)
CALL GSLWSC (2.)
C
C Draw the arrows.
C
DO 101 I=1,4
CALL DARROW (ALAT(I),ALON(I),BLAT(I),BLON(I),ALTH,AWTH)
101 CONTINUE
C
C Set the polyline color and line width back to the defaults.
C
CALL PLOTIT (0.,0.,2)
CALL GSPLCI (1)
CALL GSLWSC (1.)
C
C Advance the frame.
C
CALL FRAME
C
C Close GKS.
C
CALL CLSGKS
C
C Done.
C
STOP
C
END
SUBROUTINE DARROW (ALAT,ALON,BLAT,BLON,ALTH,AWTH)
C
C Draw the body of the arrow.
C
CALL DRSGCR (ALAT,ALON,BLAT,BLON)
C
C Find the latitude and longitude of one end point of the arrowhead
C and draw it.
C
CALL GTLTLN (BLAT,BLON,ALAT,ALON,+AWTH/2.,+ALTH,CLAT,CLON)
CALL DRSGCR (BLAT,BLON,CLAT,CLON)
C
C Find the latitude and longitude of the other end point of the
C arrowhead and draw it.
C
CALL GTLTLN (BLAT,BLON,ALAT,ALON,-AWTH/2.,+ALTH,CLAT,CLON)
CALL DRSGCR (BLAT,BLON,CLAT,CLON)
C
C Done.
C
RETURN
C
END
SUBROUTINE GTLTLN (ALAT,ALON,BLAT,BLON,CLAT,CLON,DLAT,DLON)
C
C (GTLTLN = GeT LaTitude and LoNgitude)
C
C This routine rotates the great circle route between points A and B on
C the globe so as to put A at the intersection of the prime meridian and
C the equator and B to the east of A on the equator. It then takes the
C point at latitude CLAT and longitude CLON in that system and rotates
C it back to find the point with latitude DLAT and longitude DLON that
C lies in the same position relative to the original great arc AB as the
C point with latitude CLAT and longitude CLON does relative to the
C rotated arc AB.
C
C Constants for conversion from degrees to radians and vice-versa.
C
DATA DTOR / .017453292519943 /
DATA RTOD / 57.2957795130823 /
C
C Get X, Y, and Z coordinates for point A.
C
XCOA=COS(DTOR*ALAT)*COS(DTOR*ALON)
YCOA=COS(DTOR*ALAT)*SIN(DTOR*ALON)
ZCOA=SIN(DTOR*ALAT)
C
C Get X, Y, and Z coordinates for point B.
C
XCOB=COS(DTOR*BLAT)*COS(DTOR*BLON)
YCOB=COS(DTOR*BLAT)*SIN(DTOR*BLON)
ZCOB=SIN(DTOR*BLAT)
C
C Get X, Y, and Z coordinates for point C.
C
XCOC=COS(DTOR*CLAT)*COS(DTOR*CLON)
YCOC=COS(DTOR*CLAT)*SIN(DTOR*CLON)
ZCOC=SIN(DTOR*CLAT)
C
C Rotate points A and B about the Z axis so as to put point A at
C longitude zero.
C
CALL NGRITD (3,-ALON,XCOA,YCOA,ZCOA)
CALL NGRITD (3,-ALON,XCOB,YCOB,ZCOB)
C
C Rotate points A and B about the Y axis so as to put point A at
C latitude zero and longitude zero.
C
CALL NGRITD (2,+ALAT,XCOA,YCOA,ZCOA)
CALL NGRITD (2,+ALAT,XCOB,YCOB,ZCOB)
C
C Rotate points A and B about the X axis so as to put point B on
C the equator.
C
ANGL=RTOD*ATAN2(ZCOB,YCOB)
C
CALL NGRITD (1,-ANGL,XCOA,YCOA,ZCOA)
CALL NGRITD (1,-ANGL,XCOB,YCOB,ZCOB)
C
C Now, inverse rotate the point C.
C
CALL NGRITD (1,+ANGL,XCOC,YCOC,ZCOC)
CALL NGRITD (2,-ALAT,XCOC,YCOC,ZCOC)
CALL NGRITD (3,+ALON,XCOC,YCOC,ZCOC)
C
C Get the latitude and longitude of point D and return them.
C
DLAT=RTOD*ASIN(ZCOC)
DLON=RTOD*ATAN2(YCOC,XCOC)
C
RETURN
C
END
SUBROUTINE DRSGCR (ALAT,ALON,BLAT,BLON)
C
C (DRSGCR = DRaw Shortest Great Circle Route)
C
C This routine draws the shortest great circle route joining two points
C on the globe.
C
PARAMETER (MPTS=360,SIZE=.5) ! MPTS = INT(180./SIZE) + 1
C
DIMENSION QLAT(MPTS),QLON(MPTS)
C
NPTS=MAX(1,MIN(MPTS,INT(ADSGCR(ALAT,ALON,BLAT,BLON)/SIZE)))
C
CALL MAPGCI (ALAT,ALON,BLAT,BLON,NPTS,QLAT,QLON)
C
CALL MAPIT (ALAT,ALON,0,IDIS)
C
DO 101 I=1,NPTS
CALL MAPIT (QLAT(I),QLON(I),1)
101 CONTINUE
C
CALL MAPIT (BLAT,BLON,1)
C
CALL MAPIQ
C
RETURN
C
END
FUNCTION ADSGCR (ALAT,ALON,BLAT,BLON)
C
C (ADSGCR = Angle in Degrees along Shortest Great Circle Route)
C
C This function returns the shortest great circle distance, in degrees,
C between two points, A and B, on the surface of the globe.
C
DATA DTOR / .017453292519943 /
DATA RTOD / 57.2957795130823 /
C
XCOA=COS(DTOR*ALAT)*COS(DTOR*ALON)
YCOA=COS(DTOR*ALAT)*SIN(DTOR*ALON)
ZCOA=SIN(DTOR*ALAT)
C
XCOB=COS(DTOR*BLAT)*COS(DTOR*BLON)
YCOB=COS(DTOR*BLAT)*SIN(DTOR*BLON)
ZCOB=SIN(DTOR*BLAT)
C
ADSGCR=2.*RTOD*ASIN(SQRT((XCOA-XCOB)**2+
+ (YCOA-YCOB)**2+
+ (ZCOA-ZCOB)**2)/2.)
C
RETURN
C
END
SUBROUTINE NGRITD (IAXS,ANGL,UCRD,VCRD,WCRD)
C
C Define a multiplicative constant to convert from degrees to radians.
C
DATA DTOR / .017453292519943 /
C
C This routine rotates the point with coordinates (UCRD,VCRD,WCRD) by
C the angle ANGL about the axis specified by IAXS (1 for the U axis,
C 2 for the V axis, 3 for the W axis). A right-handed coordinate
C system is assumed.
C
SINA=SIN(DTOR*ANGL)
COSA=COS(DTOR*ANGL)
C
UTMP=UCRD
VTMP=VCRD
WTMP=WCRD
C
IF (IAXS.EQ.1) THEN
VCRD=VTMP*COSA-WTMP*SINA
WCRD=WTMP*COSA+VTMP*SINA
ELSE IF (IAXS.EQ.2) THEN
UCRD=UTMP*COSA+WTMP*SINA
WCRD=WTMP*COSA-UTMP*SINA
ELSE
UCRD=UTMP*COSA-VTMP*SINA
VCRD=VTMP*COSA+UTMP*SINA
END IF
C
RETURN
C
END
_______________________________________________
ncarg-talk mailing list
ncarg-talk@ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncarg-talk
This archive was generated by hypermail 2b29 : Thu Nov 14 2002 - 14:37:29 MST