pro quick_interp_tdm2,year1,year2,out_prefix,dist,$ log=log,gs=gs,area=area,check=check,angular=angular,nostn=nostn,$ synth_prefix=synth_prefix,limit=limit,hires=hires,pts_prefix=pts_prefix,$ anomfac=anomfac,synthfac=synthfac,binfac=binfac,actfac=actfac,$ dumpbin=dumpbin,dumpmon=dumpmon,dumpglo=dumpglo,info=info,ratio=ratio ; runs Idl trigrid interpolation of anomaly files plus synthetic anomaly grids ; first reads in anomaly file data ; the adds dummy gridpoints that are ; further than distance (dist) ; from any of the observed data ; TDM: the dummy grid points default to zero, but if the synth_prefix files are present in call, ; the synthetic data from these grids are read in and used instead if n_params() lt 1 then begin print,'quick_interp,year1,year2,out_prefix,dist,' print,' binfac=binfac,gs=gs,/area,/log,/check,' print,' /hires' return endif ; close,/all ; ; set up a few defaults if n_elements(dist) eq 0 then xkm=500 else xkm=dist ; correlation decay distance if n_elements(constrain) eq 0 then constrain=0 else constrain=1 ; obsolete if n_elements(log) eq 0 then log=0 else log=1 ; obsolete if n_elements(binfac) eq 0 then binfac=1.0 ; multi factor for output bin files to write them ; as binary integer files (e.g. T would have an value ; of 10 for binfac) if n_elements(actfac) eq 0 then actfac=1.0 ; multi factor to obtain actual (i.e. =*1) values, ; required for dumping to .mon and .glo if n_elements(synthfac) eq 0 then synthfac=1.0 ; multi factor to obtain synth file actual values if n_elements(anomfac) eq 0 then anomfac=1.0 ; multi factor to obtain anom file actual values if n_elements(pts_prefix) eq 0 then pts_prefix='../fort.' ; option to have an alternative input file prefix if keyword_set(angular) eq 0 then angular=0 else angular=1 ; use angular distance interpolation rather ; than triangulation (slower!) if angular eq 1 then print,'Angular distance weighting set' limit=glimit(/all) ; sets limit to global field if n_elements(gs) eq 1 then gs=[gs,gs] ; sets grid spacing for output grid, or if gs not set if n_elements(gs) eq 0 then gs=[0.5,0.5] ; sets gridspacing to 0.5 deg lat/lon if n_elements(out_prefix) eq 0 then begin ; set up output file prefix out_prefix=' ' read,out_prefix,prompt='Enter prefix for output grid file: ' out_prefix=strip(out_prefix) endif im1=1 & im2=12 monthname=strarr(12) monthname[*]=['01','02','03','04','05','06','07','08','09','10','11','12'] get_lun,LogLun openw,LogLun,'/cru/tyn1/f014/scratch/log-qit2.dat' if keyword_set(info) then printf,-1," MEAN AV MAG STD DEV",format="(18x,a30)" relaxval=0.0 & if keyword_set(ratio) then relaxval=1.0 ;START OF ANNUAL LOOP ; for iy=year1,year2 do begin grid=fltarr(360/gs(0),180/gs(1),12)-9999 ; set up output grid if keyword_set(synth_prefix) then begin synthfile=strtrim(synth_prefix,2)+strtrim(string(iy),2) rdbin,dummygrid,synthfile,gridsize=2.5,/quiet synthmiss=where(dummygrid eq -9999,nsynthmiss) if (nsynthmiss gt 0) then dummygrid(synthmiss)=relaxval dummygrid=dummygrid/synthfac endif else begin dummygrid=fltarr(144,72,12) ; set up a "zero anomaly" grid for infilling spaces with missing data ; missing data defined as areas of grid further than the decay distance from any ; real station point dummygrid(*,*,*)=relaxval endelse for im=im1,im2 do begin jm=im if n_elements(nostn) eq 0 then begin ptsfile=pts_prefix+string(iy,".",im,".txt",form='(i4.4,a1,i2.2,a4)') UnixCommand = "cat " + ptsfile + " | wc -l" spawn, UnixCommand, NCount NLine=fix(strtrim(NCount)) NStn=NLine(0) pts1=fltarr(NStn,5) openr, lunMon, ptsfile, /get_lun for XLine = 0, NStn-1 do begin readf, lunMon, lat,lon,elv,dat,wmo, format="(2f8.2,f8.1,f13.5,i7)" pts1(XLine,0)=lat & pts1(XLine,1)=lon & pts1(XLine,2)=(dat*anomfac) & pts1(XLine,3)=elv endfor free_lun, lunMon if keyword_set(info) then begin mean=mean(pts1(*,2)) & avmag=mean(abs(pts1(*,2))) & stddev=stddev(pts1(*,2)) printf,-1,"data ",iy," month ", im, mean, avmag, stddev,format="(a5,i4,a7,i2,3f10.4)" endif endif else begin pts1=fltarr(2,5) pts1(*,0)=90.0 & pts1(*,1)=0.0 & pts1(*,2)=relaxval & pts1(*,3)=0.0 endelse ;--------------------------------------------------- ; map all points with influence radii - that is with decay distance [IDL variable is decay] ; this is done in the virtual Z device, and essentially finds all points on a 2.5 deg grid that ; fall outside station influence dummymax=max(dummygrid(*,*,(im-1))) while dummymax gt -9999 do begin if keyword_set(test) eq 0 then begin set_plot,'Z' ; set plot window to "virtual" erase,255 ; clear current plot buffer, set backgroudn to white device,set_res=[144,72] endif else begin initx set_plot,'x' window,0,xsize=144,ysize=72 endelse lim=glimit(/all) nel=n_elements(pts1(*,0))-1 map_set,limit=lim,/noborder,/isotro,/advance,xmargin=[0,0],ymargin=[0,0],pos=[0,0,1,1] a=findgen(33)*!pi*2/32 for i=0.0,nel do begin x=cos(a)*(xkm/110.0)*(1.0/cos(!pi*pts1(i,0)/180.0))+pts1(i,1) x=x<179.9 & x=x>(-179.9) y=sin(a)*(xkm/110.0)+pts1(i,0) y=y>(-89.9) & y=y<89.9 catch,error_value ; avoids a bug in IDL that throws out an occasional ; plot error in virtual window if error_value ne 0 then begin error_value=0 i=i+1 goto,skip_poly endif polyfill,x,y,color=160 skip_poly: endfor catch,/cancel image=tvrd() ; read plot of missing data as grid and then (in lines below) ; find those areas that are white (i.e. value 255) and add ; these to an array of "missing" input points (pts2) set_plot,'X' pts2=fltarr(200000,5) pts2(0:nel,*)=pts1 ; add real input data to pts2 nlat=n_elements(image(0,*)) & dlat=180.0/nlat nlon=n_elements(image(*,0)) & dlon=360.0/nlon nelold=nel for lon=0,nlon-1 do begin ; add a selection of synth/missing points to pts2 for lat=0,nlat-1 do begin ; (cannot have regular grid as it makes tri unstable) if (image(lon,lat) ne 160 $ ; if new synth, use genuines only, else use all and (dummymax eq 0.0 or dummygrid(lon,lat,(im-1)) ne 0.0)) then begin nel=nel+1 pts2(nel,0:2)=[((lat*dlat)+(dlat/2.0)- 91+(2*randomu(seed))), $ ((lon*dlon)+(dlon/2.0)-181+(2*randomu(seed))),dummygrid(lon,lat,(im-1))] endif endfor endfor pts2=pts2(0:nel,*) if (dummymax gt relaxval) then begin dummygrid(*,*,(im-1))=relaxval ; on next iteration, fill with relaxation pts1=pts2 & pts2=0 endif else begin dummygrid(*,*,(im-1))=-9999 ; no more iterations pts1=0 endelse dummymax=max(dummygrid(*,*,(im-1))) endwhile n=where(pts2(*,2) gt -9999) case 1 of keyword_set(area): begin ; ANGULAR DISTANCE WEIGHTING bounds=[-180+gs(0),-90+gs(1),180-gs(0),90-gs(1)] r=area_grid(pts2(n,1),pts2(n,0),pts2(n,2),gs*2.0,bounds,dist,angular=angular) r=rebin(r,720,360) end else: begin if keyword_set(hires) eq 0 then begin ; LOW RESOLUTION TRIANGULATION bounds=[-180+gs(0)/2.0,-90+gs(1)/2.0,180-gs(0)/2.0,90-gs(1)/2.0] triangulate,pts2(n,1),pts2(n,0),tr,b r=trigrid(pts2(n,1),pts2(n,0),pts2(n,2),tr,gs,bounds,$ xgrid=xgrid,ygrid=ygrid) endif else begin ; HIGH RESOLUTION TRINAGULATION, THEN RESAMPLE TO LOW RES bounds=[-180+gs(0)/10.0,-90+gs(1)/10,180-gs(0)/10,90-gs(1)/10] triangulate,pts2(n,1),pts2(n,0),tr,b r=trigrid(pts2(n,1),pts2(n,0),pts2(n,2),tr,gs/5.0,bounds,$ xgrid=xgrid,ygrid=ygrid) if keyword_set(check) then help,r r=rebin(r,720,360) if keyword_set(check) then help,r endelse end endcase r=r<9999 & r=r>(-9998) grid(*,*,im-1)=r ; print,rnge(pts2(*,2)),rnge(r) dumpvalid=where(r NE 0,dumpvalidn) if keyword_set(dumpglo) then begin r=r*actfac SaveGrid=-1 if (gs[0] eq 0.5) then SaveGrid=12 if (gs[0] eq 2.5) then SaveGrid=22 if (gs[0] eq 5.0) then SaveGrid=23 if (SaveGrid eq -1) then print, "@@@@@ grid not found for save to .glo @@@@@" SaveFile=strip(out_prefix+monthname(im-1)+'.'+string(iy,form='(i4)')+'.glo') SaveTitle=strip(monthname(im-1)+string(iy,form='(i5)')) SaveGlo,SaveGrid,r,CallFile=Savefile,CallTitle=SaveTitle ; plot routine - change at will! if keyword_set(info) then begin mean=mean(r(dumpvalid)) & avmag=mean(abs(r(dumpvalid))) & stddev=stddev(r(dumpvalid)) printf,-1,"grid "," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,4x,a9,3f10.4,a,i9)" endif endif if keyword_set(check) then goto,l1 endfor l1: if keyword_set(smooth) then dist_smooth,grid,dist if keyword_set(dumpmon) then begin ; @@@@@@@@@ problem with save/load mon @@@@@@@@@ dumpvalid=where(grid NE 0,dumpvalidn) mean=mean(grid(dumpvalid)) & avmag=mean(abs(grid(dumpvalid))) & stddev=stddev(grid(dumpvalid)) printf,-1,"grid ",iy," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,i4,a,3f10.4,a,i9)" savefile=strip(out_prefix+string(iy,form='(i4)')+'.mon') ; print,savefile openw, lunSave, savefile, /get_lun for xmon = 0, 11 do begin r=grid(*,*,xmon) r=r*actfac for xlat = (180/gs(0))-1, 0, -1 do begin for xlon = (360/(gs(0)*2)), (360/gs(0))-1, 6 do begin xlon0=xlon+0 & xlon1=xlon+1 & xlon2=xlon+2 & xlon3=xlon+3 & xlon4=xlon+4 & xlon5=xlon+5 printf,lunSave,r(xlon0,xlat),r(xlon1,xlat),r(xlon2,xlat), $ r(xlon3,xlat),r(xlon4,xlat),r(xlon5,xlat),format='(6e12.4)' endfor for xlon = 0, (360/(gs(0)*2))-1, 6 do begin xlon0=xlon+0 & xlon1=xlon+1 & xlon2=xlon+2 & xlon3=xlon+3 & xlon4=xlon+4 & xlon5=xlon+5 printf,lunSave,r(xlon0,xlat),r(xlon1,xlat),r(xlon2,xlat), $ r(xlon3,xlat),r(xlon4,xlat),r(xlon5,xlat),format='(6e12.4)' endfor endfor endfor free_lun, lunSave endif if (keyword_set(dumpbin)) then begin dumpvalid=where(fix(grid*binfac) NE 0,dumpvalidn) mean=mean(grid(dumpvalid)) & avmag=mean(abs(grid(dumpvalid))) & stddev=stddev(grid(dumpvalid)) printf,-1,"grid ",iy," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,i4,a,3f10.4,a,i9)" flname=strip(out_prefix+string(iy,form='(i4)')) wrbin,grid*binfac,flname ; WRITE BINARY GRID (INTEGER), no compression endif if keyword_set(check) then goto,l2 endfor l2: close,LogLun free_lun,LogLun end