This archive contains answers to questions sent to Unidata support through mid-2025. Note that the archive is no longer being updated. We provide the archive for reference; many of the answers presented here remain technically correct, even if somewhat outdated. For the most up-to-date information on the use of NSF Unidata software and data services, please consult the Software Documentation first.
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 >