[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
20020131: radius of earth bug in imar2gm.f -- second notice
- Subject: 20020131: radius of earth bug in imar2gm.f -- second notice
- Date: Thu, 31 Jan 2002 16:18:05 -0700
David,
I updated the imar2gm.f recently to add the TANC projection as well as
fix an outstating bug in PS. The PS section I fixed does use the
radius of earth from the MciDAS file. I can make the other RADIUS to re
changes. There could still be some discrepancy in the lack of eccentricity
use however.
Steve Chiswell
Unidata User SUpport
>From: David Ovens <address@hidden>
>Organization: UCAR/Unidata
>Keywords: 200201312147.g0VLlkx05878
>GEMPAK Support,
>
>I was disappointed to see that my bug fixes from Dec. 2000 for
>gemlib/im/imar2gm.f have not fully made it into GEMPAK 5.6.E.1 or
>previous GEMPAK versions. The crux of our differences now are:
>
> 1) in the actual use of the radius of the earth that
> is in the McIDAS file (re) vs GEMPAK's value (RADIUS), and
> 2) the location of re in the ianav array for projections other than
> 'MERC' -- I have re = ianav(8) for everything but 'MERC' where I
> have ianav(7), whereas you have re = ianav(7) for everything.
>
>Unless, McIDAS code has changed, I think your version is still in
>error. Below you will find an excerpt from my earlier email and
>below that my newest version of imar2gm.f.
>
>Please let me know what you think, I'd like to get this incorporated
>into the standard release. This is the only way that my McIDAS area
>files of topography jibe with the map database.
>
>David
>--
>
>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
>
>
>>From ovens Thu Dec 28 06:58:44 2000
>Subject: I found and fixed a bug in GEMPAK imar2gm.f
>To: address@hidden, justin (Justin Sharp), harry (Harry Edmon)
>Date: Thu, 28 Dec 2000 06:58:44 -0800 (PST)
>X-Mailer: ELM [version 2.5 PL2]
>Content-Length: 16228
>Status: OR
>
>Hello,
>
>I've been attempting to modify Justin Sharp's excellent program that
>makes McIdas area files out of 30-second terrain data. (A program that
>he has placed on CD with the terrain data set). It has always worked
>fine with the RECTilinear projection of McIdas, but when I tried to
>use IMGREMAP as in the following code extract
> $imgremap = "imgremap.k MYDATA/IMAGES.$area MYDATA/IMAGES.9998 PRO=MERC LATL
> ON=$center_lat $center_lon RES=1 SIZE=ALL";
>I was having trouble with the GEMPAK map being offset with the
>terrain. In McIdas, however, I could not discern any problem with the
>remapped file.
>
>Well, I have gone to the McIdas source and to the GEMPAK source and I
>have been able to fix the problem in GEMPAK-5.6.A by correcting a
>discrepancy between the Earth radius used in GEMPAK and that used in
>the McIdas file. I have modified $GEMPAK/source/gemlib/im/imar2gm.f
>as you will see below to replace RADIUS with re (obtained from the
>navigation header of the file).
>
>I have tested this with MERC, LAMB, and PS terrain/McIdas files and it
>all seems to line up fine using sfmap.
>
>David
>========================= imar2gm.f ====================
> SUBROUTINE IM_AR2GM ( imgfil, rnvblk, iret )
>C************************************************************************
>C* IM_AR2GM *
>C* *
>C* This subroutine reads the header information from a MCIDAS AREA file
> *
>C* and sets the navigation *
>C* *
>C* IM_AR2GM ( IMGFIL, RNVBLK, IRET ) *
>C* *
>C* Input parameters: *
>C* IMGFIL CHAR* Image file *
>C* *
>C* Output parameters: *
>C* RNVBLK (LLNNAV) REAL Navigation block *
>C* IRET INTEGER Return code *
>C* 0 = normal return *
>C* -2 = Error opening/reading file*
>C* -3 = Invalid image file format *
>C* -4 = Invalid image file header *
>C* -5 = Invalid image navigation *
>C** *
>C* Log:
> *
>C* J. Cowie/COMET 2/95 *
>C* C. Lin/EAI 6/95 Change input from file name to lunmf *
>C* Change GG_* to IM_* for projection setup*
>C* J. Cowie/COMET 10/95 Set xy image scaling to 1.0 *
>C* J. Cowie/COMET 4/96 Add calibration type, check for 'RAW ' *
>C* T. Lee/GSC 7/96 Renamed IM_ARHD; Cleaned up & eliminated*
>C* map projection calls *
>C* D. Keiser/GSC 8/96 Added RETURN after call to GSATMG *
>C* T. Lee/GSC 9/96 Used RADIUS and floating point *
>C* D.W.Plummer/NCEP 10/96 Fix RADR processing bug, rem to RADIUS *
>C* T. Lee/GSC 11/96 Declared timfil; Used meter everywhere *
>C* J. Cowie/COMET 1/97 Changed IMGDEF common variables names *
>C* S. Chiswell/UNIDATA 1/97 Add extra scaling chck for RECT images *
>C* J. Cowie/COMET 2/97 Fixed RADR off-by-one error *
>C* J. Cowie/COMET 12/97 Changed jyear to 1900 + current year *
>C* J. Cowie/COMET 12/97 Added imradf, removed immode, imraw *
>C* S. Danz/AWC 05/99 Fix to RECT off-by-one error *
>C* S. Jacobs/NCEP 6/00 Used correct image name in GSATMG call *
>C* S. Danz/AWC 10/00 Fix RECT setting center longitude *
>C* J. Cowie/COMET 12/00 Added MCRADR projection type *
>C* D. Ovens/UW 12/00 Fix MERC,LAMB,PS Earth Radius mismatches*
>C* D.W.Plummer/NCEP 9/01 Set imbswp=0 when no byte swapping *
>C* S. Chiswell/Unidata 10/01 Fix PS dlat calculation *
>C************************************************************************
> INCLUDE 'GEMPRM.PRM'
> INCLUDE 'IMGDEF.CMN'
> INCLUDE 'AREAFL.CMN'
>C*
> CHARACTER*(*) imgfil
> REAL rnvblk (*)
>C*
> CHARACTER navtyp*4, proj*3
> INTEGER iarray (704), idtarr (3)
>C*
> EQUIVALENCE (iarray (1), iadir (1))
> EQUIVALENCE (iarray (65), ianav (1))
>C------------------------------------------------------------------------
> iret = 0
> sign = 1.
>C
>C* Open the file and read AREA DIR and NAV blocks
>C
> CALL FL_DOPN ( imgfil, 704, .false., lunmf, iret )
> CALL FL_READ ( lunmf, 1, 704, iarray, ier )
> IF ( ier .ne. 0 ) THEN
> iret = -1
> RETURN
> END IF
> CALL FL_CLOS ( lunmf, ier )
>C
>C* Check to see if the header items needs byte-swapping.
>C
> IF ( iadir (2) .ne. 4 ) THEN
>C
> imbswp = 1
>C
>C* Swap bytes in AREA DIR except comment.
>C
> ier = MV_SWP4 ( 24, iadir ( 1 ), iadir ( 1 ) )
> ier = MV_SWP4 ( 19, iadir ( 33 ), iadir ( 33 ) )
> ier = MV_SWP4 ( 3, iadir ( 54 ), iadir ( 54 ) )
> ier = MV_SWP4 ( 7, iadir ( 58 ), iadir ( 58 ) )
>C
>C* Swap bytes in NAV block except 1st word (nav label).
>C
> ier = MV_SWP4 ( 640-1, ianav ( 2 ) , ianav ( 2 ) )
>C
> ELSE
>C
> imbswp = 0
>C
> END IF
>C
>C* Determine calibration type Calibrations: 'RAW'
>C* 'BRIT'
>C* 'TEMP'
>C
> CALL ST_ITOC ( iadir (53), 1, cmcalb, ier )
>C
>C* Parse AREA data header. Check for radar image.
>C
> imsorc = iadir ( 3 )
> IF ( imsorc .eq. 7 ) then
> imradf = 1
> imtype = iadir ( 22 )
> ELSE
> imradf = 0
> imtype = iadir ( 19 )
> END IF
>C
> jyear = ( iadir ( 4 ) / 1000 ) + 1900
> jday = MOD ( iadir ( 4 ), 1000 )
> CALL TI_JTOI ( jyear, jday, idtarr, ier )
> imdate = idtarr ( 1 ) * 10000 + idtarr ( 2 ) * 100 + idtarr ( 3 )
> imtime = iadir ( 5 )
> iullne = iadir ( 6 )
> iulele = iadir ( 7 )
> imnlin = iadir ( 9 )
> imnpix = iadir ( 10 )
> imdpth = iadir ( 11 )
> rmyres = iadir ( 12 )
> rmxres = iadir ( 13 )
> imnchl = iadir ( 14 )
> imprsz = iadir ( 15 )
> imdcsz = iadir ( 49 )
> imclsz = iadir ( 50 )
> imlvsz = iadir ( 51 )
> imvald = iadir ( 36 )
> icaloff = iadir ( 63 )
> if(icaloff.gt.0) CALL IM_CALIB(imgfil,icaloff,iret)
> rmxysc = 1.0
>C
>C* Determine offset to start of data
>C
> IF ( iadir ( 34 ) .gt. 0 ) THEN
> imdoff = iadir ( 34 )
> ELSE
> iret = -4
> RETURN
> END IF
> imldat = imnlin * ( imnpix * imdpth * imnchl + imprsz )
>C
>C* Set image subset to full image for now.
>C
> imleft = 1
> imtop = 1
> imrght = imnpix
> imbot = imnlin
>C
>C* Set navigation.
>C
>C* Standard latitude.
>C
> angdd = ianav ( 4 ) / 10000
> angmm = ianav ( 4 ) / 100 - angdd * 100
> angss = ianav ( 4 ) - angdd * 10000 - angmm * 100
> stdlat = angdd + angmm / 60. + angss / 3600.
> phi0 = stdlat
> phi0r = phi0 * DTR
> IF ( phi0r .lt. 0. ) sign = - 1.
>C
>C* Normal longitude. If west positive, make west negative.
>C
> angdd = ianav ( 6 ) / 10000
> angmm = ianav ( 6 ) / 100 - angdd * 100
> angss = ianav ( 6 ) - angdd * 10000 - angmm * 100
> clon = angdd + angmm / 60. + angss / 3600.
> IF ( ianav ( 10 ) .ge. 0 ) clon = - clon
>C
>C* Set pixel/grid spacing and earth radius.
>C
> dx = ianav ( 5 ) * iadir ( 13 )
>cdo re = ianav ( 7 )
> re = ianav ( 8 )
> if (re .le. 6200000.) re = RADIUS
>C
>C* Set pixel/grid dimension.
>C
> rgx = imnpix
> rgy = imnlin
>C
>C* Location of pole point (rxp, ryp); (1,1) at lower-left corner.
>C
> rxp = float ( ianav ( 3 ) - iadir ( 7 ) ) / iadir ( 13 ) + 1.
> ryp = rgy - float ( ianav ( 2 ) - iadir ( 6 ) ) / iadir ( 12 )
>C
> CALL ST_ITOC ( ianav, 1, navtyp, ier )
>C
>C* ******************************************
>C* *** Satellite projection (MCIDAS GOES) ***
>C* ******************************************
>C
> IF ( navtyp .eq. 'GOES' .or. navtyp .eq. 'GVAR' .or.
> + navtyp .eq. 'RADR') THEN
>C
> CALL GSATMG ( imgfil, iadir, ianav, imleft, imtop, imrght,
> + imbot, ier )
> RETURN
>C
>C* ***************************
>C* *** Polar Stereographic ***
>C* ***************************
>C
> ELSE IF ( navtyp .eq. 'PS' ) THEN
>C
>C* Compute lat/lon of the lower-left corner.
>C
>
>C squished projection, need (12) LINERES
> dy = ianav ( 5 ) * rmyres
>
>
> dxp = 1. - rxp
> dyp = 1. - ryp
> alpha = 1. + SIN ( ABS ( phi0r ) )
> rm = dy * SQRT ( dxp * dxp + dyp * dyp ) / alpha
>cdo dlat1 = sign * ( HALFPI - 2. * ATAN ( rm / RADIUS ) ) * RTD
> dlat1 = sign * ( HALFPI - 2. * ATAN ( rm / re ) ) * RTD
> IF ( dyp .ne. 0. ) THEN
> dyp = - dyp * sign
> thta = ATAN2 ( dxp, dyp ) * RTD
> dlon1 = clon + thta
> CALL PRNLON ( 1, dlon1, ier )
> ELSE
> dlon1 = clon
> END IF
>C
>C* Compute lat/lon of the upper-right corner.
>C
> dxp = rgx - rxp
> dyp = rgy - ryp
> rm = dy * SQRT ( dxp * dxp + dyp * dyp ) / alpha
>cdo dlat2 = sign * ( HALFPI - 2. * ATAN ( rm / RADIUS ) ) * RTD
> dlat2 = sign * ( HALFPI - 2. * ATAN ( rm / re ) ) * RTD
> IF ( dyp .ne. 0. ) THEN
> dyp = - dyp * sign
> thta = ATAN2 ( dxp, dyp ) * RTD
> dlon2 = clon + thta
> CALL PRNLON ( 1, dlon2, iret )
> ELSE
> dlon2 = clon
> END IF
>C
>C* Set map projection angles. Angle 2 is the central longitude.
>C
> IF ( stdlat .ge. 0. ) THEN
> dangl1 = 90.
> ELSE
> dangl1 = - 90.
> END IF
> dangl2 = clon
> dangl3 = 0.
> proj = 'STR'
>C
>C* ****************
>C* *** Mercator ***
>C* ****************
>C
> ELSE IF ( navtyp .eq. 'MERC' ) THEN
>C
>C* Reset earth radius and set eccentricity.
>C
> re = ianav ( 7 )
> if (re .le. 6200000.) re = RADIUS
> ecc = ianav ( 8 ) / 1000000.
>C
>C* Compute lat/lon of the lower-left corner.
>C
> dxp = 1. - rxp
> dyp = 1. - ryp
> rm = dx * dyp
>cdo rcos = RADIUS * COS ( phi0r )
> rcos = re * COS ( phi0r )
> arg = EXP ( rm / rcos )
> dlat1 = ( 2. * ATAN ( arg ) - HALFPI ) * RTD
> dlon1 = clon + ( dx * dxp / rcos ) * RTD
> CALL PRNLON ( 1, dlon1, ier )
>C
>C* Compute lat/lon of the upper-right corner point.
>C
> dxp = rgx - rxp
> dyp = rgy - ryp
> rm = dx * dyp
> arg = EXP ( rm / rcos )
> dlat2 = ( 2. * ATAN ( arg ) - HALFPI ) * RTD
> dlon2 = clon + ( dx * dxp / rcos ) * RTD
> CALL PRNLON ( 1, dlon2, ier )
>C
>C* Set the map projection and angles.
>C
> dangl1 = 0.
> dangl2 = clon
> dangl3 = 0.
> proj = 'MER'
>C
>C* *************************
>C* *** Lambert Conformal ***
>C* *************************
>C
> ELSE IF ( navtyp .eq. 'LAMB' ) THEN
>C
>C* Standard latitude #1
C
> std1 = stdlat
>C
>C* Standard latitude #2
>C
> angdd = ianav ( 5 ) / 10000
> angmm = ianav ( 5 ) / 100 - angdd * 100
> angss = ianav ( 5 ) - angdd * 10000 - angmm * 100
> std2 = angdd + angmm / 60. + angss / 3600.
>C
>C* Normal longitude. If west positive, make west negative.
>C
> angdd = ianav ( 7 ) / 10000
> angmm = ianav ( 7 ) / 100 - angdd * 100
> angss = ianav ( 7 ) - angdd * 10000 - angmm * 100
> clon = angdd + angmm / 60. + angss / 3600.
> IF ( ianav ( 10 ) .ge. 0 ) clon = - clon
>C
>C* Compute pixel/grid spacing and colatitudes.
>C
> dx = ianav ( 6 ) * iadir ( 13 )
> psi1 = HALFPI - ABS ( std1 ) * DTR
> psi2 = HALFPI - ABS ( std2 ) * DTR
>C
>C* Compute cone constant.
>C
> IF ( psi1 .eq. psi2 ) THEN
> ccone = COS ( psi1 )
> ELSE
> tmp1 = ALOG ( SIN ( psi2 ) / SIN ( psi1 ) )
> tmp2 = ALOG ( TAN ( psi2 / 2. ) / TAN ( psi1 / 2. ) )
> ccone = tmp1 / tmp2
> END IF
>C
>C* Compute lat/lon of the lower-left corner. Sign = 1/-1
>C* denotes NH/SH.
>C
> dxp = 1. - rxp
> dyp = 1. - ryp
> rm = dx * SQRT ( dxp * dxp + dyp * dyp )
>cdo tmp = ccone / ( RADIUS * SIN ( psi1 ) )
> tmp = ccone / ( re * SIN ( psi1 ) )
> arg = ( rm * tmp ) ** ( 1. / ccone ) * TAN ( psi1 / 2. )
> dlat1 = sign * ( HALFPI - 2. * ATAN ( arg ) ) * RTD
>C
> IF ( dyp .ne. 0. ) THEN
> dyp = - dyp
> thta = ATAN2 ( dxp, dyp ) * RTD / ccone
> dlon1 = clon + thta
> CALL PRNLON ( 1, dlon1, ier )
> ELSE
> dlon1 = clon
> END IF
>C
>C* Compute lat/lon of the upper-right corner.
>C
> dxp = rgx - rxp
> dyp = rgy - ryp
> rm = dx * SQRT ( dxp * dxp + dyp * dyp )
> arg = ( rm * tmp ) ** ( 1. / ccone ) * TAN ( psi1 / 2. )
> dlat2 = sign * ( HALFPI - 2. * ATAN ( arg ) ) * RTD
> IF ( dyp .ne. 0. ) THEN
> dyp = - dyp
> thta = ATAN2 ( dxp, dyp ) * RTD / ccone
> dlon2 = clon + thta
> CALL PRNLON ( 1, dlon2, ier )
> ELSE
> dlon2 = clon
> END IF
>C
>C* Northern or Southern Hemisphere projection
>C
> IF ( ( std1 + std2 ) .gt. 0. ) THEN
> proj = 'LCC'
> ELSE
> proj = 'SCC'
> END IF
> dangl1 = std1
> dangl2 = clon
> dangl3 = std2
>C
>C
>C* *******************
>C* *** Rectilinear ***
>C* *******************
>C
> ELSE IF ( navtyp .eq. 'RECT' ) THEN
>C
>C* Row/column and lat/lon of reference point.
>C
> rrow = ianav ( 2 )
> rcol = ianav ( 4 )
> dlat = ianav ( 3 ) / 10000.
> dlon = ianav ( 5 ) / 10000.
>C
>C* Determine longitude convention.
>C
> IF ( ianav ( 11 ) .ge. 0 ) sign = - 1.
> dlon = sign * dlon
>C
>C* Lat/lon resolution (degrees/imageline).
>C* Check for scaling factor. If it is not present use the
>C* default scaling.
>C
> IF ( ianav ( 14 ) .eq. 0 ) THEN
> reslat = ianav ( 6 ) / 10000.
> ELSE
> reslat = ianav ( 6 ) / ( 10. ** ianav (14) )
> END IF
> IF ( ianav ( 15 ) .eq. 0 ) THEN
> reslon = ianav ( 7 ) / 10000.
> ELSE
> reslon = ianav ( 7 ) / ( 10. ** ianav (15) )
> END IF
>C
>C* Determine upper-left lat/lon
>C* (McIDAS image full coordinates 1,1).
>C
> uldlat = dlat + reslat * ( rrow - iullne )
> uldlon = dlon + reslon * sign * ( rcol - iulele )
>C
>C* Lat/lon for lower-left corner.
>C
> dlat1 = uldlat - ( reslat * (imnlin * rmyres - 1))
> dlon1 = uldlon
>C
>C* Lat/lon for upper-right corner.
>C
> dlat2 = uldlat
> dlon2 = uldlon - sign * ( reslon * (imnpix * rmxres - 1) )
>C
>C* Set angle1 to the lon/lat aspect ratio.
>C* Set angle2 to the central longitude.
>C
> dangl1 = reslon / reslat
> dangl2 = uldlon -
> + sign * ( reslon * ( (imnpix * rmxres) / 2 - 1 ) )
> dangl3 = 0.
> proj = 'MCD'
>C
>C* Not valid
>C
> ELSE
> iret = -5
> RETURN
> END IF
>C
>C* Set navigation blocks.
>C
> CALL GR_MNAV ( proj, imnpix, imnlin, dlat1, dlon1, dlat2,
> + dlon2, dangl1, dangl2, dangl3, .true.,
> + rnvblk, ier )
>C
>C* Set map projection.
>C
> CALL GSMPRJ ( proj, dangl1, dangl2, dangl3,
> + dlat1, dlon1, dlat2, dlon2, ier )
> IF ( ier .ne. 0 ) iret = -5
>C
> RETURN
> END
>