pro fcwx ;declare procedure close, /all ;close files if any are open ;DECLARE DEFAULTS ;you can place this in a startup file so that you ;don't have to do this every time you run IDL device,true_color=24 device, retain=2, decomposed=0, set_character_size=[6,9] device, bypass_translation =0 device, get_visual_depth=depth print, 'Display depth: ', depth print, 'Color table size: ',!d.table_size loadct, 39 ; 0 = black, 255=white, etc !p.charsize=1.5 !p.color = 0 ; set draw color to black !p.background = 255 ; set background color to white !x.style = 1 !y.style = 1 !p.font = -1 !p.charsize = 1.3 !p.charthick = 1.5 ;------------------------------------ ;BEGIN BY READING IN DATA FROM FILE ;------------------------------------ ;The data is arranged in the following way ; ;stationID LAT LON elevation staNAME ;year month day maxT minT precip snowfall snowdepth ;year month day maxT minT precip snowfall snowdepth ;year month day maxT minT precip snowfall snowdepth ;... ... ... ... ... ... ... ... ;------------------------------------ ;THERE ARE MANY WAYS TO READ IN ASCII DATA. ;FIRST A LONG WAY AND THEN A SHORT WAY ;the long way will be used to show programming techniques ;the short way to show how to take advantage of IDL's routines ;short versus long toggle long=0 ;long=0 for short; long=1 for long ;declare file directory='/Volumes/atmos-www/dept/classes/programmingCourse/' filein=directory+'data/fort_collins_wx.txt' ;declare a string with file name ;------------------------------------ ; LONG WAY ;------------------------------------ IF (long) then begin maxlines=1000 ;declare maximum number of lines ;This will be cut-down later ;begin declaring our reading variables as integers ;IDL is NOT case sensitive header=' ' ;read the first line as a header (don't need those numbers) _year=0 ;dummy _month=0 ;dummy _day=0 ;dummy _maxT=0 ;dummy _minT=0 ;dummy _precip=0 ;dummy _snow=0 ;dummy _snowD=0 ;dummy year=FLTARR(maxlines) ;generate an floating point arrays containing month=FLTARR(maxlines) ;the maximum number of elements (1000) day=FLTARR(maxlines) ;floating point arrays are used here to avoid maxT=FLTARR(maxlines) ;potential mathematical errors of rounding. minT=FLTARR(maxlines) ;this is more expensive, but safer for science (opinion) precip=FLTARR(maxlines) snow=FLTARR(maxlines) snowD=FLTARR(maxlines) ; openr, lun, filein ,/get_lun openr, 13, filein ;open file for reading ;The 13 refers to the logical unit number. ;This tells the computer which file to ;access. ;BEGIN READING THE FILE for i=0,maxlines-1 do begin ;loop from 0 to 9999 (IDL starts indexing with 0) ;this line throws a flag to IDL that we have an input/output error ;if this occurs, it will exit the loop and go to fileend: below on_ioerror, fileend ;if statement to check for i=0. if this is true, then ;IDL will read the header line from file 13. if (i eq 0) then begin readf, 13, header ;if i /= 0, then IDL will read the data for each line (i) endif else begin ;read data values from file 13 readf, 13, _year, _month, _day, _maxT, _minT, $ _precip, _snow, _snowD ;store data in the FLTARR's. Notice the indexing is i-1. ;This accounts for the i=0 statement year[i-1]=_year month[i-1]=_month day[i-1]=_day maxT[i-1]=_maxT minT[i-1]=_minT precip[i-1]=_precip snow[i-1]=_snow snowD[i-1]=_snowD endelse endfor ;goto statement for i/o error and close the file fileend: close, 13 ;free_lun, lun ;based on the i/o error, store the ival to know the number of lines ;in the file. nvals=i-1 ;i-1 to account for header ;shrink arrays from maxlines to nvals year=year[0:nvals-1] month=month[0:nvals-1] day=day[0:nvals-1] maxT=maxT[0:nvals-1] minT=minT[0:nvals-1] precip=precip[0:nvals-1] snow=snow[0:nvals-1] snowD=snowD[0:nvals-1] endif else begin ;------------------------------------ ; SHORT WAY ;------------------------------------ ;read all data into a structure {data} ;data_start=1 skips the header line (remember idl starts at 0) data=READ_ASCII(filein,data_start=1) ;grab fields year=data.field1[0,*] month=data.field1[1,*] day=data.field1[2,*] maxT=data.field1[3,*] minT=data.field1[4,*] precip=data.field1[5,*] snow=data.field1[6,*] snowD=data.field1[7,*] ;this will produce a [1,nvals] array ;lets also reform this to a [nvals] vector nvals=n_elements(year) year=reform(year,[nvals]) month=reform(month,[nvals]) day=reform(day,[nvals]) maxT=reform(maxT,[nvals]) minT=reform(minT,[nvals]) precip=reform(precip,[nvals]) snow=reform(snow,[nvals]) snowD=reform(snowD,[nvals]) endelse ;------------------------------------ ; END READING DATA ;------------------------------------ ;LETS DEAL WITH MISSING DATA (9999) using WHERE missval=9999 bad=where(year EQ missval,badcount) ;find WHERE year equals 9999 if (badcount gt 0) then year[bad]=!values.f_nan ;if we have MD, set to NaN bad=where(month EQ missval,badcount) if (badcount gt 0) then month[bad]=!values.f_nan bad=where(day EQ missval,badcount) if (badcount gt 0) then day[bad]=!values.f_nan bad=where(maxT EQ missval,badcount) if (badcount gt 0) then maxT[bad]=!values.f_nan bad=where(minT EQ missval,badcount) if (badcount gt 0) then minT[bad]=!values.f_nan bad=where(precip EQ missval,badcount) if (badcount gt 0) then precip[bad]=!values.f_nan ; if (bad ne -1) then precip[bad]=!values.f_nan bad=where(snow EQ missval,badcount) if (badcount gt 0) then snow[bad]=!values.f_nan bad=where(snowD EQ missval,badcount) if (badcount gt 0) then snowD[bad]=!values.f_nan ;DEAL WITH TRACE PRECIPITATION DATA (9998) AND STORE THEM traceval=9998 trace=where(precip EQ traceval,ntrace) if (ntrace gt 0) then precip[trace]=0. ;------------------------------------ ; BEGIN PLOTTING ; ; lets begin with a simple plotting exercise to see how temperature ; has varied over our 1-year of data. Lets make a plot of mean temp ; with the min and max denotes as "error" bars ;------------------------------------ window, xsize=1000, ysize=600 ;calculate mean temp meanT=(maxT+minT)/2. title='Daily Temperature' xtitle='Day' ytitle='Temperature (F)' ;get our temperature range and allow extra "room" ymin=min(minT) ymax=max(maxT) yrange=[ymin-5.,ymax+5.] ;plot mean temperature as a thick line plot, meanT, yrange=yrange, thick=5, charsize=2, $ title=title, ytitle=ytitle, xtitle=xtitle stop ;make max/min region to be filled for i=0,nvals-2 do begin ;note nvals-2 because of i+1 below ;account for missing if ((finite(meanT[i]) eq 1) and $ (finite(maxT[i]) eq 1) and $ (finite(meanT[i+1]) eq 1) and $ (finite(maxT[i+1]) eq 1)) then begin ;fill MAX region RED polyfill, [i,i,i+1,i+1], [maxT[i],meanT[i],maxT[i+1],meanT[i+1]], $ /data, color=250 ;fill MIN region BLUE polyfill, [i,i,i+1,i+1], [minT[i],meanT[i],minT[i+1],meanT[i+1]], $ /data, color=50 endif endfor ;overplot mean temperature for cleaner look oplot, meanT, thick=5 ;add legend ;this legend.pro is not standard with IDL ;I have downloaded this from the web and modified (see website) legend, ['Max Range','Min Range','Mean Temp'],colors=[250,50,0],linestyle=[0,0,0],$ thick=[1,1,5],position=[.13,.89],/normal,charsize=2 stop ;the stop allows you to keep the variables active from the command line end ;end procedure