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.
------- Forwarded Message >To: address@hidden, >To: address@hidden (Justin Sharp), >To: address@hidden (Harry Edmon) >From: David Ovens <address@hidden> >Subject: I found and fixed a bug in GEMPAK imar2gm.f >Organization: UCAR/Unidata >Keywords: 200012281724.eBSHOMo09342 To: address@hidden Justin Sharp, U.W. (address@hidden) Harry Edmon, U.W. (address@hidden) 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 LATLON=$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. I am also considering an optional projection "PROJ = sat2mer" that would remap a satellite file of any projection to a Mercator projection inside GEMPAK. I'll let you know if I ever get around to doing this.... 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* Chiz/Unidata 7/00 Added swap test initialization and * C* im_calib call if calibration exists * C* S. Danz/AWC 10/00 Fix RECT setting center longitude * C* D. Ovens/UW 12/00 Fix MERC,LAMB,PS Earth Radius mismatches* 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 imbswp = 0 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. 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 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' ) 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* Set pixel/grid spacing and earth radius. C dx = ianav ( 5 ) * iadir ( 13 ) re = ianav ( 8 ) if (re .le. 6200000.) re = RADIUS C C* Compute lat/lon of the lower-left corner. C dxp = 1. - rxp dyp = 1. - ryp alpha = 1. + SIN ( ABS ( phi0r ) ) rm = dx * SQRT ( dxp * dxp + dyp * dyp ) / alpha c 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 = dx * SQRT ( dxp * dxp + dyp * dyp ) / alpha c 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* Set pixel/grid spacing and earth radius and eccentricity. C dx = ianav ( 5 ) * iadir ( 13 ) 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 c 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* Set pixel/grid spacing and earth radius and eccentricity. C dx = ianav ( 5 ) * iadir ( 13 ) re = ianav ( 8 ) if (re .le. 6200000.) re = RADIUS 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 ) c 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* *** Radar *** C* ************* C ELSE IF ( navtyp .eq. 'RADR' ) THEN C C* Get radar site info. Assume radar site is center of image. C rrow = ianav ( 2 ) rcol = ianav ( 3 ) C C* Reference latitude. C dlat = stdlat C C* Reference longitude - convert to GEMPAK convention. C angdd = ianav ( 5 ) / 10000 angmm = ianav ( 5 ) / 100 - angdd * 100 angss = ianav ( 5 ) - angdd * 10000 - angmm * 100 dlon = angdd + angmm / 60. + angss / 3600. dlon = - dlon C C* Pixel resolution, north rotation (not used). C pxresx = ianav ( 6 ) pxresy = pxresx xrot = ianav ( 7 ) / 1000. * DTR C C* Calculate lat/lon pixel resolutions in radians C c rlatpx = ( RADIUS * DTR ) / pxresy c rlonpx = ( RADIUS * DTR ) / pxresx rlatpx = ( re * DTR ) / pxresy rlonpx = ( re * DTR ) / pxresx C C* Determine lower-left lat/lon. C yline = imnlin + iullne - 1 xelem = iulele C dyline = rrow - yline dxelem = rcol - xelem dlat1 = dlat + dyline / rlatpx dlon1 = dlon - dxelem / ( rlonpx * COS ( dlat1 * DTR ) ) C C* Determine upper-right lat/lon coords. C yline = iullne xelem = imnpix + iulele - 1 C dyline = rrow - yline dxelem = rcol - xelem dlat2 = dlat + dyline / rlatpx dlon2 = dlon - dxelem / ( rlonpx * COS ( dlat2 * DTR) ) C C* Set the map projection C dangl1 = dlat dangl2 = dlon dangl3 = 0. proj = 'STR' 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 ========================= imar2gm.f ==================== -- 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 ------- End of Forwarded Message