> I am trying to plot data on a staggered grid as I've tried to illustrate
> below:
>
> . . . . . . . . .
>
>
> . . . . . . . . .
>
>
> . . . . . . . . .
>
>
> In essence, the data are arranged in a diamond formation.
>
> Form the NG manual, I'd presume that the CPSPS2 routine (or it's updated)
> equivalent would be the one to use.
> Ah, but wait!!! For some unknown reason, we do not have ANY of the routines
> for interpolating data - we don't seem to have any of the following
> routines installed:
>
> 1. CPSPS1 for sparse gridded data,
> 2. CPSPS2 for irregularly spaced gridded data,
> 3. IDSFFT for non-gridded data
> ...CPMPXY I've never tried to write, so I don't know if there's a default
> version.
>
> Essentially, the version of NG installed on our machine has not been used
> to do any interpolation from irregular or non-gridded data because those
> routines don't seem to exist here. Is it possible to get these from you now?
>
> For info, we're running
> NCAR Graphics - UNIX Version 3.1.1
> Copyright (C) 1991 - All Rights Reserved
> University Corporation for Atmospheric Research
>
> with different versions of ctrans and idt supplied by you:
> ctrans - Version 3.2.1
> idt - Version 3.2.1
In version 3.1.1 of NCAR Graphics, you should find the routine CPSPRS (the
original of what became CPSPS1). CPSPS1 and CPSPS2 are not there. IDSFFT
*may* be there, but only as part of an example.
Getting the current versions of these routines is not recommended; they
would have to be pretty heavily modified to run with what you now have.
I think your best bet is to use a version of CPMPXY tailored to handle a
staggered grid. I have attached a program showing how to do this for a
specific case. I have tested this program with the latest version of NCAR
Graphics, but I can't think of anything that would prevent it from running
with version 3.1.1 (famous last words); try it out and let me know what
happens.
Dave Kennison
Program follows:
PROGRAM TESTIT
C
C ZDAT is an array of data to be contoured. Assume that, for odd values
C of J, ZDAT(I,J) is a function value at the point (REAL(I),REAL(J)) and
C that, for even values of J, ZDAT(I,J) is a function value at the point
C (REAL(I),REAL(J)+.5). Note that, in a real user case, the values at
C the ends of the odd-numbered rows may be missing, in which case some
C value should be interpolated in some reasonable way.
C
DIMENSION ZDAT(11,11)
C
C Declare work arrays required by CONPACK.
C
DIMENSION RWRK(2000),IWRK(1000)
C
C Define data for CONPACK such that we'll be able to look at the final
C plot and see if things worked all right. (Basically, the field values
C define a shallow bowl; the contours should all be roughly circular.)
C CONPACK will generate contours on this grid of data as if it were not
C staggered, but the contour lines will be mapped, just before plotting,
C by the version of CPMPXY provided below and this will map the contour
C lines to the staggered grid.
C
DO 102 J=1,11
YCRD=REAL(J)
XOFF=.5*REAL(MOD(J-1,2))
DO 101 I=1,11
XCRD=REAL(I)+XOFF
ZDAT(I,J)=((XCRD-6.)**2+(YCRD-6.)**2)/50.
101 CONTINUE
102 CONTINUE
C
C Open GKS.
C
CALL OPNGKS
C
C Turn GKS clipping off.
C
CALL GSCLIP (0)
C
C Define how the "user" coordinate system is to be mapped into the
C plotter frame.
C
CALL SET (.1,.9,.1,.9,1.,11.,1.,11.,1)
C
C Tell CONPACK not to do another SET call, but to depend on the one
C just done.
C
CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
C
C Tell CONPACK what X/Y values to associate with the extreme values
C of the data-array subscripts.
C
CALL CPSETR ('XC1 - X COORDINATE FOR I=1', 1.)
CALL CPSETR ('XCM - X COORDINATE FOR I=M',11.)
CALL CPSETR ('YC1 - Y COORDINATE FOR J=1', 1.)
CALL CPSETR ('YCN - Y COORDINATE FOR J=N',11.)
C
C Tell CONPACK how to map the X and Y coordinates that it uses to draw
C the contour lines. Telling it to use mapping number 3 requires that
C we also supply a version of CPMPXY to define that mapping (which see,
C below). Changing the value to 0 (to turn mapping off) or to 4 (to
C request the identity mapping) will make it clear that mapping 3 is as
C desired. (It produces nice round contours and the grid has jagged
C edges.)
C
CALL CPSETI ('MAP - MAPPING FLAG',3)
C
C Tell CONPACK to draw the edge of the grid.
C
CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1)
CALL CPSETI ('CLU - CONTOUR LEVEL USE FLAG',1)
C
C Initialize CONPACK.
C
CALL CPRECT (ZDAT,11,11,11,RWRK,2000,IWRK,1000)
C
C Draw 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
C
C PLEASE NOTE: The following version of CPMPXY is a modified copy of
C the current default version of CPMPXY and therefore has in it some
C stuff that wasn't required with version 3.1.1. Nevertheless, it
C should work with version 3.1.1. If you do use something like this
C with later versions of NCAR Graphics, though, it would be better to
C modify mapping 3 to do both the forward and reverse transformations.
C
SUBROUTINE CPMPXY (IMAP,XINP,YINP,XOTP,YOTP)
C
C By default, CONPACK draws contour lines using X coordinates in the
C range from XAT1 to XATM and Y coordinates in the range from YAT1 to
C YATN. Setting the flag 'MAP' non-zero causes those coordinates to be
C transformed by calling the subroutine CPMPXY. By default, 'MAP' = 1
C selects the EZMAP transformations and 'MAP' = 2 selects the polar
C coordinate transformations. The user of CONPACK may replace this
C subroutine as desired to transform the final coordinates and thus to
C transform the objects drawn.
C
C NOTE: As of 4/25/91, the default CPMPXY calls the new EZMAP routine
C MAPTRA instead of MAPTRN. The new routine returns 1.E12 for points
C which project outside the EZMAP perimeter.
C
C NOTE: As of 1/14/92, the default CPMPXY has been changed so that,
C when IMAP is negated, the inverse mapping is requested: (XINP,YINP)
C is a point in the current user coordinate system; (XOTP,YOTP) is
C returned and is the point which would be carried into (XINP,YINP) by
C the mapping numbered ABS(IMAP).
C
C An additional convention has been adopted which will allow CONPACK to
C find out whether a given inverse transformation is available. A call
C of the form
C
C CALL CPMPXY (0,REAL(IMAP),RFLG,DUM1,DUM2)
C
C will return information in RFLG about the mapping numbered IMAP, as
C follows:
C
C RFLG forward mapping defined inverse mapping defined
C ---- ----------------------- -----------------------
C 0. no no
C 1. yes no
C 2. no yes
C 3. yes yes
C
C Versions of CPMPXY that have not been updated to include these new
C features should continue to work for a period of time, but ought to
C be updated eventually.
C
C ---------------------------------------------------------------------
C
C Handle a request by the caller for information about the capabilities
C of this version of CPMPXY. Note that, if you modify CPMPXY to do
C other mappings, you should update the following code to correctly
C reflect the capabilities of the modified routine.
C
IF (IMAP.EQ.0) THEN ! MODIFIED CODE
IF (INT(XINP).EQ.1) THEN ! MODIFIED CODE
YINP=3. ! MODIFIED CODE
ELSE IF (INT(XINP).EQ.2) THEN ! MODIFIED CODE
YINP=3. ! MODIFIED CODE
ELSE IF (INT(XINP).EQ.3) THEN ! MODIFIED CODE
YINP=1. ! MODIFIED CODE
ELSE IF (INT(XINP).GE.4) THEN ! MODIFIED CODE
YINP=3. ! MODIFIED CODE
ELSE ! MODIFIED CODE
YINP=0. ! MODIFIED CODE
END IF ! MODIFIED CODE
C
C Handle the EZMAP case ...
C
ELSE IF (ABS(IMAP).EQ.1) THEN
IF (IMAP.GT.0) THEN
CALL MAPTRA (YINP,XINP,XOTP,YOTP)
IF (ICFELL('CPMPXY',1).NE.0) RETURN
ELSE
CALL MAPTRI (XINP,YINP,YOTP,XOTP)
IF (ICFELL('CPMPXY',2).NE.0) RETURN
END IF
C
C ... the polar coordinate case ...
C
ELSE IF (ABS(IMAP).EQ.2) THEN
IF (IMAP.GT.0) THEN
XOTP=XINP*COS(.017453292519943*YINP)
YOTP=XINP*SIN(.017453292519943*YINP)
ELSE
XOTP=SQRT(XINP*XINP+YINP*YINP)
YOTP=57.2957795130823*ATAN2(YINP,XINP)
END IF
C
C ... the staggered-grid case ...
C
ELSE IF (ABS(IMAP).EQ.3) THEN ! NEW CODE
IF (MOD(YINP+1.,2.).LE.1.) THEN ! NEW CODE
XOTP=XINP+.5*MOD(YINP+1.,2.) ! NEW CODE
ELSE ! NEW CODE
XOTP=XINP+1.-.5*MOD(YINP+1.,2.) ! NEW CODE
END IF ! NEW CODE
YOTP=YINP ! NEW CODE
C
C ... and everything else.
C
ELSE
XOTP=XINP
YOTP=YINP
END IF
C
C Done.
C
RETURN
C
END
This archive was generated by hypermail 2b29 : Wed Jun 28 2000 - 09:44:55 MDT