On March 23, Peizhong asked two questions, the second of which was as
follows (I don't know the answer to the first question):
> Another question: When I am trying to do map contouring, I have to
> convert the data from lat/lon into X/Y at first. Is there an easy way
> that I can plot map contouring directory from the data in lat/lon
> coordinator?
I think the attached program does what he wants to do. It may require a
bit of studying. (The bottom line is that it's pretty easy to do, once
you know how.)
In the example program, the array ZDAT is dimensioned 71 by 51. The element
ZDAT(I,J) is associated with the longitude -171+I (in degrees) and latitude
4+J (in degrees), so the data covers the area bounded by the meridians at
170W and 100W and the parallels at 5N and 55N. For purposes of illustration,
ZDAT is filled with dummy data.
First, we do the appropriate EZMAP calls to draw a map of the globe, in
orthographic projection, centered on the region covered by the data array.
Then, we
CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
to tell CONPACK not to call SET for itself, but to depend on the SET call
done by EZMAP. (The arguments of SET specify how "user coordinates" are to
be mapped to the plotter frame.) Then, we
CALL CPSETI ('MAP - MAPPING FUNCTION',1)
to tell CONPACK to transform the coordinates associated with the data array
using mapping function 1 (which, by default, is a mapping from lon/lat to
the user coordinate system defined by the SET call done by EZMAP). We also
CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12)
to tell CONPACK what value will be returned by the mapping routine when
a given lat/lon is out of sight around the edge of the globe (1.E12 by
convention), so that it can react properly.
Then we use the four calls
CALL CPSETR ('XC1 - LONGITUDE AT I = 1',-170.)
CALL CPSETR ('XCM - LONGITUDE AT I = M',-100.)
CALL CPSETR ('YC1 - LATITUDE AT J = 1' , 5.)
CALL CPSETR ('YCN - LATITUDE AT J = N' , 55.)
to tell CONPACK what values are to be associated with the extreme values
of the indices of the data array. (Linear interpolation determines what
values are to be associated with intervening values of the indices.)
Then, we do the calls necessary to initialize the drawing of the contour
plot and to draw the contour lines. The effect is that CONPACK generates
contour lines defined by X/Y coordinates where X is between -170 and -100
and Y is between 5 and 55, transforms those lines using a mapping routine
that interprets X values as longitudes and Y values as latitudes to get
coordinates in the current user system and then draws those contour lines
using a user-system-to-plotter-frame mapping that positions them correctly
relative to the map already drawn.
Run the program and take a look at the metafile it produces.
Program follows:
PROGRAM TESTIT
C
C This program demonstrates simple contouring on a globe.
C
C Declare an array to hold the data to be contoured.
C
DIMENSION ZDAT(71,51)
C
C Declare the required real and integer workspaces for CONPACK.
C
DIMENSION RWRK(5000),IWRK(1000)
C
C Open GKS.
C
CALL OPNGKS
C
C Generate an array of test data.
C
CALL GENDAT (ZDAT,71,71,51,21,21,-1.,1.)
C
C Initialize EZMAP and draw a map background.
C
CALL MAPPOS (.05,.95,.01,.91)
CALL MAPROJ ('OR - ORTHOGRAPHIC',40.,-135.,0.)
CALL MAPSET ('MA - MAXIMAL AREA',0.,0.,0.,0.)
C
CALL MAPDRW
C
C Tell CONPACK to do no SET call (EZMAP has done it).
C
CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
C
C Tell CONPACK to use more contour levels.
C
CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',32)
C
C Turn on the drawing of the grid edge ("contour line number -1") and
C thicken it somewhat.
C
CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1)
CALL CPSETI ('CLU - CONTOUR LEVEL USE FLAG',1)
CALL CPSETR ('CLL - CONTOUR LEVEL LINE WIDTH',2.)
C
C Tell CONPACK to use EZMAP for mapping and what the out-of-range
C signal is.
C
CALL CPSETI ('MAP - MAPPING FUNCTION',1)
CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12)
C
C Tell CONPACK what range of coordinates to use.
C
CALL CPSETR ('XC1 - LONGITUDE AT I = 1',-170.)
CALL CPSETR ('XCM - LONGITUDE AT I = M',-100.)
CALL CPSETR ('YC1 - LATITUDE AT J = 1' , 5.)
CALL CPSETR ('YCN - LATITUDE AT J = N' , 55.)
C
C Initialize the drawing of the contour plot.
C
CALL CPRECT (ZDAT,71,71,51,RWRK,5000,IWRK,1000)
C
C Force the selection of contour levels so that associated quantities
C may be tweaked.
C
CALL CPPKCL (ZDAT,RWRK,IWRK)
C
C Increase the line width for labelled levels.
C
CALL CPGETI ('NCL - NUMBER OF CONTOUR LEVELS',NCLV)
C
DO 101 ICLV=1,NCLV
CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',ICLV)
CALL CPGETI ('CLU - CONTOUR LEVEL USE FLAG',ICLU)
IF (ICLU.EQ.3) THEN
CALL CPSETI ('CLL - CONTOUR-LINE LINE WIDTH',2)
END IF
101 CONTINUE
C
C Draw the contour lines.
C
CALL CPCLDR (ZDAT,RWRK,IWRK)
C
C Advance the frame.
C
CALL FRAME
C
C Close GKS.
C
CALL CLSGKS
C
C Done.
C
STOP
C
END
SUBROUTINE GENDAT (DATA,IDIM,M,N,MLOW,MHGH,DLOW,DHGH)
C
C This is a routine to generate test data for two-dimensional graphics
C routines. Given an array "DATA", dimensioned "IDIM x 1", it fills
C the sub-array ((DATA(I,J),I=1,M),J=1,N) with a two-dimensional field
C of data having approximately "MLOW" lows and "MHGH" highs, a minimum
C value of exactly "DLOW" and a maximum value of exactly "DHGH".
C
C "MLOW" and "MHGH" are each forced to be greater than or equal to 1
C and less than or equal to 25.
C
C The function used is a sum of exponentials.
C
DIMENSION DATA(IDIM,1),CCNT(3,50)
C
FOVM=9./FLOAT(M)
FOVN=9./FLOAT(N)
C
NLOW=MAX0(1,MIN0(25,MLOW))
NHGH=MAX0(1,MIN0(25,MHGH))
NCNT=NLOW+NHGH
C
DO 101 K=1,NCNT
CCNT(1,K)=1.+(FLOAT(M)-1.)*FRAN()
CCNT(2,K)=1.+(FLOAT(N)-1.)*FRAN()
IF (K.LE.NLOW) THEN
CCNT(3,K)=-1.
ELSE
CCNT(3,K)=+1.
END IF
101 CONTINUE
C
DMIN=+1.E36
DMAX=-1.E36
DO 104 J=1,N
DO 103 I=1,M
DATA(I,J)=.5*(DLOW+DHGH)
DO 102 K=1,NCNT
TEMP=-((FOVM*(FLOAT(I)-CCNT(1,K)))**2+
+ (FOVN*(FLOAT(J)-CCNT(2,K)))**2)
IF (TEMP.GE.-20.) DATA(I,J)=DATA(I,J)+
+ .5*(DHGH-DLOW)*CCNT(3,K)*EXP(TEMP)
102 CONTINUE
DMIN=AMIN1(DMIN,DATA(I,J))
DMAX=AMAX1(DMAX,DATA(I,J))
103 CONTINUE
104 CONTINUE
C
DO 106 J=1,N
DO 105 I=1,M
DATA(I,J)=(DATA(I,J)-DMIN)/(DMAX-DMIN)*(DHGH-DLOW)+DLOW
105 CONTINUE
106 CONTINUE
C
RETURN
C
END
FUNCTION FRAN ()
C
C This routine generates "random" numbers for GENDAT.
C
DIMENSION RSEQ (100)
C
SAVE ISEQ
C
DATA RSEQ / .749,.973,.666,.804,.081,.483,.919,.903,.951,.960 ,
+ .039,.269,.270,.756,.222,.478,.621,.063,.550,.798 ,
+ .027,.569,.149,.697,.451,.738,.508,.041,.266,.249 ,
+ .019,.191,.266,.625,.492,.940,.508,.406,.972,.311 ,
+ .757,.378,.299,.536,.619,.844,.342,.295,.447,.499 ,
+ .688,.193,.225,.520,.954,.749,.997,.693,.217,.273 ,
+ .961,.948,.902,.104,.495,.257,.524,.100,.492,.347 ,
+ .981,.019,.225,.806,.678,.710,.235,.600,.994,.758 ,
+ .682,.373,.009,.469,.203,.730,.588,.603,.213,.495 ,
+ .884,.032,.185,.127,.010,.180,.689,.354,.372,.429 /
C
DATA ISEQ / 0 /
C
ISEQ=MOD(ISEQ,100)+1
FRAN=RSEQ(ISEQ)
C
RETURN
C
END
This archive was generated by hypermail 2b29 : Wed Jun 28 2000 - 09:44:55 MDT