[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: 20020227: reading/writing McIDAS area files for use in GEMPAK



David,

There is a point of singularity in the transformation matrix
when rotating the globe. For example, when using ORT/40;-120;0,
the point -40;60 is the opposite side of the center of projection.
Although this point is not "viewable", the transformation can jump from
0 to 180, or 90 to -90 degrees rapidly as the point is approached
leading to the spurrious lines jumping across the globe.

I have attatched an updated version of:
$GEMPAK/source/gplt/transform/gmtowq.f

You can update your version of gplt by placing the attatched routine
in the gplt/transform directory and compiling with:

cd $GEMPAK/source/gplt/transform
make clean
make all

cd $GEMPAK/source/gplt/main
make clean
make all
make install
make clean

At this point, all programs that launch gplt from the command line will
be accessing the new gplt executable (you should make sure that you don't
have any gplt processes running at the time you install the new
executable).

To update the GUI's, you have to relink them since they include the
gplt.a library directly. You can do this with:

cd $GEMPAK/source/programs/gui/nmap2
make clean
make all
make install
make clean

cd $NAWIPS/nprogs
make clean
make all
make installl
make clean

cd $GARPHOME
make clean
make all
make install
make clean

I did some quick loops and things appear to look good.

Steve Chiswell
Unidata User Support



>
>
>
> >From: David Ovens <address@hidden>
> >Organization: UCAR/Unidata
> >Keywords: 200202280659.g1S6xYx22632
>
> >Steve,
> >
> >Turns out the orthographic projection is even buggier than I thought
> >with my last email.  I cannot produce a clean series of images with
> >any PROJ combination for any of the global GEMPAK files that I tried.
> >In
> >  http://www.atmos.washington.edu/~ovens/testimages/
> >you can check out these loops created with PROJ=ort/30;-135;0:
> >  loop_avn.html
> >  loop_eta104.html  -- no problems for this limited area data set
> >  loop_mrf3_m.html
> >  loop_nogaps.html
> >  loop_thin.html
> >  loop_ukmet.html
> >to see the sporadic spurious contour lines and the ever-present
> >strangeness of filled contours near the bottom of the globe.  Near the
> >bottom, it looks as though we are looking at contours on the other
> >side of the globe in a small wedge.
> >
> >David
> >
> >Unidata Support wrote:
> >>
> >>
> >> David,
> >>
> >> I looked at your gif inages in the testimages directory below.
> >> Can you also put a sample of your AREA file with the poked
> >> values in that directory so that I can look at the map drawing?
> >>
> >> I did some playing around with the existing projections that
> >> already exist in GEMPAK. I believe the Orthographic projection
> >> accomplishes what you originally described.
> >>
> >> For example, I created a globe-like plot using the following in GPMAP:
> >>  GAREA    = 0;-135;0;-135
> >>  PROJ     = ort/35;-135;0
> >>  LATLON   = 6/1/1/1;1/20;20/10;-135
> >>
> >>
> >> The latlon drawing is a little sensitive (eg buggy) when the longitude
> >> drawn on the back side of the earth is 180 degrees oposite the
> >> central longitude - so I picked the longitude increment of 20 (instead of 
> >> 15
> > )
> >> to avoid that problem.
> >>
> >> Steve Chiswell
> >> Unidata User Support
> >>
> >>
> >>
> >> >From: David Ovens <address@hidden>
> >> >Organization: UCAR/Unidata
> >> >Keywords: 200202130953.g1D9rkx03539
> >>
> >> >Steve and Tom,
> >> >
> >> >If you look in http://www.atmos.washington.edu/~ovens/testimages/
> >> >you will find 3 GIF images and 3 commented *lwu_pokes files that show
> >> >how the bogus McIDAS AREA file used in the corresponding plot was
> >> >generated.  Turns out that once I found the right words to POKE, this
> >> >was possible and both McIDAS-X and GEMPAK can handle it (for the most
> >> >part).  One thing I am wondering if you could still help on is the
> >> >apparent bug in the mapping.  Look for Greenland in the various plots
> >> >and you'll see that only avn3_2002021300_000.gif has its complete
> >> >outline.  Similar disappearances of the lat/lon lines also occur near
> >> >the pole.
> >> >
> >> >I realize that what I've done to generate the AREA file is pretty
> >> >klugey, but I still wonder if I can get trick GEMPAK into creating
> >> >avn1_2002021300_000.gif with a full map and lat/lon lines.  If I add
> >> >more bogus data to the file to correspond to the POKEd words (5, 6, 8,
> >> >and 9), would this make the maps and lat/lon labels around the pole
> >> >reappear?
> >> >
> >> >David
> >> >--
> >> >
> >> >David Ovens               e-mail: address@hidden
> >> >(206) 685-8108          plan: Real-time MM5 forecasting for Pacific 
> >> >Northwe
> > st
> >> >Research Meteorologist
> >> >Dept of Atmospheric Sciences, Box 351640
> >> >University of Washington
> >> >Seattle, WA  98195
> >> >
> >> >
> >> >Unidata Support wrote:
> >> >>
> >> >> >From: David Ovens <address@hidden>
> >> >> >Organization: University of Washington
> >> >> >Keywords: 200202121913.g1CJDNx21764 McIDAS AREA
> >> >>
> >> >> David,
> >> >>
> >> >> Chiz already responded to you on this from a GEMPAK perspective, but
> >> >> I figured I would throw in my 2 cents worth on the McIDAS side.
> >> >>
> >> >> >I am hoping to create a bogus McIdas area file -- basically a GOES-10
> >> >> >image rotated up so Seattle is about at the center of the image -- so
> >> >> >that I can get a different satellite projection in GEMPAK (if you want
> >> >> >an example of what my output would look like, see
> >> >> > http://www.atmos.washington.edu/~ovens/loops/wxloop.cgi?joes_okc_h500+a
> > ll
> >> >> >).  Anyhow, I think that all I need to do is figure out how to read an
> >> >> >existing McIdas file (GOES-10) and write it back out changing only the
> >> >> >4th word in the navigation section from "102042" to "474500" or
> >> >> >something like that.  Is there FORTRAN (preferably), C, or PERL code
> >> >> >out there that could help me with this?
> >> >>
> >> >> This can easily be done in McIDAS.  You can use the LWU command to list
> >> >> out and change values in any file.  So, for instance, let's suppose
> >> >> that the image you want to play around with is in file AREA1234.
> >> >> Use LWU from a McIDAS session to list out the first 64 words of the
> >> >> file:
> >> >>
> >> >> LWU LIST AREA1234 0 63
> >> >>
> >> >> Now, suppose you wanted to change the value of word 123 to 474500:
> >> >>
> >> >> LWU POKE AREA1234 474500 123
> >> >>
> >> >> Remember that the index into the file that LWU will use is zero based.
> >> >>
> >> >> As a side comment, the GVAR navigation model is so complex that I don't
> >> >> think that changing the value of one word in the NAV block will do the
> >> >> rotation that you want.  But hey, I have never tried this so go for
> >> >> it!
> >> >>
> >> >> After changing a NAV value, you will need to redisplay the image
> >> >> and then draw a map on it to see the effect:
> >> >>
> >> >> IMGDISP MYDATA/IMAGES.1234 REFRESH='EG;MAP'
> >> >>
> >> >> If this fails, it means that you can not simply change one NAV parameter
> >> >> in the GVAR navigation and be successful.
> >> >>
> >> >> One navigation that does lend itself to simple rotations in McIDAS is
> >> >> that for METEOSAT.  The subsatellite longitude is called out explicitly
> >> >> in METEOSAT imagery.  Changing it effectively rotates the image to
> >> >> the new longitude.  What is not called out, however, is the subsatellite
> >> >> latitude as that is assumed to be the equator.  The only problem for
> >> >> you using a METEOSAT image is that GEMPAK does not support METEOSAT
> >> >> navigation.
> >> >>
> >> >> >By the way, I have looked at maknav and imgremap in McIdas-X and I
> >> >> >understand the LAMB, MERC, PS, RADR, RECT, and MOLL projections, it's
> >> >> >possible that SIN, TANC, or NOAA does what I want, but I don't
> >> >> >understand those.  Basically, I didn't see how I could do a new
> >> >> >satellite projection if I did not already have a file with the
> >> >> >projection that I want to create.
> >> >>
> >> >> IMGREMAP can be used to write image data into a new AREA file
> >> >> specifying the projection that you want.  The projections
> >> >> supported by IMGERMAP for new output files are:
> >> >>
> >> >>      =LAMB slat1 slat2 slon | Lambert Conformal, standard latitude
> >> >>                               and standard longitude (slat1 def=30,
> >> >>                               slat2 def=50, slon def=center longitude)
> >> >>      =MERC slat | Mercator projection and standard latitude (slat def=0)
> >> >>      =MOLL slon | Mollweide projection and standard longitude (slon 
> >> >> def=0
> > )
> >> >>      =PS slat slon | Polar Stereographic projection, standard latitude
> >> >>                      and standard longitude
> >> >>                      (slat def=60, slon def=center longitude)
> >> >>      =RADAR rot | Radar projection and rotation angle (rot def=0)
> >> >>      =RECT | Rectilinear projection
> >> >>      =SIN | Sinusoidal Equal Area projection
> >> >>      =TANC slat slon | Tangent Cone projection, standard latitude and
> >> >>
> >> >> If you want the navigation of the resultant image to be that for GVAR,
> >> >> then your only option is to remap into an existing image.
> >> >>
> >> >> I am curious to hear the results of your changing NAV parameters in
> >> >> a GVAR image!
> >> >>
> >> >> >Thanks for any help or suggestions.
> >> >> >
> >> >> >David
> >> >> >David Ovens            e-mail: address@hidden
> >> >> >(206) 685-8108          plan: Real-time MM5 forecasting for Pacific 
> >> >> >Nort
> > hwe
> >> > st
> >> >> >Research Meteorologist
> >> >> >Dept of Atmospheric Sciences, Box 351640
> >> >> >University of Washington
> >> >> >Seattle, WA  98195
> >> >>
> >> >> Tom Yoksas
> >> >> *************************************************************************
> > ***
> >> >> Unidata User Support                                    UCAR Unidata 
> >> >> Prog
> > ram
> >> >> (303)497-8644                                                  P.O. Box 
> >> >> 3
> > 000
> >> >> address@hidden                                   Boulder, CO 80
> > 307
> >> >> -------------------------------------------------------------------------
> > ---
> >> >> Unidata WWW Service                        http://www.unidata.ucar.edu/
> >
> >> >> *************************************************************************
> > ***
> >> >>
> >> >
> >> >
> >> >From address@hidden Tue Feb 12 15:21:48 2002
> >> >Received: from blizzard.atmos.washington.edu 
> >> >(blizzard.atmos.washington.edu
> >  [1
> >> > 28.95.175.12])
> >> >  by unidata.ucar.edu (UCAR/Unidata) with ESMTP id g1CMLlx08683
> >> >  for <address@hidden>; Tue, 12 Feb 2002 15:21:47 -0700 (MST)
> >> >Organization: UCAR/Unidata
> >> >Keywords: 200202122221.g1CMLlx08683
> >> >Received: (from ovens@localhost)
> >> >  by blizzard.atmos.washington.edu (8.11.6+Sun/8.11.6) id g1CMLld06158
> >> >  for address@hidden; Tue, 12 Feb 2002 14:21:47 -0800 (PST)
> >> >From: David Ovens <address@hidden>
> >> >Message-Id: <address@hidden>
> >> >Subject: Re: 20020212: sat projection in GEMPAK
> >> >To: address@hidden
> >> >Date: Tue, 12 Feb 2002 14:21:46 -0800 (PST)
> >> >In-Reply-To: <address@hidden> from "Unidata Sup
> > por
> >> > t" at Feb 12, 2002 01:59:24 PM
> >> >X-Mailer: ELM [version 2.5 PL2]
> >> >MIME-Version: 1.0
> >> >Content-Type: text/plain; charset=us-ascii
> >> >Content-Transfer-Encoding: 7bit
> >> >
> >> >Steve,
> >> >
> >> >Thanks for the idea, but the professor for whom I doing this really
> >> >wants that circular, satellite view that NCARGraphics is able to
> >> >provide.  I can still use the NCARGraphics code, but it takes about
> >> >5-10 times as long to run and is prone to failing every week or two
> >> >(something to do with attempting to do contours at the limb, I
> >> >think).  If I can get a fake satellite image, then GEMPAK can do it
> >> >for me.
> >> >
> >> >Too bad GEMPAK doesn't have the ability to do SAT projections all by
> >> >itself.
> >> >
> >> >I will come up with a solution and let you and Tom know what it is.
> >> >
> >> >David
> >> >
> >> >
> >> >Unidata Support wrote:
> >> >>
> >> >> >From: Unidata User Support <address@hidden>
> >> >> >Organization: Unidata Program Center/UCAR
> >> >> >Keywords:
> >> >>
> >> >>
> >> >> David,
> >> >>
> >> >> As I figured, Tom was able to suggest an easy way to create the
> >> >> satellite projection for an AREA file. As an asside, if you weren't 
> >> >> tryin
> > g
> >> >> to project a satellite image, but rather just wanted to give the
> >> >> illusion of the disk then you could use other GEMPAK projections
> >> >> such as rotated cylindrical projections (depending on how much you 
> >> >> wanted
> >  to
> >> >> zoom in on Washington). Such as:
> >> >> GAREA    = -10;-180;20;-30
> >> >> PROJ     = ced/45;-120;90
> >> >>
> >> >>
> >> >> Steve CHiswell
> >> >> *************************************************************************
> > ***
> >> >> Unidata User Support                                    UCAR Unidata 
> >> >> Prog
> > ram
> >> >> (303)497-8644                                                  P.O. Box 
> >> >> 3
> > 000
> >> >> address@hidden                                   Boulder, CO 80
> > 307
> >> >> -------------------------------------------------------------------------
> > ---
> >> >> Unidata WWW Service                        http://www.unidata.ucar.edu/
> >
> >> >> *************************************************************************
> > ***
> >> >>
> >> >
> >> >
> >> >--
> >> >
> >> >David Ovens               e-mail: address@hidden
> >> >(206) 685-8108          plan: Real-time MM5 forecasting for Pacific 
> >> >Northwe
> > st
> >> >Research Meteorologist
> >> >Dept of Atmospheric Sciences, Box 351640
> >> >University of Washington
> >> >Seattle, WA  98195
> >> >
> >>
> >> ****************************************************************************
> >> Unidata User Support                                    UCAR Unidata 
> >> Program
> >> (303)497-8644                                                  P.O. Box 
> >> 3000
> >> address@hidden                                   Boulder, CO 80307
> >> ----------------------------------------------------------------------------
> >> Unidata WWW Service                        http://www.unidata.ucar.edu/
> >> ****************************************************************************
> >>
> >
> >
> >--
> >
> >David Ovens          e-mail: address@hidden
> >(206) 685-8108          plan: Real-time MM5 forecasting for Pacific Northwest
> >Research Meteorologist
> >Dept of Atmospheric Sciences, Box 351640
> >University of Washington
> >Seattle, WA  98195
> >
>
>
        SUBROUTINE GMTOWQ ( amtrx, invers, polnr,
     +                      np, rlat, rlon, slat, slon, iret )
C************************************************************************
C* GMTOWQ                                                               *
C*                                                                      *
C* This subroutine converts a location from actual latitude/longitude   *
C* coordinates to latitude and longitude in a rotated coordinate, if    *
C* INVERS = .false.  The inverse transformation is done if INVERS =     *
C* .true.  Actual lat/lon coordinates are M coordinates; the lat/lon    *
C* coordinates in the rotated frame are W coordinates (map) or Q        *
C* coordinates (grid).                                                  *
C*                                                                      * 
C* This routine performs the transformation needed for the any map/grid *
C* projection.  The rotation matrix is AMTRX.  The rotation matrices    *
C* for map and grid projections are set in the update routines UPD for  *
C* map/grid projections.                                                *
C*                                                                      *
C* The transformation from M -> W/Q proceeds as follows:                *
C*                                                                      *
C*   1)  Transform actual lat/lon to cartesian coordinates.             *
C*                                                                      *
C*   2)  Rotate the cartesian system three times:                       *
C*       a) Through POLON = ANGLE2 to reposition xy axes                *
C*       b) Through POLAT = ANGLE1 to reposition xz axes                *
C*       c) Through ROTAT = ANGLE3 to reposition yz or xy axes          *
C*                                                                      *
C*   3)  Transform from rotated cartesian to spherical W/Q coordinates. *
C*                                                                      *
C* The inverse from W/Q -> M reverses these steps.                      *
C*                                                                      *
C* The longitude rotation angle is only used for the special case when  *
C* a point falls on a pole in the rotated coordinate.  In that case,    *
C* longitude is indefinite, but for cylindrical projections a definite  *
C* longitude is required when a pole corresponds to a corner point      *
C* That is why this rotation is done.                                   *
C*                                                                      *
C*                                                                      *
C* GMTOWQ ( AMTRX, INVERS, POLNR, NP, RLAT, RLON, SLAT, SLON, IRET )    *
C*                                                                      *
C* Input parameters:                                                    *
C*      AMTRX           REAL            Transformation matrix           *
C*      INVERS          LOGICAL         Inverse transform flag          *
C*      POLNR           REAL            Longitude rotation angle        *
C*      NP              INTEGER         Number of points                *
C*      RLAT (NP)       REAL            Input latitudes                 *
C*      RLON (NP)       REAL            Input longitudes                *
C*                                                                      *
C* Output parameters:                                                   *
C*      SLAT (NP)       REAL            Output latitudes                *
C*      SLON (NP)       REAL            Output longitudes               *
C*      IRET            INTEGER         Return code                     *
C**                                                                     *
C* Log:                                                                 *
C* K. Brill/EMC          3/96                                           *
C* T. Lee/GSC            7/96           Fixed pole points plot for RADAR*
C* S. Chiswell/Unidata   2/02           Check for SLAT~=-90 in ATAN2    *
C************************************************************************
        INCLUDE         'ERROR.PRM'
        INCLUDE         'GEMPRM.PRM'
        PARAMETER       ( EPS = .03 )
