http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_f.f ---------------------------------------------------------------------- diff --git a/climatology/clim/gaussInterp_f.f b/climatology/clim/gaussInterp_f.f new file mode 100644 index 0000000..9366870 --- /dev/null +++ b/climatology/clim/gaussInterp_f.f @@ -0,0 +1,219 @@ +c +c gaussInterp routine -- Gaussian weighted smoothing in lat, lon, and time +c +c Based on Ed Armstrong's routines. Designed to be called from python using f2py. +c +c +c Gaussian weighting = exp( vfactor * (((x - x0)/sx)^2 + ((y - y0)/sy)^2 + ((t - t0)/st)^2 )) +c +c where deltas are distances in lat, lon and time and sx, sy, st are one e-folding sigmas. +c +c Cutoffs for neighbors allowed in the interpolation are set by distance in lat/lon (see wlat/wlon); +c for time all epochs are included. + + + subroutine gaussInterp_f(var, mask, + & time, lat, lon, + & outlat, outlon, + & wlat, wlon, + & slat, slon, stime, + & vfactor, + & missingValue, + & verbose, + & vinterp, vweight, status, + & ntime, nlat, nlon, noutlat, noutlon) + + implicit none + integer*4 ntime, nlat, nlon +c ! prepared 3D array of variable data over time, lon, & lat + real*4 var(ntime, nlat, nlon) +c ! variable to be interpolated, over ntime epochs + integer*1 mask(ntime, nlat, nlon) +c ! pixel quality mask + integer*4 time(ntime) +c ! time epochs to gaussian-interpolate over + real*4 lat(nlat) + ! latitude corrdinate vector + real*4 lon(nlon) + ! longitude corrdinate vector +cf2py intent(in) var, mask, time, lat, lon !! annotations for f2py processor + integer*4 noutlat, noutlon + real*4 outlat(noutlat), outlon(noutlon) +c ! lat/lon grid to interpolate to +cf2py intent(in) outlat, outlon + real*4 wlat, wlon +c ! window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees) +c ! if an integer, then it is the number of neighbors on EACH side +cf2py intent(in) wlat, wlon + real*4 slat, slon, stime +c ! sigma for gaussian downweighting with distance in lat, lon, & time +cf2py intent(in) slat, slon, stime + real*4 vfactor +c ! factor in front of gaussian expression + real*4 missingValue +c ! value to mark missing values in interp result + integer*4 verbose +c ! integer to set verbosity level +cf2py intent(in) vfactor, missingValue + real*4 vinterp(noutlat, noutlon) +c ! interpolated variable using gaussians, missing values not counted + real*4 vweight(noutlat, noutlon) +c ! weight of values interpolated (can be zero), if so should become missing value + integer*4 status +c ! negative status indicates error +cf2py intent(out) vinterp, vweight, status + + integer*4 clamp + integer*4 imin, imax, jmin, jmax, itmp + integer*4 iin, jin, kin + integer*4 iMidTime + integer*4 i, j, iwlat, iwlon + + real*4 wlat2, wlon2, lat1, lon1, dlat, dlon, fac + real*4 varmin, varmax, val + real*8 midTime + logical*1 nLatNeighbors, nLonNeighbors + + write(6, *) 'Echoing inputs ...' + write(6, *) 'ntime, nlat, nlon:', ntime, nlat, nlon + write(6, *) 'noutlat, noutlon:', noutlat, noutlon + write(6, *) 'wlat, wlon:', wlat, wlon + write(6, *) 'slat, slon, stime:', slat, slon, stime + write(6, *) 'vfactor:', vfactor + write(6, *) 'missingValue:', missingValue + + status = 0 + if (int(wlat) .eq. wlat) then + iwlat = int(wlat) + nLatNeighbors = .TRUE. + write(6, *) 'Using', iwlat, + & 'neighbors on each side in the lat direction.' + else + nLatNeighbors = .FALSE. + end if + if (int(wlon) .eq. wlon) then + iwlon = int(wlon) + nLonNeighbors = .TRUE. + write(6, *) 'Using', iwlon, + & 'neighbors on each side in the lon direction.' + else + nLonNeighbors = .FALSE. + end if + + wlat2 = wlat / 2. + wlon2 = wlon / 2. + lat1 = lat(1) + lon1 = lon(1) + dlat = lat(2) - lat(1) + dlon = lon(2) - lon(1) + iMidTime = int(ntime/2 + 0.5) + midTime = time(iMidTime) + + if (verbose .gt. 3) then + write(6, *) 'time:', time + write(6, *) 'lat:', lat + write(6, *) 'lon:', lon + write(6, *) 'outlat:', outlat + write(6, *) 'outlon:', outlon +c write(6, *) 'mask(iMidTime):', mask(iMidTime,:,:) + write(6, *) 'var(iMidTime):', var(iMidTime,:,:) + end if + + do i = 1, noutlat + if (verbose .gt. 1) write(6, *) outlat(i) + do j = 1, noutlon + vinterp(i,j) = 0.0 + vweight(i,j) = 0.0 + if (verbose .gt. 3) then + write(6, *) '(i,j) = ', i, j + write(6, *) '(outlat,outlon) = ', outlat(i), outlon(j) + end if + + if (nLatNeighbors) then + imin = clamp(i - iwlat, 1, nlat) + imax = clamp(i + iwlat, 1, nlat) + else + imin = clamp(int((outlat(i) - wlat2 - lat1)/dlat + 0.5), + & 1, nlat) + imax = clamp(int((outlat(i) + wlat2 - lat1)/dlat + 0.5), + & 1, nlat) + end if + if (imin .gt. imax) then +c ! input latitudes descending, so swap + itmp = imin + imin = imax + imax = itmp + end if + + if (nLonNeighbors) then + jmin = clamp(j - iwlon, 1, nlon) + jmax = clamp(j + iwlon, 1, nlon) + else + jmin = clamp(int((outlon(j) - wlon2 - lon1)/dlon + 0.5), + & 1, nlon) + jmax = clamp(int((outlon(j) + wlon2 - lon1)/dlon + 0.5), + & 1, nlon) + end if + + if (verbose .gt. 2) then + write(6, *) '(imin,imax,minlat,maxlat) = ', + & imin, imax, lat(imin), lat(imax) + write(6, *) '(jmin,jmax,minlon,maxlon) = ', + & jmin, jmax, lon(jmin), lon(jmax) + end if + + if (verbose .gt. 4 .and. i .eq. noutlat/2 + & .and. j .eq. noutlon/2) then + write(6, *) 'Echoing factors ...' + end if + + do kin = 1, ntime + varmin = 0. + varmax = 0. + do iin = imin, imax + do jin = jmin, jmax + if (mask(kin,iin,jin) .eq. 0) then + fac = exp( vfactor * + & (((outlat(i) - lat(iin))/slat)**2 + & + ((outlon(j) - lon(jin))/slon)**2 + & + ((midTime - time(kin))/stime)**2)) + + val = var(kin,iin,jin) + if (verbose .gt. 4 .and. i .eq. noutlat/2 + & .and. j .eq. noutlon/2) then + write(6, *) kin, iin, jin, time(kin), + & lat(iin), lon(jin), val, fac, val*fac + end if + + vinterp(i,j) = vinterp(i,j) + val * fac + vweight(i,j) = vweight(i,j) + fac + if (verbose .gt. 3) then + varmin = min(varmin, val) + varmax = max(varmax, val) + end if + end if + end do + end do + if (verbose .gt. 3) then + write(6, *) '(itime, varmin, varmax) = ', + & kin, varmin, varmax + end if + end do + if (vweight(i,j).ne. 0.0) then + vinterp(i,j) = vinterp(i,j) / vweight(i,j) + else + vinterp(i,j) = missingValue + end if + end do + end do + return + end + + + integer*4 function clamp(i, n, m) + integer*4 i, n, m + clamp = i + if (i .lt. n) clamp = n + if (i .gt. m) clamp = m + end +
http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_f.mk ---------------------------------------------------------------------- diff --git a/climatology/clim/gaussInterp_f.mk b/climatology/clim/gaussInterp_f.mk new file mode 100644 index 0000000..548ce33 --- /dev/null +++ b/climatology/clim/gaussInterp_f.mk @@ -0,0 +1 @@ +f2py -c -m gaussInterp_f gaussInterp_f.f http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/gaussInterp_slow.py ---------------------------------------------------------------------- diff --git a/climatology/clim/gaussInterp_slow.py b/climatology/clim/gaussInterp_slow.py new file mode 100644 index 0000000..6a7bb8f --- /dev/null +++ b/climatology/clim/gaussInterp_slow.py @@ -0,0 +1,130 @@ +# +# gaussInterp_slow routine -- Gaussian weighted smoothing in lat, lon, and time +# +# Based on Ed Armstrong's routines. Pure python implementation. +# +# +# Gaussian weighting = exp( vfactor * (((x - x0)/sx)^2 + ((y - y0)/sy)^2 + ((t - t0)/st)^2 )) +# +# where deltas are distances in lat, lon and time and sx, sy, st are one e-folding sigmas. +# +# Cutoffs for neighbors allowed in the interpolation are set by distance in lat/lon (see dlat/dlon); +# for time all epochs are included. +# + +import sys +import numpy as np +from math import exp +from numba import jit, int32 + +VERBOSE = 0 + + +def gaussInterp_slow(var, # bundle of input arrays: masked variable, coordinates + varNames, # list of names in order: primary variable, coordinates in order lat, lon, time + outlat, outlon, # output lat/lon coordinate vectors + wlat, wlon, # window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees) + slat, slon, stime, # sigma for gaussian downweighting with distance in lat, lon (deg), & time (days) + vfactor=-0.6931, # factor in front of gaussian expression + missingValue=-9999., # value to mark missing values in interp result + verbose=VERBOSE, # integer to set verbosity level + optimization='python'): # Mode of optimization, using 'fortran' or 'cython' or 'python' + '''Gaussian interpolate in lat, lon, and time to a different lat/lon grid, and over a time window to the center time. +Bundle of arrays (var) contains a 3D masked variable and coordinate arrays for lat, lon, and time read from netdf/hdf files. +Returns the 2D interpolated variable (masked) and a status for failures. + ''' + v = var[varNames[0]][:] + vmask = np.ma.getmask(v)[:] + vtime = var[varNames[1]][:] + lat = var[varNames[2]][:] + lon = var[varNames[3]][:] + + vinterp, vweight, status = \ + gaussInterp_(v, vmask, vtime, lat, lon, + outlat, outlon, wlat, wlon, slat, slon, stime, vfactor, missingValue) + + vinterp = np.ma.masked_where(vweight == 0.0, vinterp) + return (vinterp, vweight, status) + +#@jit(nopython=False) +def gaussInterp_(var, # variable & mask arrays with dimensions of time, lon, lat + vmask, + vtime, # coordinate vectors for inputs + lat, + lon, + outlat, # coordinate vectors for grid to interpolate to + outlon, + wlat, wlon, # window of lat/lon neighbors to gaussian weight, expressed in delta lat (degrees) + slat, slon, stime, # sigma for gaussian downweighting with distance in lat, lon (deg), & time (days) + vfactor, # factor in front of gaussian expression + missingValue): # value to mark missing values in interp result + '''Gaussian interpolate in lat, lon, and time to a different lat/lon grid, and over a time window to the center time. +Returns the 2D interpolated variable (masked), the weight array, and a status for failures. + ''' + vinterp = np.zeros( (outlat.shape[0], outlon.shape[0]), dtype=var.dtype ) # interpolated variable, missing values not counted + vweight = np.zeros( (outlat.shape[0], outlon.shape[0]), dtype=var.dtype ) # weight of values interpolated (can be zero) + status = 0 # negative status indicates error + + ntime = vtime.shape[0] + nlat = lat.shape[0] + nlon = lon.shape[0] + + noutlat = outlat.shape[0] + noutlon = outlon.shape[0] + + midTime = vtime[int(ntime/2 + 0.5)] + wlat2 = wlat / 2. + wlon2 = wlon / 2. + lat0 = lat[0] + lon0 = lon[0] + dlat = lat[1] - lat[0] + dlon = lon[1] - lon[0] + + for i in xrange(noutlat): + print >>sys.stderr, outlat[i] + for j in xrange(noutlon): + if VERBOSE: print >>sys.stderr, '\n(i,j) = %d, %d' % (i, j) + if VERBOSE: print >>sys.stderr, '\n(outlat,outlon) = %f, %f' % (outlat[i], outlon[j]) + + imin = clamp(int((outlat[i] - wlat2 - lat0)/dlat + 0.5), 0, nlat-1) + imax = clamp(int((outlat[i] + wlat2 - lat0)/dlat + 0.5), 0, nlat-1) + if imin > imax: (imin, imax) = (imax, imin) # input latitudes could be descending + if VERBOSE: print >>sys.stderr, '(imin, imax) = %d, %d' % (imin, imax) + if VERBOSE: print >>sys.stderr, '(minlat, maxlat) = %f, %f' % (lat[imin], lat[imax]) + jmin = clamp(int((outlon[j] - wlon2 - lon0)/dlon + 0.5), 0, nlon-1) + jmax = clamp(int((outlon[j] + wlon2 - lon0)/dlon + 0.5), 0, nlon-1) + if VERBOSE: print >>sys.stderr, '(jmin, jmax) = %d, %d' % (jmin, jmax) + if VERBOSE: print >>sys.stderr, '(minlon, maxlon) = %f, %f' % (lon[jmin], lon[jmax]) +# stencil = np.zeros( (ntime, imax-imin+1, jmax-jmin+1) ) + + for kin in xrange(ntime): + for iin in xrange(imin, imax+1): + for jin in xrange(jmin, jmax+1): + if not vmask[kin,iin,jin]: + fac = exp( vfactor * + (((outlat[i] - lat[iin])/slat)**2 + + ((outlon[j] - lon[jin])/slon)**2 + + ((midTime - vtime[kin])/stime)**2)) +# stencil[kin, iin-imin, jin-jmin] = fac + val = var[kin, iin, jin] + if VERBOSE > 1: print >>sys.stderr, kin, iin, jin, vtime[kin], lat[iin], lon[jin], val, fac, val*fac + + vinterp[i,j] = vinterp[i,j] + val * fac + vweight[i,j] = vweight[i,j] + fac + + + if vweight[i,j] != 0.0: + vinterp[i,j] = vinterp[i,j] / vweight[i,j] +# if VERBOSE > 1: print >>sys.stderr, 'stencil:\n', stencil +# if VERBOSE: print >>sys.stderr, 'stencil max:\n', np.max(np.max(stencil)) +# if VERBOSE: print >>sys.stderr, 'stencil min:\n', np.min(np.min(stencil)) + else: + vinterp[i,j] = missingValue + + return (vinterp, vweight, status) + +#@jit( int32(int32,int32,int32), nopython=False) +def clamp(i, n, m): + if i < n: return n + if i > m: return m + return i http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/interp.f ---------------------------------------------------------------------- diff --git a/climatology/clim/interp.f b/climatology/clim/interp.f new file mode 100644 index 0000000..8977798 --- /dev/null +++ b/climatology/clim/interp.f @@ -0,0 +1,302 @@ +#include "interp.h" +c $Revision: 2.12 $ +c subprogram interp.f +c reads hdf file, extract data, call bisum(), calculates +c and writes out interpolated maps + + + subroutine interp( temp, sst_inter, sst_seq, sl, + & x, y, z, f, + & xx, yy, zz, vsum, ixmin, ixmax, + & iymin, iymax, izmin, izmax, ndatamax, + & x_length, y_length, imax, jmax, kmax ) + + integer x_length,y_length + integer dummy + +c-- dimension the arrays + real*4 temp(x_length,y_length) + real*4 sl(imax,jmax,kmax) + real*4 sst_inter(kmax), sst_seq(imax) + real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax) + real*4 xx(imax),yy(jmax),zz(kmax) + real*4 vsum(2,imax,jmax,kmax) + real*4 ULlon, ULlat, LRlon, LRlat + integer*4 ixmin(ndatamax),ixmax(ndatamax) + integer*4 iymin(ndatamax),iymax(ndatamax) + integer*4 izmin(ndatamax),izmax(ndatamax) + +c-- hdf data array must be hardwired, no dync alloc :( + byte in_data( XSIZE, YSIZE ) + + character*80 plabel + character*50 infile + character*30 outfile + character*150 datafile + character*30 basename + character*35 tmpfile + character*10 mapformat + character*15 outformat + character*3 cday + character*4 cyr + character*200 zcatcmd + + common /parms/ imin,jmin,kmin, + & xwin2,ywin2,zwin2, + & xh,yh,zh, + & dx,dy,dz, + & xlo,xhi,ylo,yhi,zlo, + & dxd,dyd + + common /fileio/ infile,outfile,istart, + & iend,mapformat,outformat + + common /hdfindex/ icolLeft,icolRight,irowUpper,irowLower + + common /basefile/ basename + +c-- open output file + open(unit=20,status='new',form=outformat,file=outfile, + & err=1000,iostat=ierr ) + plabel='sst interpolation' + +c-- initialize vsum to '0's + do k1=1,2 + do k2=1,kmax + do k3=1,jmax + do k4=1,imax + vsum(k1,k4,k3,k2)=0. + enddo + enddo + enddo + enddo + +c-- initialize arrays of long, lat and time intervals + do ii=imin,imax + xx(ii)=xlo+float(ii-1)*dx ! Long to interpolate to + enddo + do jj=jmin,jmax + yy(jj)=ylo+float(jj-1)*dy ! Lat to interpolate to + enddo + do kk=kmin,kmax + zz(kk)=zlo+float(kk-1)*dz ! Time to interpolate to + enddo + + do n=1,ndatamax + x(n)=0.0 + y(n)=0.0 + z(n)=0.0 + f(n)=0.0 + ixmin(n) = 0 + ixmax(n) = 0 + iymin(n) = 0 + iymax(n) = 0 + izmin(n) = 0 + izmax(n) = 0 + enddo + +c-- slope (cal) and offset to convert DN to SST + cal= 0.15 + offset = 3.0 + +c-- Open input file list + open(unit=15,status='old',file=infile,access='direct', + & recl=RECSIZE,form='formatted',err=1500,iostat=ierr) + +c-- read a infile one record at a time and process . . . + do k=istart,iend + ipts=0 + read( 15,'(a)',rec=k ) datafile +c print *, datafile + + icompress = 0 +c-- hack to get trailing 'Z' (if compressed) +c do islash = 150, 1, -1 +c if ( datafile(islash:islash) .eq. 'Z' ) then +c icompress = 1 +c goto 101 +c endif +c enddo +c101 continue + +c-- 12/27/99: C call for basename implemented +c-- variable basename returned via common statement + call getbase( datafile ) + print *,'\n file: ', basename + +c-- cyr and cday parsed from basename +c-- (e.g., parsed from '1988332h54dd-gdm.hdf') + cday=basename( 5:7 ) + cyr=basename( 1:4 ) + +c-- if unix compressed, zcat to a temporary filename + if( basename(22:22) .eq. 'Z' ) then + icompress = 1 + tmpfile = '/tmp/' // basename(1:20) + zcatcmd = 'zcat ' // datafile( 1:len(datafile)-1 ) + & // ' > ' // tmpfile + call system( zcatcmd ) + datafile = tmpfile + endif + +c-- convert iday character to integer + read(cday,'(i3)') iday + read(cyr,'(i4)') iyr +c write(6,*)'cday = ',cday +c write(6,*)'cyr = ',cyr + +c***> HDF call: d8gimg() + retn=d8gimg(datafile,in_data,x_length, + & y_length,dummy) + +c-- read hdf DN values and convert to SST + do i=icolLeft, icolRight + do j=irowUpper, irowLower + xlat=-89.75+dyd*float(j-1) + +c-- center output in Pacific Ocean +c if( i .lt. x_length ) ix = i-(x_length/2) +c if( i .le. (x_length/2) ) ix = i+(x_length/2) + +c-- center output in Atlantic (default) + ix = i + +c-- convert signed byte to float +c if ( in_data(i,j) .lt. 0 .and. abs(xlat) .lt. 70. ) then + if ( in_data(i,j).lt. 0 ) then + temp(ix,j)=float( in_data(i,j) ) + 256 + else + temp(ix,j)=float( in_data(i,j) ) + endif + + + +C MULTIPLY THE PATHFINDER DIGITAL NUMBER BY THE CALIBRATION NUMBER (0.15) +C AND ADD THE OFFSET (-3.0) TO GET DEGREES CELSIUS + + if ( temp(ix,j).gt.0 ) then + ipts=ipts+1 + f(ipts)=( cal*temp(ix,j) ) - offset + x(ipts) = ( dxd*float(ix-1) ) - 180. + dxd/2 + y(ipts) = ( 90. - (dyd*float(j-1)) ) - dyd/2 + z(ipts)=float(iday)+float(iyr-1985)*365. + if(iyr.gt.1988) z(ipts)=z(ipts)+1 + if(iyr.gt.1992) z(ipts)=z(ipts)+1 + if(iyr.gt.1996) z(ipts)=z(ipts)+1 + if(iyr.gt.2000) z(ipts)=z(ipts)+1 + if(iyr.gt.2004) z(ipts)=z(ipts)+1 + endif + enddo + enddo + + nfnew=ipts + print *, ' no of pts in file ',' = ', nfnew + +c-- calculate interpolation weights and vsum +c-- arrays passed directly... a common statement +c-- does not seem to work. + call binsum( nfnew, x, y, z, f, + & xx, yy, zz, vsum, + & ixmin, ixmax, iymin, + & iymax, izmin, izmax, + & imax, jmax, kmax, ndatamax ) + + if( icompress .eq. 1 ) call system( 'rm -f ' // tmpfile ) +c-- ..... read next hdf file from infile + enddo + +c-- all input files processed; calculate interpolated SST + do kk=1,kmax + do jj=1,jmax + do ii=1,imax + sl(ii,jj,kk)=0. + if (vsum(2,ii,jj,kk).gt.0) then + sl(ii,jj,kk)=vsum(1,ii,jj,kk)/vsum(2,ii,jj,kk) + endif + enddo + enddo + enddo + +c-- write output as map "interleaved" or map "sequential" +c-- "interleaved" is the original implementation +c-- both formats preceded by header + +c-- "interleaved" refers to each row in the output +c-- file representing a unique lon, lat position with columns +c-- of associated sst values (maps) +c-- i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst +c-- row two: lon2 lat1 sst1 sst2 sst3....lastsst + +c-- "sequential" refers to each map written as rows and +c-- columns of lat, lon with the array element representing the +c-- sst at that geo-position for that map. +c-- each map is written sequentially in the file + +c-- geo-positions of UL and LR corners + ULlon = -(180 - ( ((icolLeft-1) * dxd) + dxd/2 )) + ULlat = (90 - ( ((irowUpper-1) * dyd) + dyd/2 )) + LRlon = -(180 - ( (icolRight * dxd) - dxd/2 )) + LRlat = (90 - ( (irowLower * dyd) - dyd/2 )) + +c-- version number, "f" refers to fortran version + version = 'f2.9' + +c-- write the 3 record header + if( outformat .eq. 'formatted' ) then + write(20,'(a)'), plabel + write(20,*) imax,jmax,kmax + write(20,*) dx,dy,dz + write(20,*) ULlon,ULlat,LRlon,LRlat + elseif( outformat .eq. 'unformatted' ) then + write(20) imax,jmax,kmax + write(20) dx,dy,dz + write(20) ULlon,ULlat,LRlon,LRlat + endif + + if( mapformat .eq. 'interleave' ) then + print *, '\n map output is interleave' + do jj=1,jmax + do ii=1,imax + do kk=1,kmax + sst_inter(kk)=sl(ii,jj,kk) + enddo + if( outformat .eq. 'formatted' ) then + write(20,*) ii,jj,sst_inter + else + write(20) ( sst_inter(i), i=1,kmax ) + endif + enddo + enddo + + else if( mapformat .eq. 'sequential' ) then + print *, '\n map output is sequential' + do kk=1,kmax + do jj=1,jmax + do ii=1,imax + sst_seq(ii)=sl(ii,jj,kk) + enddo + if( outformat .eq. 'formatted' ) then + write(20,*) jj,kk,sst_seq + else + write(20) ( sst_seq(i), i=1,imax ) + endif + enddo + enddo + endif + + print *, '\ndone. . . ' + close( 15,status='keep',err=2000,iostat=ierr ) + close( 20,status='keep',err=2500,iostat=ierr ) + stop + + 1000 print *, 'Error opening output: ', outfile, 'error num: ', ierr + goto 102 + 1500 print *, 'Error opening input: ', infile, 'error num: ', ierr + goto 102 + 2000 print *, 'Error closing input: ', infile, 'error num: ', ierr + goto 102 + 2500 print *, 'Error closing output: ', outfile, 'error num: ', ierr + goto 102 + 102 continue + + end http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/jobClimatology2.py ---------------------------------------------------------------------- diff --git a/climatology/clim/jobClimatology2.py b/climatology/clim/jobClimatology2.py new file mode 100644 index 0000000..4b3933e --- /dev/null +++ b/climatology/clim/jobClimatology2.py @@ -0,0 +1,22 @@ +# +# jobClimatology2.py -- mrjob class for computing climatologies by gaussian interpolation +# + +from mrjob.job import MRJob +from climatology2 import climByAveragingPeriods + + +class MRClimatologyByGaussInterp(MRJob): + + def mapper(self, _, line): + yield "chars", len(line) + yield "words", len(line.split()) + yield "lines", 1 + + def reducer(self, key, values): + yield key, sum(values) + + +if __name__ == '__main__': + MRWordFrequencyCount.run() + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/jobTest.py ---------------------------------------------------------------------- diff --git a/climatology/clim/jobTest.py b/climatology/clim/jobTest.py new file mode 100644 index 0000000..5b1b0f7 --- /dev/null +++ b/climatology/clim/jobTest.py @@ -0,0 +1,18 @@ + +from mrjob.job import MRJob + + +class MRWordFrequencyCount(MRJob): + + def mapper(self, _, line): + yield "chars", len(line) + yield "words", len(line.split()) + yield "lines", 1 + + def reducer(self, key, values): + yield key, sum(values) + + +if __name__ == '__main__': + MRWordFrequencyCount.run() + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/README ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/README b/climatology/clim/orig/C/README new file mode 100644 index 0000000..7e50ce0 --- /dev/null +++ b/climatology/clim/orig/C/README @@ -0,0 +1,6 @@ +2 Mar 00 + +`test1-100day' was produced by gaussinterp (C version) using /sst/vol7/PATHFINDER/avhrrlists/test.list +and /sst/vol7/PATHFINDER/parmfiles/interp.test.parms. First 100 files in test.list were processed. + + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/binsum.c ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/binsum.c b/climatology/clim/orig/C/binsum.c new file mode 100644 index 0000000..724245c --- /dev/null +++ b/climatology/clim/orig/C/binsum.c @@ -0,0 +1,125 @@ +/* $Revision: 1.6 $ */ +/* calculate interpolation weights and vsum */ + +#include <math.h> +#include "interp.h" + +#define min(x,y) ( (x) < (y) ? (x) : (y) ) +#define max(x,y) ( (x) > (y) ? (x) : (y) ) +#define XFAC -0.6931 + + binsum( int nf ) { + + register int i, ii, jj, kk, n; + int lon_180_index; + float fac; + float long_diff_deg, long_diff_rad; + float west_long_endpoint, east_long_endpoint; + + /* index in array xx (longitude interp bins) for long. 180 + in extended portion of array beyond imax */ + if( all_zones ) + lon_180_index = parms.imax + long_extend / 2; + + /* determine spatial/temporal interpolation bins for every point */ + for( n=0; n < nf; n++ ) { + + /* longitude bins ranges */ + /* Check for the discontinuities along long. 180. If crossing long. 180 + set ixmin and ixmax to reflect indicies in "extended" portion of + longitude bins in array xx */ + + west_long_endpoint = x[n] - parms.xwin2; + east_long_endpoint = x[n] + parms.xwin2; + + ixmin[n] = (int) floorf( (x[n]-parms.xwin2-parms.xlo)/parms.dx ); + ixmax[n] = (int) floorf( (x[n]+parms.xwin2-parms.xlo)/parms.dx ); + + /* crossing west of 180 in search and global interp */ + if( west_long_endpoint < -180. && east_long_endpoint > -180. && + all_zones ) { + + ixmin[n] = lon_180_index - abs( ixmin[n] ); + ixmax[n] = lon_180_index + ixmax[n]; +/* printf( " ixmin, ixmax %d %d in extended array\n", ixmin[n], ixmax[n] ); */ + } + + /* crossing east of 180 in search and global interp */ + else if( west_long_endpoint < 180. && east_long_endpoint > 180. && + all_zones ) { + + ixmin[n] = lon_180_index - ( parms.imax - ixmin[n] ); + ixmax[n] = lon_180_index + abs( parms.imax - ixmax[n] ); +/* printf( "ixmin, ixmax %d %d in extended array\n", ixmin[n], ixmax[n] ); */ + } + + /* no crossing of 180 in search */ + else { + ixmin[n] = max( (int) floorf( (x[n]-parms.xwin2-parms.xlo)/parms.dx ), 0 ); + ixmax[n] = min( (int) floorf( (x[n]+parms.xwin2-parms.xlo)/parms.dx ), parms.imax - 1 ); + } + + /* lat bin ranges */ + iymin[n] = max( (int) floorf( (y[n]-parms.ywin2-parms.ylo)/parms.dy ), 0 ); + iymax[n] = min( (int) floorf( (y[n]+parms.ywin2-parms.ylo)/parms.dy ), parms.jmax - 1 ); + + /* temporal bin ranges */ + izmin[n] = max( (int) floorf( (z[n]-parms.zwin2-parms.zlo)/parms.dz ), 0 ); + izmax[n] = min( (int) floorf( (z[n]+parms.zwin2-parms.zlo)/parms.dz ), parms.kmax - 1 ); + + } + + + /* weight the SST and sum */ + for( n=0; n < nf; n++ ) { + for( kk=izmin[n]; kk <= izmax[n]; kk++ ) { + for( jj=iymin[n]; jj <= iymax[n]; jj++ ) { + for( i=ixmin[n]; i <= ixmax[n]; i++ ) { + + /* Derive the guassian weights. Because + longitudes can cross 180, e.g., the ranges + can be 175 to -175, use the acos(cos( long_diff ))) */ + long_diff_deg = x[n] - xx[i]; + long_diff_rad = acos( cos(deg2rad(long_diff_deg)) ); + + fac = exp( XFAC * ( pow( (rad2deg(long_diff_rad) + /parms.xh),2 ) + + pow( ((y[n]-yy[jj])/parms.yh),2 ) + + pow( ((z[n]-zz[kk])/parms.zh),2 ) )); + + /* if crossing long. 180 must drop weighted SST + in correct longitude bin between 0 and imax */ + + /* range crossing east of 180 */ + if( i >= lon_180_index && all_zones ) { + ii = i - lon_180_index; + +/* printf(" orig bin is %d, correct bin is %d\n", i, ii ); */ + } + /* range crossing west of 180 */ + else if( (i < lon_180_index && i >= parms.imax) && all_zones ) { + ii = i - long_extend / 2; + +/* printf(" orig bin is %d, correct bin is %d\n", i, ii ); */ + } + /* no crossing */ + else { + ii = i; + } + + vsum[0][ii][jj][kk] += f[n]*fac; + vsum[1][ii][jj][kk] += fac; + } + } + } + } + } + + + + + + + + + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/clouderosion.c ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/clouderosion.c b/climatology/clim/orig/C/clouderosion.c new file mode 100644 index 0000000..6dec692 --- /dev/null +++ b/climatology/clim/orig/C/clouderosion.c @@ -0,0 +1,33 @@ +/* $Revision:$ */ +/* Apply the Casey cloud erosion to the satellite data. + This algorithm flags as cloudy all data with one + nearest neighbor of a Pathfinder flagged cloud pixel. */ + +#include "interp.h" +#include "hdf.h" + +clouderosion( float in_data, int i, int j ) { + int row, col; + + /* for each image pixel loop over + nearest neighbors; if any neightbor is + a cloud, set temp[j][i] to cloud value (0) + */ + for( col = i-1; col <= i+1; col++ ) { + for( row = j-1; row <= j+1; row++ ) { + printf( "col is %d, row is %d\n", col, row ); + /* make sure we are searching within image array bounds */ + if( col >= 0 && col < parms.x_length && + row >= 0 && row < parms.y_length ) { + printf( "valid col is %d, valid row is %d\n\n", col, row ); + if( in_data[col][row] == 0 ) { + temp[j][i] = 0.; + printf( "in_data[%d][%d] is %c\n", col, row, in_data[col][row] ); + printf( " setting temp[%d][%d] to zero\n\n\n", j,i ); + break; + } + } + } + } + +} http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/gaussinterp.readme ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/gaussinterp.readme b/climatology/clim/orig/C/gaussinterp.readme new file mode 100644 index 0000000..4a0e017 --- /dev/null +++ b/climatology/clim/orig/C/gaussinterp.readme @@ -0,0 +1,159 @@ +05 Oct 00 +ed armstrong NASA/JPL/Caltech + +`gaussinterp' is a program designed to spatially and temporally +interpolate SST pathfinder data using gaussian methodology. This +description is for the C version. To run the program from the command +line: + +% gaussinterp <infile> <beg_rec> <parmfile> <outfile> + +where + infile = file containing list of files to process + (ascii; fixed length records; temporally sorted from + old to recent) + beg_rec = first record (file) in infile processed + parmfile = input parameter file + outfile = interpolated output file (ascii or binary) + + +`infile' must be ascii with a fixed record length of 150 bytes. Use +`pad.pl' to convert variable or other fixed record length files to +150. `infile' records should be in monotonically ascending in time +(year and julian day). beg_rec is the record number of the temporal +start of the interpolation. The day for the file corresponding to this +record number is also the first map center date. The next map center is +'first_map_center' + time_step (parameter dz; see parameter file +description) and so forth. `outfile' can be either an ascii or binary +file, sequential or interleaved format; see below for more +information. `parmfile' is an ascii file read by the program that +contains many interpolation and output parameters. + + +** Parameter file description ** + +~~~~ example of parameter file: ~~~~~ +200000 +720 360 +10 +3.0 3.0 20.0 +1.0 1.0 10.0 +1.0 1.0 8.0 +-180. 180. +-90. 90. +sequential +unformatted +1 +1 +6 +~~~~ end example ~~~~ + +record 1: ndatamax = maximum points in each file + 2: xlength, y_length = HDF data array: cols and rows + 3: kmax = number of "maps" or time "steps" + 4: xwin2, ywin2, zwin2 = longitude, latitude "search" + range (decimal degrees), time search range (days) + 5: xh, yh, zh = e folding scaling for long, lat (decimal + degrees) and time (days) + 6: dx, dy, dz = longitude, latitude interpolation + step (decimal degrees), and + time step (days) + 7: xlo, xhi = minimum, maximum longitude (-180 to 180 decimal degrees) + 8: ylo, yhi = minimum, maximum latitude (-90 to 90 decimal degrees) + 9: 'interleave' or 'sequential' map format + 10: 'formatted' (ascii) or 'unformatted' (binary) output file + 11: cloud erosion flag: 1=cloud erosion; 0=no cloud erosion + 12: all pixel data flag: 1="all pixel data"; 0="best pixel" data + 13: quality flag value: quality flag value (0-7) to use for "all pixel" data; + ignored if "best pixel" data used + + +Some notes on the parameter file: + + * ndatamax should be x_length*y_length but in reality is much + less since much of the data are missing due to clouds. + + * the user can subset into the image array by choosing proper + xlo, xhi, ylo, yhi. xlo must always be less than xhi; the same + for ylo and yhi. + + * interleave format refers to each row in the output file + representing a unique lon, lat position with columns of + associated sst values (maps) + i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst + row two: lon2 lat1 sst1 sst2 sst3....lastsst + + * sequential format refers to each map written as rows and + columns of lat, lon with the array element representing the sst + at that geo-position for that map. Each map is written + sequentially in the output file. + + * generally unformatted (binary) output should be chosen as + this results in a smaller output file size. Using binary format + each pixel is written as a 1 byte char (0-255). + + * cloud erosion refers to the Casey cloud erosion algorithm + used in the NSIPP (Casey) Pathfinder climatology. + + +** Header description ** + +For each output format a header is written: + ASCII format: + record 1: some text + record 2: imax, jmax, kmax + imax = width of output in pixels + jmax = height of output in pixels + kmax = number of "maps" or images + record 3: xwin2, ywin2, zwin2 + spatial and temporal search "radii" + in decimal degrees and days + record 4: xh, yh, zh + e folding scales in decimal degrees + record 5: dx,dy,dz + interpolation intervals + in decimal degrees + record 6: ULx, ULy, LRx, LRy + upper left and lower right longitude + and latitudes in decimal degrees + record 7: date.startdate, date.enddate + first and last map center dates in + YYYYDDD (e.g., 1995010) + + BINARY (C) [word is 4 bytes, types are indicated]: + word 1 [char]: magic number + words 2-4 [int]: imax, jmax, kmax + words 5-7 [float]: xwin2, ywin2, zwin2 + words 8-10 [float]: xh, yh, zh + words 11-13 [float]: dx,dy,dz + words 14-17 [float]: ULx, ULy, LRx, LRy + word 18 [int]: first map center date + word 19 [int]: last map center date + [total header size is 76 bytes] + +** Additional information ** + +- The parameter file can be edited to change/update all entries but +care must be taken to avoid typos. The program is compiled for a fixed +record length for the input datafile. These can be modified in +`interp.h' and the program recompiled. + +- In the case of global interpolation where xlo=-180 and xhi=180, the +program will correctly interpolate across longitude 180 East (all_zones +is true). That is, there is no "cutoff" in the interpolation if the +search window for a particular satellite pixel crosses this longitude. +In all other cases this is not so. For example, if xlo=-180 and +xhi=-150, interpolations in the vicinity of longitude -180 are "cutoff" +at this boundary, i.e., points west of -180 do not contribute to the +interpolation. Since this is a regional subset, the interpolations are +also "cutoff" in the vicinity of longitude -150 (data east of -150 do +not contribute). By extension this is also true for a subset in +latitude. + +- The first line in the usage dump (type gaussinterp with no arguments) +will specify if the executable is from Fortran or C code. + +- Sample parameter files can be found in +/sst/vol7/PATHFINDER/interp/parmfiles . Datafile lists can be found in +/sst/vol7/PATHFINDER/interp/avhrrlists . + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/gaussinterp_C_code.tar ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/gaussinterp_C_code.tar b/climatology/clim/orig/C/gaussinterp_C_code.tar new file mode 100644 index 0000000..5db4b0e Binary files /dev/null and b/climatology/clim/orig/C/gaussinterp_C_code.tar differ http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/interp.c ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/interp.c b/climatology/clim/orig/C/interp.c new file mode 100644 index 0000000..c057875 --- /dev/null +++ b/climatology/clim/orig/C/interp.c @@ -0,0 +1,448 @@ +/* $Revision: 1.20 $ */ + +/* SYNOPSIS + called by setupinterp(). This subroutine reads each + HDF file, passes data to binsum() for interpolation + calculation and writes output image(s). +*/ + +#include <stdio.h> +#include <string.h> +#include <libgen.h> +#include <hdf.h> +#include "interp.h" + +/* C binary file id = 3890345459 (0xE7E1F5F3) */ +#define MAGIC 3890345459L + + interp() { + uint8 in_data[parms.y_length][parms.x_length]; + uint8 quality_data[parms.y_length][parms.x_length]; + uint8 aerosol_corr[parms.y_length][parms.x_length]; + + intn retn, haspal; + int32 width, height; + + float ULlon, ULlat, LRlon, LRlat; + float cal, offset, xlat; + int k1,k2,k3,k4,ii,jj,kk,n,i,j,k; + int ipts, iyr, iday; + int jj_reverse; + int row, col; + unsigned long int magic_num; + + FILE *fid_infile, *fid_outfile, *fid_datafile; + FILE *fid_aerosol_list, *fid_aerosol_file; + + char plabel[40]; + char datafile[RECSIZE]; + char aerosol_file[200]; + char *basefile; + char cday[3]; + char cyr[4]; + char tmpstr[22]; + char aerosol_tmpstr[30]; + char tmpfile[35]; + char aerosol_tmpfile[45]; + char zcatcmd[200]; + char rmcmd[45]; + char compress[1]; + char version[4]; + char sstDN; + /* short int sstDN; */ + + /* open output file */ + fid_outfile = fopen( fileio.outfile, "w" ); + + + /* Set arrays of interpolation bins. + Long and lat bin centers set to 1/2 + interpolation step */ + + /* Long to interpolate to: + 1st for loop sets interpolation bins + between the ranges of xlo to xhi. + If all_zones is true, 2nd and 3rd for + loops set bins for those "extended" interpolations + that must cross longitude 180. */ + for( ii=parms.imin; ii < parms.imax; ii++ ) + xx[ii] = parms.xlo + (float)(ii)*parms.dx + parms.dx/2; + + if( all_zones ) { + /* west of long. 180 */ + for( ii=0; ii < long_extend/2; ii++ ) + xx[ii+parms.imax] = 180 - (long_extend/2 * parms.dx) + (float)(ii)*parms.dx + + parms.dx/2; + + /* east of long. 180; sets center point */ + for( ii=0; ii <= long_extend/2; ii++ ) + xx[ii+parms.imax+long_extend/2] = -180 + (float)(ii)*parms.dx + + parms.dx/2; + } + + /* for( ii = 0; ii <= parms.imax + long_extend; ii++) + printf( " ii: %d xx[ii] = %f \n", ii, xx[ii] ); */ + + + /* Lat to interpolate to */ + for( jj=parms.jmin; jj < parms.jmax; jj++ ) + yy[jj] = parms.ylo + (float)(jj)*parms.dy + parms.dy/2; + + /* Time to interpolate to */ + for( kk=parms.kmin; kk < parms.kmax; kk++ ) + zz[kk] = parms.zlo + (float)(kk)*parms.dz; + + /* slope (cal) and offset to convert DN to SST */ + cal = 0.15; + offset = 3.0; + + /* open input data file list */ + if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){ + printf( "\nCan not open %s\n", fileio.infile ); + exit(1); + } + + /* position to proper record in file */ + fseek( fid_infile, ( (fileio.istart-1)*RECSIZE ), 0 ); + + + + /* open input aerosol file list if necessary */ + if(aerosol) { + if( (fid_aerosol_list = fopen( fileio.aerosol_list, "r" )) == NULL ){ + printf( "\nCan not open %s\n", fileio.aerosol_list ); + exit(1); + } + /* position to proper record in file */ + fseek( fid_aerosol_list, ( (fileio.istart-1)*RECSIZE ), 0 ); + + } + + + /* read infile one record at a time and process . . . */ + for( k=fileio.istart; k <= fileio.iend; k++ ){ + + ipts = 0; + strcpy( compress, " " ); + strcpy( tmpstr, " " ); + + fgets( datafile, sizeof(datafile)+1, fid_infile ); + + /* strip off newline if necessary */ + if( datafile[ strlen(datafile) ] == '\n' ){ + datafile[ strlen(datafile) ] = '\0'; + } + /* printf("datafile is %s \n", datafile ); */ + + /* check to make sure datafile actually exists + if not ..... quit program */ + fid_datafile = fopen( datafile, "r" ); + if( fid_datafile == NULL ) { + printf( "could not open:\n" ); + printf( "%s \n", datafile ); + printf( ". . . quitting program \n" ); + exit(2); + } + fclose( fid_datafile ); + + basefile = basename( &datafile[0] ); + printf( "\nfile: %s\n", basefile ); + + /* e.g, parse year and jul day from 1985004h54dd-gdm.hdf */ + strcpy( tmpstr, basefile ); + strncpy( cyr, &tmpstr[0], 4 ); + iyr = atoi( cyr ); + strncpy( cday, &tmpstr[4], 3 ); + iday = atoi( cday ); + /* printf(" iday is %d\n", iday ); + printf("iyr is %d\n", iyr ); + */ + + /* if unix compressed, zcat to a temporary filename */ + strncpy( compress, &tmpstr[21], 1 ); + if( strcmp( compress, "Z" ) == 0 ) { + strcpy( tmpstr, basefile ); + strcpy( tmpfile, "/tmp/" ); + strncat( tmpfile, &tmpstr[0], 20 ); + + strcpy( zcatcmd, "zcat "); + strcat( zcatcmd, datafile ); + strcat( zcatcmd, " > " ); + strcat( zcatcmd, tmpfile ); + + (void) system( zcatcmd ); + strcpy( datafile, tmpfile ); + } + + /* do the same for the aerosol file if necessary */ + if( aerosol ) { + strcpy( compress, " " ); + strcpy( aerosol_tmpstr, " " ); + + fgets( aerosol_file, sizeof(aerosol_file)+1, fid_aerosol_list ); + /* printf("aerosol_file is %s \n", aerosol_file ); */ + + /* strip off newline if necessary */ + if( aerosol_file[ strlen(aerosol_file) ] == '\n' ){ + aerosol_file[ strlen(aerosol_file) ] = '\0'; + } + + + basefile = basename( &aerosol_file[0] ); + strcpy( aerosol_tmpstr, basefile ); + printf( "aerosol file: %s\n", basefile ); + + + /* if unix compressed, zcat to a temporary filename */ + strncpy( compress, &aerosol_tmpstr[27], 1 ); + if( strcmp( compress, "Z" ) == 0 ) { + strcpy( aerosol_tmpstr, basefile ); + strcpy( aerosol_tmpfile, "/tmp/" ); + strncat( aerosol_tmpfile, &aerosol_tmpstr[0], 26 ); + + strcpy( zcatcmd, "zcat "); + strcat( zcatcmd, aerosol_file ); + strcat( zcatcmd, " > " ); + strcat( zcatcmd, aerosol_tmpfile ); + + (void) system( zcatcmd ); + strcpy( aerosol_file, aerosol_tmpfile ); + } + } + + /* read the aerosol file into an array */ + if( aerosol ) { + fid_aerosol_file = fopen( aerosol_file, "r" ); + fread( aerosol_corr, sizeof(uint8), parms.y_length * parms.x_length, fid_aerosol_file ); + fclose( fid_aerosol_file ); + } + + /* read HDF SST data dimensions */ + retn = DFR8getdims( datafile, &width, &height, &haspal ); + /* printf( "HDF width and height: %d %d \n", width, height ); */ + + if( width != parms.x_length || height != parms.y_length ) { + printf( "\nbounds mismatch after reading HDF file dimensions !!\n" ); + printf( "width: %d, parms.x_length: %d \n", width, parms.x_length ); + printf( "height: %d, parms.y_length: %d \n", height, parms.y_length ); + printf( "skipping %s . . . \n\n", datafile ); + continue; + } + + /* ------------------------------------------------- */ + /* Read the SST image of the HDF file. + HDF call: DFR8getimage() */ + /* ------------------------------------------------- */ + + retn=DFR8getimage(datafile, (VOIDP)in_data, parms.x_length, + parms.y_length, NULL); + + /* ------------------------------------------------- */ + /* If "all pixel" data choosen read the quality + flag image */ + /* ------------------------------------------------- */ + + if( sst_type ) + retn=DFR8getimage(datafile, (VOIDP)quality_data, + parms.x_length, parms.y_length, NULL); + + + + /* ------------------------------------------------- */ + /* Read hdf file DN values and convert to SST. + If "all pixel" data read only those pixels + with a quality flag greater than or equal to + quality_flag */ + /* ------------------------------------------------- */ + + for( i=hdf.icolLeft; i <= hdf.icolRight; i++ ) { + for( j=hdf.irowUpper; j <= hdf.irowLower; j++ ) { + + /* only conider aerosol correction greater than 5 cnt */ + if( aerosol && aerosol_corr[j][i] > 5 && in_data[j][i] > 0 ) + in_data[j][i] = in_data[j][i] + aerosol_corr[j][i]; + + /* perform Casey cloud erosion on non-zero SST data */ + /* if nearest neighbor is a cloud, set pixel to cloud (0) */ + if( in_data[j][i] > 0 && cloud_erosion ) { + for( col = i-1; col <= i+1; col++ ) { + for( row = j-1; row <= j+1; row++ ) { + /* make sure we are searching within + image array bounds */ + if( col >= 0 && col < parms.x_length && + row >= 0 && row < parms.y_length && + in_data[row][col] == 0 ) { + in_data[j][i] = 0; + goto end_erosion; + } + } + } + } + end_erosion: + + + /* if all_pixel (sst_type = 1) filter for quality level, + or if best_sst accept data values greater than 0 */ + + /*printf( "quality flag is %d \n", quality_flag); + printf (" sst_type is %d \n", sst_type); + */ + + if( (sst_type && quality_data[j][i] >= quality_flag) || + (!sst_type && in_data[j][i] > 0) ) { + /* printf (" quality_data = %d \n", quality_data[j][i]); */ + + f[ipts] = cal * (float)in_data[j][i] - offset; + x[ipts] = -180. + parms.dxd * (float)(i) + parms.dxd/2; + y[ipts] = 90. - parms.dyd * (float)(j) - parms.dyd/2; + z[ipts] = (float)( iday+( (iyr-1985)*365 ) ); + + if(iyr > 1988) z[ipts]=z[ipts]+1; + if(iyr > 1992) z[ipts]=z[ipts]+1; + if(iyr > 1996) z[ipts]=z[ipts]+1; + if(iyr > 2000) z[ipts]=z[ipts]+1; + if(iyr > 2004) z[ipts]=z[ipts]+1; + if(iyr > 2008) z[ipts]=z[ipts]+1; + ipts++; + } + } + } + + /* calculate interpolation weights and vsum */ + + printf(" num of pts in file: %d \n", ipts ); + binsum( ipts ); + + if( strcmp( compress, "Z" ) == 0 ) { + strcpy( rmcmd, "rm -f " ); + strcat( rmcmd, tmpfile ); + system( rmcmd ); + if( aerosol ) { + strcpy( rmcmd, "rm -f " ); + strcat( rmcmd, aerosol_tmpfile ); + system( rmcmd ); + } + } + /* ..... read next hdf file from infile */ + + } + fclose( fid_infile ); + if( aerosol ) + fclose( fid_aerosol_list ); + + /* all input files processed; calculate interpolated SST maps */ + for( kk=0; kk < parms.kmax; kk++ ) { + for( jj=0; jj < parms.jmax; jj++ ) { + for( ii=0; ii < parms.imax; ii++ ) { + if ( vsum[1][ii][jj][kk] > 0 ) { + + /* reverse jj indicies to flip N to S */ + jj_reverse = parms.jmax - jj - 1; + sl[ii][jj_reverse][kk] = + vsum[0][ii][jj][kk] / vsum[1][ii][jj][kk]; + } + } + } + } + + /* write output as map "interleaved" or map "sequential". + both formats preceded by header. + see gaussinterp.readme for more info. + */ + + /* geo-positions of UL and LR corners */ + ULlon = parms.xlo + parms.dx/2; + ULlat = parms.yhi - parms.dy/2; + LRlon = parms.xhi - parms.dx/2; + LRlat = parms.ylo + parms.dy/2; + + /* version number, "c" refers to C version */ + strcpy( version, "c1.3" ); + strcpy( plabel, "sst interpolation"); + + /* write the header */ + if( strcmp( fileio.outformat, "formatted" ) == 0 ) { + fprintf( fid_outfile, "%s -- %s\n ", version, plabel ); + fprintf( fid_outfile, "%d %d %d\n", + parms.imax,parms.jmax,parms.kmax ); + fprintf( fid_outfile, "%f %f %f\n", + parms.xwin2,parms.ywin2,parms.zwin2 ); + fprintf( fid_outfile, "%f %f %f\n", + parms.xh,parms.yh,parms.zh ); + fprintf( fid_outfile, "%f %f %f\n", + parms.dx,parms.dy,parms.dz ); + fprintf( fid_outfile, "%f %f %f %f\n", + ULlon,ULlat,LRlon,LRlat ); + fprintf( fid_outfile, "%ld %ld\n", + date.startdate, date.enddate ); + } + else if( strcmp( fileio.outformat, "unformatted" ) == 0 ) { + magic_num = MAGIC; + fwrite( &magic_num, sizeof(magic_num),1,fid_outfile ); + fwrite( &parms.imax, sizeof(parms.imax),1,fid_outfile ); + fwrite( &parms.jmax, sizeof(parms.jmax),1,fid_outfile ); + fwrite( &parms.kmax, sizeof(parms.kmax),1,fid_outfile ); + fwrite( &parms.xwin2, sizeof(parms.xwin2),1,fid_outfile ); + fwrite( &parms.ywin2, sizeof(parms.ywin2),1,fid_outfile ); + fwrite( &parms.zwin2, sizeof(parms.zwin2),1,fid_outfile ); + fwrite( &parms.xh, sizeof(parms.xh),1,fid_outfile ); + fwrite( &parms.yh, sizeof(parms.yh),1,fid_outfile ); + fwrite( &parms.zh, sizeof(parms.zh),1,fid_outfile ); + fwrite( &parms.dx, sizeof(parms.dx),1,fid_outfile ); + fwrite( &parms.dy, sizeof(parms.dy),1,fid_outfile ); + fwrite( &parms.dz, sizeof(parms.dz),1,fid_outfile ); + fwrite( &ULlon, sizeof(ULlon),1,fid_outfile ); + fwrite( &ULlat, sizeof(ULlat),1,fid_outfile ); + fwrite( &LRlon, sizeof(LRlon),1,fid_outfile ); + fwrite( &LRlat, sizeof(LRlat),1,fid_outfile ); + fwrite( &date.startdate, sizeof(date.startdate),1,fid_outfile); + fwrite( &date.enddate, sizeof(date.enddate),1,fid_outfile); + /* binary header is 76 bytes */ + } + + /* write the maps */ + if( strcmp( fileio.mapformat, "interleave" ) == 0 ) { + printf( "\n map output is interleave\n" ); + for( jj=0; jj < parms.jmax; jj++ ) { + for( ii=0; ii < parms.imax; ii++ ) { + for( kk=0; kk < parms.kmax; kk++ ) { + sst_inter[kk] = sl[ii][jj][kk]; + + /* scale sst back to 0-255 DNs */ + if( sst_inter[kk] != 0 ) + sstDN = (char) nint( (sst_inter[kk] + offset) / cal ); + else + sstDN = 0; /* O's in array sst_inter are missing data */ + + if( strcmp( fileio.outformat, "formatted" ) == 0 ) + fprintf( fid_outfile, "%d ", sstDN ); + else + fwrite( &sstDN, sizeof(char),1,fid_outfile ); + + } + } + } + } else if( strcmp( fileio.mapformat, "sequential" ) == 0 ) { + printf( "\n map output is sequential\n" ); + for( kk=0; kk < parms.kmax; kk++ ) { + for( jj=0; jj < parms.jmax; jj++ ) { + for( ii=0; ii < parms.imax; ii++ ) { + sst_seq[ii] = sl[ii][jj][kk]; + + /* scale sst back to 0-255 DNs */ + if( sst_seq[ii] != 0 ) + sstDN = (char) nint( (sst_seq[ii] + offset) / cal ); + else + sstDN = 0; /* O's in array sst_seq are missing data */ + + if( strcmp( fileio.outformat, "formatted" ) == 0 ) + fprintf( fid_outfile, "%d ", sstDN ); + else + fwrite( &sstDN, sizeof(char),1,fid_outfile ); + + } + } + } + } + fclose( fid_outfile ); +} http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/makefile ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/makefile b/climatology/clim/orig/C/makefile new file mode 100644 index 0000000..ab004b3 --- /dev/null +++ b/climatology/clim/orig/C/makefile @@ -0,0 +1,33 @@ +# makefile for guassinterp (C version). + +# C compiler +CC = cc +#HDF = /usr/local/hdf +HDF = /usr/local/hdf4.1r3-n32 + +CFLAGS = -n32 -I$(HDF)/include -O -w +#LIBS = -L$(HDF)/lib -lmfhdf -ldf -ljpeg -lz -lgen -lm +LIBS = -L$(HDF)/lib -ldf -ljpeg -lz -lgen -lm +BIN = /sst/vol7/PATHFINDER/bin/sgi +PGM = gaussinterp + +SRC = setupinterp.c interp.c binsum.c +OBJ = setupinterp.o interp.o binsum.o + +$(PGM): $(OBJ) + $(CC) $(CFLAGS) -o $(PGM) $(OBJ) $(LIBS) + +setupinterp.o: + $(CC) $(CFLAGS) $(DEFINE) -c setupinterp.c +interp.o: + $(CC) $(CFLAGS) $(DEFINE) -c interp.c +binsum.o: + $(CC) $(CFLAGS) $(DEFINE) -c binsum.c + + +install: $(PGM) + cp $(PGM) $(BIN) + +clean: + rm -f $(OBJ) $(PGM) + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/C/setupinterp.c ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/C/setupinterp.c b/climatology/clim/orig/C/setupinterp.c new file mode 100644 index 0000000..9ca67b3 --- /dev/null +++ b/climatology/clim/orig/C/setupinterp.c @@ -0,0 +1,431 @@ +/* NAME gaussinterp (C version) + 16 Jan 00 earmstrong NASA/JPL + $Revision: 1.16 $ + + DESCRIPTION + Gaussian interpolation program modified to map SST (HDF) data. + This program is a C port of the fortran version of gaussinterp. + + SYNOPSIS + setupinterp() allocates space for arrays, calls interp() + which does the interpolation + + USAGE + % gaussinterp <infile> <beg_rec> <parmfile> <outfile> */ + + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <libgen.h> +#include <math.h> +#include "interp.h" +/* see interp.h for global structs and arrays */ + +main( int argc, char *argv[] ) { + + int iday, iyr, calloc_bytes, rec_step; + long max_infile_size; + + int value; + FILE *fid_infile, *fid_parm; + + char datafile[RECSIZE]; + char *basefile; + char *parmfile; + char *startrec, *endrec; + char cday[3], cyr[4], cdate[7]; + char *cstartdate; + char *cenddate; + char tmpline[80]; + +/* check command line arguments */ + if( argc != 5 && argc != 6 ) { + usage(); + exit(1); + } + fileio.infile = argv[1]; + startrec = argv[2]; + parmfile = argv[3]; + fileio.outfile = argv[4]; + if( argc == 6 ) + fileio.aerosol_list = argv[5]; + + /* open parameter file and read line by line */ + if( (fid_parm = fopen( parmfile, "r" )) == NULL ){ + printf( "\nCan not open %s\n", parmfile ); + exit(1); + } + + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &parms.ndatamax ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d %d", &parms.x_length, &parms.y_length ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &parms.kmax ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%f %f %f", &parms.xwin2, &parms.ywin2, &parms.zwin2 ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%f %f %f", &parms.xh, &parms.yh, &parms.zh ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%f %f %f", &parms.dx, &parms.dy, &parms.dz ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%f %f", &parms.xlo, &parms.xhi ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%f %f", &parms.ylo, &parms.yhi ); + fgets( fileio.mapformat, sizeof(fileio.mapformat), fid_parm ); + fgets( fileio.outformat, sizeof(fileio.outformat), fid_parm ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &cloud_erosion ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &sst_type ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &quality_flag ); + fgets( tmpline, sizeof(tmpline), fid_parm ); + sscanf( tmpline, "%d", &aerosol ); + + fileio.mapformat[ strlen(fileio.mapformat)-1 ] = '\0'; + fileio.outformat[ strlen(fileio.outformat)-1 ] = '\0'; + + fclose( fid_parm ); + + if( strcmp(fileio.mapformat, "interleave") != 0 && + strcmp(fileio.mapformat, "sequential") != 0 ) { + printf( "check parameter file for interleave or sequential maps !\n" ); + exit(1); + } + if( strcmp(fileio.outformat, "formatted") != 0 && + strcmp(fileio.outformat, "unformatted") != 0 ) { + printf( "check parameter file for formatted or unformatted output !\n" ); + exit(1); + } + + /* Calculate number of interpolation "steps" based on + long and lat ranges divided by interpolation interval. + Origin is at lon=-180, lat=90 */ + parms.imax = floorf((parms.xhi + 180) / parms.dx) - floorf((parms.xlo + 180) / parms.dx ); + parms.jmax = floorf((parms.yhi - 90) / parms.dy) - floorf((parms.ylo - 90) / parms.dy ); + printf( "imax, jmax, kmax: %d %d %d \n", + parms.imax, parms.jmax, parms.kmax ); + + parms.imin = 0; + parms.jmin = 0; + parms.kmin = 0; + + /* calculate degrees per grid point */ + parms.dxd = 360. / (float)parms.x_length; + parms.dyd = 180. / (float)parms.y_length; + + /* find columns and rows for subsetting into hdf image array + remember @ hdfarray(0,0) --> lon=-180, lat=90 */ + hdf.icolLeft = floorf( (parms.xlo + 180.) / parms.dxd ); + hdf.icolRight = floorf( ((parms.xhi + 180.) / parms.dxd) - 1. ); + hdf.irowUpper = floorf( fabsf(parms.yhi - 90.) / parms.dyd ); + hdf.irowLower = floorf( (fabsf(parms.ylo - 90.) / parms.dyd) - 1. ); + + /* printf( "array corners: %d %d %d %d \n", hdf.icolLeft, hdf.irowUpper, + hdf.icolRight ,hdf.irowLower); */ + + /* set global interpolation flag if necessary */ + if( parms.xlo == -180. && parms.xhi == 180. ) { + all_zones = 1; /* true */ + printf( "all_zones is true: interpolating across -180 West\n" ); + } + + if( cloud_erosion ) + printf( "performing cloud erosion filtering on SST data \n" ); + if( sst_type ) + printf( "working with 'all pixel' data: quality flag is %d \n", quality_flag ); + if( aerosol ) + printf( "performing aerosol correction \n" ); + + /* allocate space for arrays */ + calloc_bytes = 0; + if( (temp = fmal2d( parms.x_length, parms.y_length )) == NULL ) + printf( "not enough memory for array temp\n" ), exit(-1); + else calloc_bytes = parms.x_length * parms.y_length * sizeof(float); + + if( (in_data = cmal2d( parms.x_length, parms.y_length )) == NULL ) + printf( "not enough memory for array in_data\n" ), exit(-1); + else calloc_bytes += parms.x_length * parms.y_length * sizeof(char); + + if( (quality_data = cmal2d( parms.x_length, parms.y_length )) == NULL ) + printf( "not enough memory for array quality_data\n" ), exit(-1); + else calloc_bytes += parms.x_length * parms.y_length * sizeof(char); + + if( (sl = fmal3d( parms.imax, parms.jmax, parms.kmax )) == NULL ) + printf( "not enough memory for array sl\n" ), exit(-1); + else calloc_bytes += parms.imax * parms.jmax * parms.kmax * sizeof(float); + + if( (vsum = fmal4d( 2, parms.imax, parms.jmax, parms.kmax )) == NULL ) + printf( "not enough memory for array vsum\n" ), exit(-1); + else calloc_bytes += parms.imax * parms.jmax * parms.kmax * 2 * sizeof(float); + + if( (sst_seq = (float *)calloc( parms.imax, sizeof(float) )) == NULL ) + printf( "not enough memory for array sst_seq\n" ), exit(-1); + else calloc_bytes += parms.imax * sizeof(float); + + if( (sst_inter = (float *)calloc( parms.kmax, sizeof(float) )) == NULL ) + printf( "not enough memory for array sst_inter\n" ), exit(-1); + else calloc_bytes += parms.kmax * sizeof(float); + + if( (x = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL ) + printf( "not enough memory for array x\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(float); + + if( (y = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL ) + printf( "not enough memory for array y\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(float); + + if( (z = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL ) + printf( "not enough memory for array z\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(float); + + if( (f = (float *)calloc( parms.ndatamax, sizeof(float) )) == NULL ) + printf( "not enough memory for array f\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(float); + + /* Extend array xx by long_extend to handle crossing + longitude 180 during global interpolation. + Extend by 4*search radius (plus 1 in xx dimension + for -180 West crossing point */ + if( all_zones ) { + long_extend = nint( 4 * parms.xwin2 / parms.dx ); + /* make sure long_extend is even for equal number of + bins on either side of -180 W */ + if( long_extend % 2 != 0 ) + long_extend++; + if( (xx = (float *)calloc( parms.imax + long_extend + 1, sizeof(float) )) == NULL ) + printf( "not enough memory for array xx\n" ), exit(-1); + else calloc_bytes += (parms.imax + long_extend + 1) * sizeof(float); + } else { + if( (xx = (float *)calloc( parms.imax, sizeof(float) )) == NULL ) + printf( "not enough memory for array xx\n" ), exit(-1); + else calloc_bytes += parms.imax * sizeof(float); + } + + if( (yy = (float *)calloc( parms.jmax, sizeof(float) )) == NULL) + printf( "not enough memory for array yy\n" ), exit(-1); + else calloc_bytes += parms.jmax * sizeof(float); + + if( (zz = (float *)calloc( parms.kmax, sizeof(float) )) == NULL ) + printf( "not enough memory for array zz\n" ), exit(-1); + else calloc_bytes += parms.kmax * sizeof(float); + + if( (ixmin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array ixmin\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + if( (ixmax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array ixmax\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + if( (iymin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array iymin\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + if( (iymax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array iymax\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + if( (izmin = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array izmin\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + if( (izmax = (int *)calloc( parms.ndatamax, sizeof(int) )) == NULL ) + printf( "not enough memory for array izmax\n" ), exit(-1); + else calloc_bytes += parms.ndatamax * sizeof(int); + + + printf( "done calloc arrays\n" ); + printf( "total bytes used: %d (%5.1f Mb)\n\n", calloc_bytes, calloc_bytes/1e6 ); + + + /* Determine zlo which corresponds to record number in infile for + first map center and number of days since 1 Jan 1985 */ + fileio.istart = atoi( startrec ); + + /* Open input file list and read first record */ + if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){ + printf( "\nCan not open %s\n", fileio.infile ); + exit(1); + } + + /* position to proper record in file */ + fseek( fid_infile, (fileio.istart-1)*RECSIZE, SEEK_SET ); + fgets( datafile, sizeof(datafile), fid_infile ); + fclose( fid_infile ); + basefile = basename( &datafile[0] ); + + strcpy( cday, " " ); /* initialize strings */ + strcpy( cyr, " " ); + /* e.g, parse year and jul day from 1985004h54dd-gdm.hdf */ + strncpy( cyr, &basefile[0], 4 ); + iyr = atoi( cyr ); + strncpy( cday, &basefile[4], 3 ); + iday = atoi( cday ); + + /* set parms.zlo: days since 1 Jan 1985 */ + parms.zlo = (float)iday+(float)(iyr-1985)*365.; + if(iyr > 1988) parms.zlo = parms.zlo + 1.; + if(iyr > 1992) parms.zlo = parms.zlo + 1.; + if(iyr > 1996) parms.zlo = parms.zlo + 1.; + if(iyr > 2000) parms.zlo = parms.zlo + 1.; + if(iyr > 2004) parms.zlo = parms.zlo + 1.; + if(iyr > 2008) parms.zlo = parms.zlo + 1.; + printf( " zlo is %f\n", parms.zlo ); + printf( " first map centered on day: %d in year: %d\n", iday, iyr ); + + /* Set the first center dates */ + cstartdate = strcat( cdate, cyr ); + cstartdate = strcat( cdate, cday ); + date.startdate = atoi( cstartdate ); + printf( " start date is %d \n", date.startdate ); + + /* reopen file to determine last map center date */ + fopen( fileio.infile, "r" ); + fileio.iend = fileio.istart + (parms.kmax-1) * parms.dz; + fseek(fid_infile, (fileio.iend-1)*RECSIZE, SEEK_SET ); + fgets( datafile, sizeof(datafile), fid_infile ); + fclose( fid_infile ); + basefile = basename( &datafile[0] ); + + /* parse year and jul day */ + strcpy( cdate, "" ); + strncpy( cyr, &basefile[0], 4 ); + strncpy( cday, &basefile[4], 3 ); + cenddate = strcat( cdate, cyr ); + cenddate = strcat( cdate, cday ); + date.enddate = atoi( cenddate ); + printf( " end date is %d\n", date.enddate ); + + /* Now determine temporal interpolation (record) ranges */ + /* set istart and iend so the first and last maps will span a complete + temporal interpolation range if possible. e.g., if first map is centered + on 1 Jan 1986 (this is date of first record to read from the + input datafile [fileio.istart]), with a search radius of 5 days (zwin2), + set istart back five days to 27 Dec 1985 and iend to 6 Jan 1986 */ + + /* set initial step back and forward size to temporal search radius */ + rec_step = (int) parms.zwin2; + + /* reopen input datafile */ + if( (fid_infile = fopen( fileio.infile, "r" )) == NULL ){ + printf( "\nCan not open %s\n", fileio.infile ); + exit(1); + } + + /* determine max bytes in infile */ + fseek( fid_infile, 0L, SEEK_END ); + max_infile_size = ftell( fid_infile ); + + /* step back into infile; use return of fseek() to determine if + file pointer is set before beginning of file */ + while( fseek(fid_infile, (fileio.istart-1-rec_step)*RECSIZE, SEEK_SET) != 0 ) + rec_step--; + fileio.istart -= rec_step; + + /* step forward into infile */ + rec_step = parms.zwin2; + fseek( fid_infile, (fileio.iend-1+rec_step)*RECSIZE, SEEK_SET ); + + /* if the file pointer is set beyond the EOF: step back */ + while( ftell(fid_infile) >= max_infile_size ) { + rec_step--; + fseek( fid_infile, (fileio.iend-1+rec_step)*RECSIZE, SEEK_SET ); + } + fileio.iend += rec_step; + + printf( " records: istart %d and iend %d \n\n", fileio.istart, fileio.iend ); + fclose( fid_infile ); + + /* call interp() */ + interp(); + + free(temp), free(sl), free(vsum), free(sst_seq), free(sst_inter); + free(x), free(y), free(z), free(f), free(xx), free(yy), free(zz); + free(ixmin), free(ixmax), free(iymin), free(iymax); + free(izmin), free(izmax), free(in_data), free(quality_data); + } + +/* functions to allocate space for float arrays */ +/* 2d dimension */ +float **fmal2d( int ix, int iy ){ + + float **ival; + int i; + + ival = (float **)calloc( ix, sizeof(float *) ); + for( i=0; i < ix; i++ ) + ival[i] = (float *)calloc( iy, sizeof(float) ); + + return(ival); +} + + /* 3d dimension */ + float ***fmal3d( int ix, int iy, int iz ){ + + float ***ival; + int i, j; + + ival = (float ***)calloc( ix, sizeof(float **) ); + for( i=0; i < ix; i++ ) + ival[i] = (float **)calloc( iy, sizeof(float *) ); + + for( i=0; i < ix; i++ ) + for( j=0; j < iy; j++ ) + ival[i][j] = (float *)calloc( iz, sizeof(float) ); + + return(ival); + } + +/* 4d dimension */ +float ****fmal4d( int ix, int iy, int iz, int ia ){ + + float ****ival; + int i, j, k; + + ival = (float ****)calloc( ix, sizeof(float ***) ); + for( i=0; i < ix; i++ ) + ival[i] = (float ***)calloc( iy, sizeof(float **) ); + for( i=0; i < ix; i++ ) + for( j=0; j < iy; j++) + ival[i][j] = (float **)calloc( iz, sizeof(float *) ); + for( i=0; i < ix; i++ ) + for( j=0; j < iy; j++) + for( k=0; k < iz; k++) + ival[i][j][k] = (float *)calloc( ia, sizeof(float) ); + + return(ival); +} + +/* function to allocate space for char array */ +/* 2d dimension */ +char **cmal2d( int ix, int iy ) { + + char **ival; + int i; + + ival = (char **)calloc( ix, sizeof(char *) ); + for( i=0; i < ix; i++ ) + ival[i] = (char *)calloc( iy, sizeof(char) ); + + return(ival); +} + + + +usage() { + printf( " ~~ gaussian interpolation of pathfinder SST data (C version) ~~\n\n " ); + printf( "USAGE: gaussinterp <infile> <beg_rec> <parmfile> <outfile>\n" ); + printf( " -- infile: file containing list to process\n" ); + printf( " -- beg_rec: first record (file) in infile processed\n" ); + printf( " -- parmfile: input parameter file (e.g., interp.parms)\n" ); + printf( " -- outfile: interpolated output file (ascii or binary)\n" ); + printf( "\n e.g., gaussinterp list.dat 100 interp.parms output.dat\n" ); + printf( "[process records of list.dat starting from record 100; interpolated \ +result is output.dat]\n" ); + printf( "note: see interp.parms for example format of parameter file\n" ); + printf( " see gaussinterp.readme for more general information\n" ); + return(0); +} http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/armstrong_interp_code.tar ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/Fortran/armstrong_interp_code.tar b/climatology/clim/orig/Fortran/armstrong_interp_code.tar new file mode 100755 index 0000000..4e125b3 Binary files /dev/null and b/climatology/clim/orig/Fortran/armstrong_interp_code.tar differ http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/binsum.f ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/Fortran/binsum.f b/climatology/clim/orig/Fortran/binsum.f new file mode 100644 index 0000000..96b65af --- /dev/null +++ b/climatology/clim/orig/Fortran/binsum.f @@ -0,0 +1,64 @@ +c $Revision: 2.4 $ +c subprogram binsum + +c-- calculate interpolation weights and vsum + + subroutine binsum( nf, x, y, z, f, + & xx, yy, zz, vsum, + & ixmin, ixmax, iymin, + & iymax, izmin, izmax, + & imax, jmax, kmax, ndatamax ) + + real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax) + real*4 xx(imax),yy(jmax),zz(kmax) + real*4 vsum(2,imax,jmax,kmax) + integer*4 ixmin(ndatamax),ixmax(ndatamax) + integer*4 iymin(ndatamax),iymax(ndatamax) + integer*4 izmin(ndatamax),izmax(ndatamax) + + common /parms/ imin,jmin,kmin, + & xwin2,ywin2,zwin2, + & xh,yh,zh, + & dx,dy,dz, + & xlo,xhi,ylo,yhi,zlo, + & dxd,dyd + + data xfac/-0.6931/ + + do n=1,nf + ixmin(n)=max(nint((x(n)-xwin2-xlo+0.5*dx)/dx),1) + ixmax(n)=min(nint((x(n)+xwin2-xlo+0.5*dx)/dx),imax) + iymin(n)=max(nint((y(n)-ywin2-ylo+0.5*dy)/dy),1) + iymax(n)=min(nint((y(n)+ywin2-ylo+0.5*dy)/dy),jmax) + izmin(n)=max(nint((z(n)-zwin2-zlo+0.5*dz)/dz),1) + izmax(n)=min(nint((z(n)+zwin2-zlo+0.5*dz)/dz),kmax) +c print *, x(n),y(n),z(n), f(n) +c print *,' ixmin, ixmax', ixmin(n), ixmax(n) +c print *,' iymin, iymax', iymin(n), iymax(n) +c print *,' izmin, izmax', izmin(n), izmax(n) + enddo + + + + do n=1,nf + do kk=izmin(n),izmax(n) + do jj=iymin(n),iymax(n) + do ii=ixmin(n),ixmax(n) + +c- - this is the algorithm coded for weights + + fac=exp( xfac*(((x(n)-xx(ii))/xh)**2 + & + ((y(n)-yy(jj))/yh)**2 + & + ((z(n)-zz(kk))/zh)**2) ) + +c print *, 'x, xx, y, yy, z, zz, fac f', +c & x(n), xx(ii), y(n), yy(jj), z(n), zz(kk), fac, f(n) + + vsum(1,ii,jj,kk)=vsum(1,ii,jj,kk)+f(n)*fac + vsum(2,ii,jj,kk)=vsum(2,ii,jj,kk)+fac + enddo + enddo + enddo + enddo + return + end http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/interp.f ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/Fortran/interp.f b/climatology/clim/orig/Fortran/interp.f new file mode 100644 index 0000000..8977798 --- /dev/null +++ b/climatology/clim/orig/Fortran/interp.f @@ -0,0 +1,302 @@ +#include "interp.h" +c $Revision: 2.12 $ +c subprogram interp.f +c reads hdf file, extract data, call bisum(), calculates +c and writes out interpolated maps + + + subroutine interp( temp, sst_inter, sst_seq, sl, + & x, y, z, f, + & xx, yy, zz, vsum, ixmin, ixmax, + & iymin, iymax, izmin, izmax, ndatamax, + & x_length, y_length, imax, jmax, kmax ) + + integer x_length,y_length + integer dummy + +c-- dimension the arrays + real*4 temp(x_length,y_length) + real*4 sl(imax,jmax,kmax) + real*4 sst_inter(kmax), sst_seq(imax) + real*4 x(ndatamax),y(ndatamax),z(ndatamax),f(ndatamax) + real*4 xx(imax),yy(jmax),zz(kmax) + real*4 vsum(2,imax,jmax,kmax) + real*4 ULlon, ULlat, LRlon, LRlat + integer*4 ixmin(ndatamax),ixmax(ndatamax) + integer*4 iymin(ndatamax),iymax(ndatamax) + integer*4 izmin(ndatamax),izmax(ndatamax) + +c-- hdf data array must be hardwired, no dync alloc :( + byte in_data( XSIZE, YSIZE ) + + character*80 plabel + character*50 infile + character*30 outfile + character*150 datafile + character*30 basename + character*35 tmpfile + character*10 mapformat + character*15 outformat + character*3 cday + character*4 cyr + character*200 zcatcmd + + common /parms/ imin,jmin,kmin, + & xwin2,ywin2,zwin2, + & xh,yh,zh, + & dx,dy,dz, + & xlo,xhi,ylo,yhi,zlo, + & dxd,dyd + + common /fileio/ infile,outfile,istart, + & iend,mapformat,outformat + + common /hdfindex/ icolLeft,icolRight,irowUpper,irowLower + + common /basefile/ basename + +c-- open output file + open(unit=20,status='new',form=outformat,file=outfile, + & err=1000,iostat=ierr ) + plabel='sst interpolation' + +c-- initialize vsum to '0's + do k1=1,2 + do k2=1,kmax + do k3=1,jmax + do k4=1,imax + vsum(k1,k4,k3,k2)=0. + enddo + enddo + enddo + enddo + +c-- initialize arrays of long, lat and time intervals + do ii=imin,imax + xx(ii)=xlo+float(ii-1)*dx ! Long to interpolate to + enddo + do jj=jmin,jmax + yy(jj)=ylo+float(jj-1)*dy ! Lat to interpolate to + enddo + do kk=kmin,kmax + zz(kk)=zlo+float(kk-1)*dz ! Time to interpolate to + enddo + + do n=1,ndatamax + x(n)=0.0 + y(n)=0.0 + z(n)=0.0 + f(n)=0.0 + ixmin(n) = 0 + ixmax(n) = 0 + iymin(n) = 0 + iymax(n) = 0 + izmin(n) = 0 + izmax(n) = 0 + enddo + +c-- slope (cal) and offset to convert DN to SST + cal= 0.15 + offset = 3.0 + +c-- Open input file list + open(unit=15,status='old',file=infile,access='direct', + & recl=RECSIZE,form='formatted',err=1500,iostat=ierr) + +c-- read a infile one record at a time and process . . . + do k=istart,iend + ipts=0 + read( 15,'(a)',rec=k ) datafile +c print *, datafile + + icompress = 0 +c-- hack to get trailing 'Z' (if compressed) +c do islash = 150, 1, -1 +c if ( datafile(islash:islash) .eq. 'Z' ) then +c icompress = 1 +c goto 101 +c endif +c enddo +c101 continue + +c-- 12/27/99: C call for basename implemented +c-- variable basename returned via common statement + call getbase( datafile ) + print *,'\n file: ', basename + +c-- cyr and cday parsed from basename +c-- (e.g., parsed from '1988332h54dd-gdm.hdf') + cday=basename( 5:7 ) + cyr=basename( 1:4 ) + +c-- if unix compressed, zcat to a temporary filename + if( basename(22:22) .eq. 'Z' ) then + icompress = 1 + tmpfile = '/tmp/' // basename(1:20) + zcatcmd = 'zcat ' // datafile( 1:len(datafile)-1 ) + & // ' > ' // tmpfile + call system( zcatcmd ) + datafile = tmpfile + endif + +c-- convert iday character to integer + read(cday,'(i3)') iday + read(cyr,'(i4)') iyr +c write(6,*)'cday = ',cday +c write(6,*)'cyr = ',cyr + +c***> HDF call: d8gimg() + retn=d8gimg(datafile,in_data,x_length, + & y_length,dummy) + +c-- read hdf DN values and convert to SST + do i=icolLeft, icolRight + do j=irowUpper, irowLower + xlat=-89.75+dyd*float(j-1) + +c-- center output in Pacific Ocean +c if( i .lt. x_length ) ix = i-(x_length/2) +c if( i .le. (x_length/2) ) ix = i+(x_length/2) + +c-- center output in Atlantic (default) + ix = i + +c-- convert signed byte to float +c if ( in_data(i,j) .lt. 0 .and. abs(xlat) .lt. 70. ) then + if ( in_data(i,j).lt. 0 ) then + temp(ix,j)=float( in_data(i,j) ) + 256 + else + temp(ix,j)=float( in_data(i,j) ) + endif + + + +C MULTIPLY THE PATHFINDER DIGITAL NUMBER BY THE CALIBRATION NUMBER (0.15) +C AND ADD THE OFFSET (-3.0) TO GET DEGREES CELSIUS + + if ( temp(ix,j).gt.0 ) then + ipts=ipts+1 + f(ipts)=( cal*temp(ix,j) ) - offset + x(ipts) = ( dxd*float(ix-1) ) - 180. + dxd/2 + y(ipts) = ( 90. - (dyd*float(j-1)) ) - dyd/2 + z(ipts)=float(iday)+float(iyr-1985)*365. + if(iyr.gt.1988) z(ipts)=z(ipts)+1 + if(iyr.gt.1992) z(ipts)=z(ipts)+1 + if(iyr.gt.1996) z(ipts)=z(ipts)+1 + if(iyr.gt.2000) z(ipts)=z(ipts)+1 + if(iyr.gt.2004) z(ipts)=z(ipts)+1 + endif + enddo + enddo + + nfnew=ipts + print *, ' no of pts in file ',' = ', nfnew + +c-- calculate interpolation weights and vsum +c-- arrays passed directly... a common statement +c-- does not seem to work. + call binsum( nfnew, x, y, z, f, + & xx, yy, zz, vsum, + & ixmin, ixmax, iymin, + & iymax, izmin, izmax, + & imax, jmax, kmax, ndatamax ) + + if( icompress .eq. 1 ) call system( 'rm -f ' // tmpfile ) +c-- ..... read next hdf file from infile + enddo + +c-- all input files processed; calculate interpolated SST + do kk=1,kmax + do jj=1,jmax + do ii=1,imax + sl(ii,jj,kk)=0. + if (vsum(2,ii,jj,kk).gt.0) then + sl(ii,jj,kk)=vsum(1,ii,jj,kk)/vsum(2,ii,jj,kk) + endif + enddo + enddo + enddo + +c-- write output as map "interleaved" or map "sequential" +c-- "interleaved" is the original implementation +c-- both formats preceded by header + +c-- "interleaved" refers to each row in the output +c-- file representing a unique lon, lat position with columns +c-- of associated sst values (maps) +c-- i.e. row one: lon1 lat1 sst1 sst2 sst3....lastsst +c-- row two: lon2 lat1 sst1 sst2 sst3....lastsst + +c-- "sequential" refers to each map written as rows and +c-- columns of lat, lon with the array element representing the +c-- sst at that geo-position for that map. +c-- each map is written sequentially in the file + +c-- geo-positions of UL and LR corners + ULlon = -(180 - ( ((icolLeft-1) * dxd) + dxd/2 )) + ULlat = (90 - ( ((irowUpper-1) * dyd) + dyd/2 )) + LRlon = -(180 - ( (icolRight * dxd) - dxd/2 )) + LRlat = (90 - ( (irowLower * dyd) - dyd/2 )) + +c-- version number, "f" refers to fortran version + version = 'f2.9' + +c-- write the 3 record header + if( outformat .eq. 'formatted' ) then + write(20,'(a)'), plabel + write(20,*) imax,jmax,kmax + write(20,*) dx,dy,dz + write(20,*) ULlon,ULlat,LRlon,LRlat + elseif( outformat .eq. 'unformatted' ) then + write(20) imax,jmax,kmax + write(20) dx,dy,dz + write(20) ULlon,ULlat,LRlon,LRlat + endif + + if( mapformat .eq. 'interleave' ) then + print *, '\n map output is interleave' + do jj=1,jmax + do ii=1,imax + do kk=1,kmax + sst_inter(kk)=sl(ii,jj,kk) + enddo + if( outformat .eq. 'formatted' ) then + write(20,*) ii,jj,sst_inter + else + write(20) ( sst_inter(i), i=1,kmax ) + endif + enddo + enddo + + else if( mapformat .eq. 'sequential' ) then + print *, '\n map output is sequential' + do kk=1,kmax + do jj=1,jmax + do ii=1,imax + sst_seq(ii)=sl(ii,jj,kk) + enddo + if( outformat .eq. 'formatted' ) then + write(20,*) jj,kk,sst_seq + else + write(20) ( sst_seq(i), i=1,imax ) + endif + enddo + enddo + endif + + print *, '\ndone. . . ' + close( 15,status='keep',err=2000,iostat=ierr ) + close( 20,status='keep',err=2500,iostat=ierr ) + stop + + 1000 print *, 'Error opening output: ', outfile, 'error num: ', ierr + goto 102 + 1500 print *, 'Error opening input: ', infile, 'error num: ', ierr + goto 102 + 2000 print *, 'Error closing input: ', infile, 'error num: ', ierr + goto 102 + 2500 print *, 'Error closing output: ', outfile, 'error num: ', ierr + goto 102 + 102 continue + + end http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/makefile ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/Fortran/makefile b/climatology/clim/orig/Fortran/makefile new file mode 100644 index 0000000..a624096 --- /dev/null +++ b/climatology/clim/orig/Fortran/makefile @@ -0,0 +1,46 @@ +# makefile for guassinterp. +# choose f77 or f90 compiler +# and appropriate flags + +# f77 compiler +FC = f77 +FFLAGS = -g +#FFLAGS = -O -n32 +DEFINE = -DLANG_F77 +HDF = /usr/local/hdf +#HDF = /usr/local/hdf4.1r3-n32 + +# f90 compiler +#FC = f90 +#FFLAGS = -n32 -bytereclen -cpp -extend_source +#DEFINE = -DLANG_F90 + +PGM = gaussinterp +LIBS = -L$(HDF)/lib -ldf -ljpeg -lz -lgen +INCLUDE = -Wf,-I/usr/local/include/hdf +BIN = /sst/vol7/PATHFINDER/bin/sgi + +SRC = setupinterp.f interp.f binsum.f getbase.c passbase.f +OBJ = setupinterp.o interp.o binsum.o getbase.o passbase.o + +$(PGM): $(OBJ) + $(FC) $(FFLAGS) -o $(PGM) $(OBJ) $(LIBS) + +setupinterp.o: + $(FC) $(FFLAGS) $(DEFINE) -c setupinterp.f +interp.o: +# $(FC) $(DEFINE) -c interp.f + $(FC) $(FFLAGS) $(DEFINE) -c interp.f +binsum.o: + $(FC) $(FFLAGS) $(DEFINE) -c binsum.f +passbase.o: + $(FC) $(FFLAGS) $(DEFINE) -c passbase.f +getbase.o: + $(CC) $(FFLAGS) -c getbase.c + +install: $(PGM) + cp $(PGM) $(BIN) + +clean: + rm -f $(OBJ) $(PGM) + http://git-wip-us.apache.org/repos/asf/incubator-sdap-nexus/blob/ff98fa34/climatology/clim/orig/Fortran/passbase.f ---------------------------------------------------------------------- diff --git a/climatology/clim/orig/Fortran/passbase.f b/climatology/clim/orig/Fortran/passbase.f new file mode 100644 index 0000000..79e2d8e --- /dev/null +++ b/climatology/clim/orig/Fortran/passbase.f @@ -0,0 +1,9 @@ + subroutine passbase( filename ) + + character*25 filename + character*25 basename + common /basefile/ basename + + basename = filename + end +
