[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
20040224: Hyperslab I/O of Varying Size: use of IMAP argument
- Subject: 20040224: Hyperslab I/O of Varying Size: use of IMAP argument
- Date: Tue, 24 Feb 2004 09:11:53 -0700
Ben,
>Date: Tue, 24 Feb 2004 09:47:14 -0600
>From: Ben Howell <address@hidden>
>Organization: University of Wisconsin, Madison/SSEC
>To: Steve Emmerson <address@hidden>
>Subject: Re: 20040220: Hyperslab I/O of Varying Size: use of IMAP argument
> Keywords: 200402181846.i1IIkk2N001457
The above message contained the following:
> Thank you very much for your suggestions.
> I extracted the lines from your message and put them into my code.
> The latitudes still have only half, i.e. 125 of 250, values followed by
> 125 zeros.
The code doesn't look obviously wrong.
Are the values being put into the "slat" vector correctly? What is the
value of the "nlines" variable?
> The longitudes however were written correctly (600 values).
> I have attached my FORTRAN code which contains, toward the end, your
> suggested
> instructions. If you still have some patience (and time) let me know if
> there's any hope.
>
> Thanks!
> :-)
> ------------------------------------------------------------------------------------------------------------------------
...
> program MAIN
>
> implicit none
>
> include 'netcdf.inc'
>
> character*31 name
> character*80 stamp, ifile, logfile, msg
> character*256 cbuf
> integer*4 lentrim, IDCDF, IRET, nlines, nelems, nbands
> integer*4 imageDate, imageTime, i, j, k, n
> integer*4 maxbands, maxlines, maxelems, start(3), count(3)
>
>
> c ... 23Feb04:
> parameter(maxbands=1,maxlines=500, maxelems=1200)
> ccccc parameter(maxbands=1,maxlines=250, maxelems=600)
>
> real*4 xdata(maxelems, maxlines, maxbands)
>
> c ... Need 2 arrays of brightness temps, one for each output configuration:
>
> real*4 btemp1(maxelems, maxlines)
> real*4 btemp2(maxlines, maxelems)
>
> real*4 dlat(maxelems,maxlines), dlon(maxelems,maxlines)
> real*4 slat(maxlines), slon(maxelems)
> integer*4 idlinesDim, idelemsDim, idbandsDim
> integer*4 idimageDate, idimageTime, idxdata, idlat, idlon
> integer*4 iyear, idayoy
> integer*4 nyrday, kyear, kmon, kalday, log
> character*3 mon
> real*4 y, ymin, yavg, ymax, hours, hoursit, fmnts
> logical found
> integer*4 ihour, mmss, imnts, isecs
>
> c ... for output file:
> integer*4 idims(2), ndims, latdim, londim
> integer*4 idbt, idbch, idrch, idtl, idtlu, idcldy
> integer*4 iddate, idtime, idate, itime
> character*80 ofile, comnts
>
> INTEGER*4 IYEREF,IMOREF,IDAREF,IHOREF,IMIREF,ISEREF
> C ... DABTIME RETURNS ELAPSED TIME IN DAYS (DOUBLE PRECISION)
> C ... EDBASE = ELAPSED TIME OF OUR BASE TIME SINCE "EARLY REF."
> C ... EDNOW = ELAPSE TIME (DAYS) FROM "EARLY REF." UNTIL SPECIFIED DATE &
> TIME.
> REAL*8 DABTIM, EDREF, EDNOW, DVALU
>
> REAL*4 ABSENT
> INTEGER*4 JDAY, NTHDAY, IOPTDATE, IOPTOUT, NLATS, NLONS
>
> c ... 17Feb04:
> REAL*4 SIGN
>
> C ... 23FEB04:
> INTEGER*4 STRIDE(2), IMAP(2)
> DATA STRIDE/1,1/
>
> DATA ABSENT/-9999.0/
>
> C ... REFERENCE DATE AND TIME:
> DATA IYEREF/1970/, IMOREF/1/, IDAREF/1/, IHOREF/0/, IMIREF/0/,
> *ISEREF/0/
>
> c ---------------------------------------------------------------------
>
> data stamp/
> *'$Id: sat_globalbt.f,v 1.2(modified) 2004/02/23 $'/
> c ---------------------------------------------------------------------
>
> data start/1,1,1/, log/13/, ndims/2/
>
> c ***
>
> cbuf=stamp
> n=lentrim(cbuf)
> write(*,*)' '
> write(*,*)stamp(1:n)
> write(*,*)' '
>
> write(*,*)'Log (output, text) file:'
> read(*,'(a)')logfile
>
> write(*,*)'Input (MCIDAS 1-band NetCDF)file:'
> read(*,'(a)')ifile
> cbuf=ifile
> n=lentrim(cbuf)
> inquire(file=ifile, exist=found)
> if(found)then
> write(msg,*)'Input file:', cbuf(1:n)
> call logmsg(log,logfile,msg)
> else
> msg='Input file not found:'
> call logmsg(log, logfile, msg)
> write(msg,*)'File:',cbuf(1:n)
> call logmsg(log, logfile, msg)
> stop 'Input file?'
> endif
>
> write(*,*)'Output (NetCDF) file:'
> read(*,'(a)')ofile
>
> write(*,*)'Choose format for date:'
> write(*,*)' 1. YYYYDDD'
> write(*,*)' 2. YYYYMMDD'
> read(*,*)ioptdate
>
> write(*,*)'Choose configuration of output:'
> write(*,*)' 1. row-wise (latitude, longitude)'
> write(*,*)' 2. column-wise (longitude, latitude)'
> read(*,*)ioptout
>
> write(*,*)'Comments(up to 80 characters):'
> read(*,'(a)')comnts
>
> c Make errors non-fatal to allow us to do error handling.
> c The default is all errors are fatal = call ncpopt(NCVERBOS+NCFATAL)
> c
> call NCPOPT(0)
>
> c Open the netCDF file as read-only NCNOWRIT; get its "ID" IDCDF:
> c
> cbuf=ifile
> n=lentrim(cbuf)
> ifile(n+1:n+1)=char(0)
> IDCDF = NCOPN(ifile, NCNOWRIT, IRET)
> if (IRET.ne.0)stop 'Error opening file'
>
> c ... get dimension IDs from their names:
> idlinesDim = NCDID(IDCDF, 'lines', IRET)
> if(IRET.ne.0)stop 'Error, lines dimension'
>
> idelemsDim = NCDID(IDCDF, 'elems', IRET)
> if(IRET.ne.0)stop 'Error, elems dimension'
>
> idbandsDim = NCDID(IDCDF, 'bands', IRET)
> if(IRET.ne.0)stop 'Error, bands dimension'
>
> c ... Get the size of each dimension
> call NCDINQ(IDCDF, idlinesDim, name, nlines, IRET)
> if ( IRET.NE.0 ) stop 'Error getting size of lines dimension'
> cbuf=name
> n=lentrim(cbuf)
> write(msg,*)'dimension name:', name(1:n)
> call logmsg(log, logfile, msg)
> write(msg,*)'dimension size:', nlines
> call logmsg(log, logfile, msg)
>
> call NCDINQ(IDCDF, idelemsDim, name, nelems, IRET)
> if ( IRET.NE.0 ) stop 'Error getting sizse of elems dimension'
> cbuf=name
> n=lentrim(cbuf)
> write(msg,*)'dimension name:', name(1:n)
> call logmsg(log, logfile, msg)
> write(msg,*)'dimension size:', nelems
> call logmsg(log, logfile, msg)
>
> call NCDINQ(IDCDF, idbandsDim, name, nbands, IRET)
> if ( IRET.NE.0 ) stop 'Error getting size of bands dimension'
> cbuf=name
> n=lentrim(cbuf)
> write(msg,*)'dimension name:', name(1:n)
> call logmsg(log, logfile, msg)
> write(msg,*)'dimension size:', nbands
> call logmsg(log, logfile, msg)
>
> c ... Check size of "band" dimension:
> if(nbands.ne.1)then
> write(msg,*)'nbands:', nbands
> call logmsg(log, logfile, msg)
> msg='Wrong band dimension. Expect 1'
> call logmsg(log, logfile, msg)
> stop 'Bands?'
> endif
>
> c ... get variable IDs from their names
> idimageDate = NCVID(IDCDF, 'imageDate', IRET)
> if(IRET.ne.0)stop 'Error, imageDate ID'
>
> idimageTime = NCVID(IDCDF, 'imageTime', IRET)
> if(IRET.ne.0)stop 'Error, imageTime ID'
>
> idlat = NCVID(IDCDF, 'latitude', IRET)
> if(IRET.ne.0)stop 'Error, latitude ID'
>
> idlon = NCVID(IDCDF, 'longitude', IRET)
> if(IRET.ne.0)stop 'Error, longitude ID'
>
> idxdata = NCVID(IDCDF, 'data', IRET)
> if(IRET.ne.0)stop 'Error, data ID'
>
> c ... Set "count" vector
> count(1)=nelems
> count(2)=nlines
> count(3)=nbands
>
> c Get data (point values first, then time series)
> call NCVGT(IDCDF, idimageDate, 1, 1, imageDate, IRET)
> call NCVGT(IDCDF, idimageTime, 1, 1, imageTime, IRET)
>
> write(*,*)'imageDate, imageTime:',imageDate, imageTime
> write(msg,*)'imageDate, imageTime:',imageDate, imageTime
> call logmsg( log, logfile, msg)
>
> iyear=imageDate/1000
> idayoy=imageDate-1000*iyear
> iyear=iyear+1900
> write(*,*)'year, day-of-year:', iyear, idayoy
> write(msg,*)'year, day-of-year:', iyear, idayoy
> call logmsg( log, logfile, msg)
>
> c ... get calendar month and day
> nyrday=idayoy + 1000*iyear
> call calday(nyrday, kyear, kmon, kalday, mon)
> write(*,*)'Calendar year, month, day:', kyear, kmon, kalday
> write(msg,*)'Calendar year, month, day:', kyear, kmon, kalday
> call logmsg(log, logfile, msg)
>
> if(ioptdate.eq.1)then
> c ... Save date in YYYYDDD format:
> jday=nthday(kyear,kmon,kalday)
> idate=kyear*1000+jday
> else
> c ... Save date in YYYYMMDD format:
> idate=kyear*10000 + kmon*100 + kalday
> endif
>
>
>
> c ... Get hours, minutes and seconds (image time):
> ihour=imageTime/10000
> mmss=imageTime-10000*ihour
> imnts=mmss/100
> isecs=mmss-100*imnts
> write(msg,*)'Image time, hour, minute, second:',
> * ihour, imnts, isecs
> call logmsg(log,logfile,msg)
>
> hoursit=float(ihour)+float(imnts)/float(60)
> * +float(isecs)/float(3600)
> write(*,*)'Image time, decimal hours:', hoursit
> write(msg,*)'Image time, decimal hours:', hoursit
> call logmsg(log,logfile,msg)
>
> c ... Get "dlat" ...
> call NCVGT(IDCDF, idlat, start, count, dlat, IRET)
> if(IRET.ne.0)stop 'Error getting dlat'
>
> c ... Get "dlon" ...
> call NCVGT(IDCDF, idlon, start, count, dlon, IRET)
> if(IRET.ne.0)stop 'Error getting dlon'
>
> c ... Store unique latitudes:
> do j=1,nlines
> slat(j)=dlat(1,j)
> enddo
> nlats=nlines
>
> c ... Correction of last latitude
> if(slat(nlines) .lt. -90.0 .or. slat(nlines) .gt. 90.0)then
> write(msg,*)'slat(nlines), before correction:', slat(nlines)
> call logmsg(log, logfile, msg)
> slat(nlines)= -90.0
> write(msg,*)'lat(nlines), after correction:', slat(nlines)
> call logmsg(log, logfile, msg)
> endif
>
> c ... Store unique longitudes:
> do i=1,nelems
>
> c ... 17Feb04, input file "GOES1993.nc" appears to have longitudes
> c as degrees east, so no conversion is necessary.
> c To be safe, check first longitude
> if(i.eq.1)then
> if(dlon(1,1).lt.0.0)then
> sign=1.0
> else
> sign= -1.0
> endif
> endif
> slon(i)=dlon(i,1)*sign
>
> enddo
> nlons=nelems
>
> c ... Ingest "data" variable into "xdata" ...
> call NCVGT(IDCDF, idxdata, start, count, xdata, IRET)
> if(IRET.ne.0)stop 'Error getting xdata'
>
> c Close the netCDF file
> call NCCLOS(IDCDF, IRET)
> if(IRET.ne.0)stop'Error closing file'
>
> k=1
> do i=1,nelems
> do j=1,nlines
>
> c ... 18Feb04, check for "missing" data:
> if(nint(xdata(i,j,k)) .eq. 330)then
> xdata(i,j,k)=-9999.0
> endif
>
> if(ioptout.eq.1)then
> btemp1(i,j)=xdata(i,j,k)
> else
> btemp2(j,i)=xdata(i,j,k)
> endif
>
> enddo
> enddo
> write(*,*)'Stored brightness temp.'
>
> c ... Get min, mean, and max for each latitiude(line)
> msg='Brightness temperature, lat, min, mean, max:'
> call logmsg(log, logfile, msg)
> write(*,*)msg
> do j=1,nlines,50
> ymin=1.0e+10
> ymax=-ymin
> yavg=0.
> do i=1,nelems
> if(ioptout.eq.1)then
> y=btemp1(i,j)
> else
> y=btemp2(j,i)
> endif
> yavg=yavg+y
> if(y.lt.ymin)ymin=y
> if(y.gt.ymax)ymax=y
> enddo
> yavg=yavg/float(nelems)
> write(*,'(3x, f5.1, 3f7.1)')
> * slat(j),ymin,yavg,ymax
> write(msg,'(3x, f5.1, 3f7.1)')
> * slat(j),ymin,yavg,ymax
> call logmsg(log, logfile, msg)
> enddo
>
> CCCCCCCCCCC END OF INPUT SECTION, START OUTPUT SECTION CCCCCCCCC
>
> c ... Create output (NetCDF) file ...
> IDCDF=nccre(ofile,NCCLOB,IRET)
> if(IRET.ne.0)then
> cbuf=ofile
> n=lentrim(cbuf)
> write(*,*)'Unable to open file:', cbuf(1:n)
> write(msg,*)'Unable to open file:', cbuf(1:n)
> call logmsg(log, logfile, msg)
> endif
>
> C ... DEFINE DIMENSIONS:
>
> LATDIM=NCDDEF(IDCDF, 'latitude', nlats, IRET)
> LONDIM=NCDDEF(IDCDF, 'longitude', nlons, IRET)
>
> C ... DEFINE VARIABLE IDS:
> IDDATE=NCVDEF(IDCDF, 'date', NCLONG, 0, 0, IRET)
> IDTIME=NCVDEF(IDCDF, 'time', NCLONG, 0, 0, IRET)
>
> IDLAT=NCVDEF(IDCDF, 'latitude', NCFLOAT, 1, LATDIM, IRET)
> IDLON=NCVDEF(IDCDF, 'longitude', NCFLOAT, 1, LONDIM, IRET)
>
> if(ioptout.eq.1)then
> IDIMS(1)=LONDIM
> IDIMS(2)=LATDIM
> else
> IDIMS(1)=LATDIM
> IDIMS(2)=LONDIM
> endif
>
> IDBT=NCVDEF(IDCDF, 'brightness_temp',NCFLOAT,NDIMS,IDIMS,IRET)
> IF(IRET.NE.0)STOP 'IDBT?'
>
>
> C ... ASSIGN ATTRIBUTES TO VARIABLES
>
> c ... Date:
>
> if(ioptdate.eq.1)then
> CALL NCAPTC(IDCDF, IDDATE, 'long_name', NCCHAR, 29,
> *' Image date, in YYYYDDD format', IRET)
> else
> CALL NCAPTC(IDCDF, IDDATE, 'long_name', NCCHAR, 30,
> *' Image date, in YYYYMMDD format', IRET)
> endif
>
> c ... Time:
> CALL NCAPTC(IDCDF, IDTIME, 'long_name', NCCHAR, 15,
> *'Image time, UTC', IRET)
> CALL NCAPTC(IDCDF,IDTIME,'units',NCCHAR,31,
> *'seconds since 1970-1-1 00:00:00',IRET)
> CALL NCAPT(IDCDF, IDTIME, 'image_time_HHMMSS_UTC', NCLONG, 1,
> *imageTime, IRET)
>
> c ... Latitude:
> CALL NCAPTC(IDCDF, IDLAT, 'long_name', NCCHAR, 8,
> *'Latitude', IRET)
> CALL NCAPTC(IDCDF, IDLAT, 'units', NCCHAR, 13,
> *'degrees_north', IRET)
> CALL NCAPTC(IDCDF, IDLAT, 'precision', NCCHAR, 4,
> *'1E-1', IRET)
>
> c ... Longitude:
> CALL NCAPTC(IDCDF, IDLON, 'long_name', NCCHAR, 9,
> *'Longitude', IRET)
> CALL NCAPTC(IDCDF, IDLON, 'units', NCCHAR, 12,
> *'degrees_east', IRET)
> CALL NCAPTC(IDCDF, IDLON, 'precision', NCCHAR, 4,
> *'1E-1', IRET)
>
> c ... Brightness temperature
> CALL NCAPTC(IDCDF, IDBT, 'long_name', NCCHAR, 32,
> *'Radiative brightness temperature', IRET)
> CALL NCAPTC(IDCDF, IDBT, 'units', NCCHAR, 14,
> *'degrees_Kelvin', IRET)
> CALL NCAPTC(IDCDF, IDBT, 'precision', NCCHAR, 4,
> *'1E-1', IRET)
> CALL NCAPT(IDCDF, IDBT, 'missing_value', NCFLOAT, 1,
> *ABSENT, IRET)
>
> C ... ADD USER COMMENTS TO "Global Attributes" ...
> CBUF=COMNTS
> N=LENTRIM(CBUF)
> IF(N.GT.79)N=79
> IF(N.LT.1)THEN
> COMNTS='No comments'
> N=11
> ENDIF
> CALL NCAPTC(IDCDF,NCGLOBAL,'Comments',NCCHAR,N,COMNTS,IRET)
>
> C ... ADD SOFTWARE VERSION TO "Global Attributes" ...
> CBUF=STAMP
> N=LENTRIM(CBUF)
> IF(N.GT.79)N=79
> stamp(n+1:n+1)=char(0)
> CALL NCAPTC(IDCDF,NCGLOBAL,'SoftwareVersion',NCCHAR,N,STAMP,
> * IRET)
>
> C ... LEAVE DEFINE MODE ...
> CALL NCENDF(IDCDF, IRET)
>
> C ... COMPUTE ELAPSED DAYS FROM EARLY REFERENCE TO 1 JAN 1970
> EDREF=DABTIM(IYEREF,IMOREF,IDAREF,IHOREF,IMIREF,ISEREF)
>
> C ... GET BASE TIME (SECONDS ELAPSED SINCE 1 JAN 1970, 0000Z)
> EDNOW=DABTIM(KYEAR,KMON,KALDAY,IHOUR,IMNTS,ISECS)
> DVALU=(EDNOW-EDREF)*DFLOAT(86400)
> ITIME=DNINT(DVALU)
>
>
> write(*,*)'Writing data to NetCDF file...'
>
> START(1)=1
> COUNT(1)=1
>
> C ... STORE DATE ...
> CALL NCVPT(IDCDF, IDDATE, START, COUNT, IDATE, IRET)
> IF(IRET.EQ.0)THEN
> WRITE(*,*)'Wrote date:', idate
> ELSE
> STOP 'ERROR WRITING DATE'
> ENDIF
>
>
> C ... STORE TIME ...
> CALL NCVPT(IDCDF, IDTIME, START, COUNT, ITIME, IRET)
> IF(IRET.EQ.0)THEN
> WRITE(*,*)'Wrote image time:', itime
> ELSE
> STOP 'ERROR WRITING IMAGE TIME'
> ENDIF
>
> C ... STORE LATITUDES:
> C ... 24FEB04, as suggested by Steve Emmerson(Unidata):
>
>
>
> C ... OUTPUT LATITUDES:
> START(1)=1
> COUNT(1)=NLATS
> c IMAP(1)=1
>
> c ... temporary:
> write(*,*)'latitudes: start(1), count(1):', start(1), count(1)
> pause
>
> c IRET=NF_PUT_VARM_REAL(IDCDF,IDLAT,START,COUNT, ! no
> c *STRIDE,IMAP,SLAT) ! no
> IRET=NF_PUT_VARA_REAL(IDCDF,IDLAT,START,COUNT, ! yes
> *SLAT) ! yes
>
> C ... OUTPUT LONGITUDES:
> START(1)=1
> COUNT(1)=NLONS
> c IMAP(1)=1 ! no
>
> c ... temporary:
> write(*,*)'longitudes: start(1), count(1):', start(1), count(1)
> pause
>
> c IRET=NF_PUT_VARM_REAL(IDCDF,IDLON,START,COUNT, ! no
> c *STRIDE,IMAP,SLON) ! no
> IRET=NF_PUT_VARA_REAL(IDCDF,IDLON,START,COUNT, ! yes
> *SLON) ! yes
>
> CCCCC CALL NCVPT(IDCDF, IDLAT, START, COUNT, SLAT, IRET)
>
> IF(IRET.EQ.0)THEN
> WRITE(*,*)'Wrote latitudes'
> ELSE
> STOP 'ERROR WRITING LATITUDES'
> ENDIF
>
> START(1)=1
> COUNT(1)=NLONS
> IMAP(1)=1
>
> C ... STORE LONGITUDES:
>
> C ... 24FEB04, as suggested by Steve Emmerson(Unidata):
> IRET=NF_PUT_VARA_REAL(IDCDF,IDLON,START,COUNT,SLON)
> ! yes
> CCCCC CALL NCVPT(IDCDF, IDLON, START, COUNT, SLON, IRET)
>
> IF(IRET.EQ.0)THEN
> WRITE(*,*)'Wrote longitudes'
> ELSE
> STOP 'ERROR WRITING LONGITUDES'
> ENDIF
>
> C ... 11FEB04:
> IF(IOPTOUT.EQ.1)THEN
> START(1)=1
> COUNT(1)=NLONS
> START(2)=1
> COUNT(2)=NLATS
> IMAP(1)=1
>
> C ... 24FEB04, as suggested by Steve Emmerson(Unidata):
> IMAP(2)=MAXELEMS ! element count between 2nd index increments
>
> C ... 23FEB04:
> IRET=NF_PUT_VARM_REAL(IDCDF,IDBT,START,COUNT,
> * STRIDE,IMAP,BTEMP1)
> CCCCC CALL NCVPT(IDCDF, IDBT, START, COUNT, BTEMP1, IRET)
>
> ELSE
> START(1)=1
> COUNT(1)=NLATS
> START(2)=1
> COUNT(2)=NLONS
> IMAP(1)=1
>
> C ... 24FEB04, as suggested by Steve Emmerson(Unidata):
> IMAP(2)=MAXLINES ! element count between 2nd index increments
>
> C ... 23FEB04:
> IRET=NF_PUT_VARM_REAL(IDCDF,IDBT,START,COUNT,
> * STRIDE,IMAP,BTEMP2)
> CCCCC CALL NCVPT(IDCDF, IDBT, START, COUNT, BTEMP2, IRET)
>
> ENDIF
> IF(IRET.EQ.0)THEN
> WRITE(*,*)'Wrote brightness temps'
> ELSE
> STOP 'ERROR WRITING BR.TEMPS'
> ENDIF
>
>
> C ... CLOSE NETCDF FILE ..
> CALL NCCLOS(IDCDF,IRET)
> IF(IRET.NE.0)THEN
> WRITE(*,*)'Bad status from NCCLOS:', IRET
> STOP 'NCCLOS error'
> ENDIF
>
> msg='*** End sat_globalbt ***'
> call logmsg(log, logfile, msg)
>
> stop '*** End sat_globalbt ***'
> end
Regards,
Steve Emmerson