C*
        REAL            amtrx (3,3)
        LOGICAL         invers
        REAL            rlat (*), rlon (*), slat (*), slon (*)
C*
        REAL            tmtrx (3,3), x (3), xp (3), phi,
     +                  rlm, cosp, xx, yy
C*
        INCLUDE         'ERMISS.FNC'
C------------------------------------------------------------------------
        iret = NORMAL
C
C*      Set the transformation matrix and look for identity matrix.
C
        IF ( .not. invers ) THEN
            DO j = 1, 3
                DO i = 1, 3
                    tmtrx (i,j) = amtrx (i,j)
                END DO
            END DO
        ELSE
C
C*          The inverse is just the transpose.
C
            DO j = 1, 3
                DO i = 1, 3
                    tmtrx (i,j) = amtrx (j,i)
                END DO
            END DO
        END IF
C
C*      Process all the points.
C
        DO ii = 1, np
            IF ( .not. ERMISS ( rlat (ii) ) .and.
     +           .not. ERMISS ( rlon (ii) ) ) THEN
C
C*              Convert from spherical to cartesian coordinates.
C
                r1 = rlat(ii)
                r2 = rlon(ii)
                phi   = rlat (ii) * DTR
                rlm   = rlon (ii) * DTR
                cosp  = COS (phi)
                xx    = COS (rlm)
                yy    = SIN (rlm)
                x (1) = cosp * xx
                x (2) = cosp * yy
                x (3) = SIN (phi)

C
C*              Rotate the coordinates and confine coordinate values.
C
                DO i = 1, 3
                    xp (i) = 0.
                    DO j = 1, 3
                        xp (i) = xp (i) + tmtrx (i,j) * x (j)
                    END DO
                    IF ( xp (i) .gt. 1. ) xp (i) = 1.
                    IF ( xp (i) .lt. -1. ) xp (i) = -1.
                END DO
C
C*              Convert from cartesian to spherical coordinates.
C*              If the rotated coordinate is at a pole; then,
C*              rotate the longitude (since it is indefinite).
C
                slat (ii) = ASIN ( xp (3) ) * RTD
                IF ( ( ABS ( slat (ii) - 90. ) .lt. EPS ) .or.
     +               ( ABS ( slat (ii) + 90. ) .lt. EPS ) ) THEN
                    IF ( slat(ii) .gt. 0 ) THEN
                       slat (ii) = 90.
                    ELSE
                       slat (ii) = -90.
                    END IF
                    rcos = COS ( polnr )
                    rsin = SIN ( polnr )
                    IF ( invers ) rsin = -rsin
                    xxp = xx * rcos + yy * rsin
                    yyp = -xx * rsin + yy * rcos
                    slon (ii) = ATAN2 ( yyp, xxp ) * RTD
                    if(slat(ii) .eq. -90) slon(ii) = ATAN2 ( yyp, -xxp ) * RTD
                ELSE
                    slon (ii) = ATAN2 ( xp (2), xp (1) ) * RTD
                END IF
C
C*              Check for the singular point in transformation matrix where
C*              either xxp or yyp become very large. This point is generally 
hidden,
C*              but prevents line from zipping across to oposite side of the 
image.
C
                splon = slon (ii) * DTR
                IF ( ( slat(ii) .lt. -89. ) .and.
     +             ( ( ABS( SIN (splon) ) .lt. EPS ) .or.
     +               ( ABS( COS (splon) ) .lt. EPS ) ) ) THEN
                    slat (ii) = RMISSD
                    slon (ii) = RMISSD
                END IF
            ELSE
                slat (ii) = RMISSD
                slon (ii) = RMISSD
            END IF
        END DO
C*
        RETURN
        END