; ============================================================= ; ; Example of reading ASCII data files of the SO2 column densities ; as produced by TEMIS/PROMOTE using IDL. ; ; For information, see the on-line product information page ; http://www.oma.be/BIRA-IASB/Molecules/SO2archive/info/?kind=dataspecidlasc ; ; (c) TEMIS / PROMOTE / BIRA-IASB ; last modified: 17 August 2007 ; ; ============================================================= ; ; Usage within IDL to compile ; ; IDL> .run so2ascread ; ; And to run: ; ; IDL> so2ascread ; ; after which the routine will ask for the name of a configuration file ; which holds the names of the ASCII files to process. ; The configuratile has thefollowing format: ; ; ; ; ; ; .... ; ; ; where the files to read are "/", that is: ; the string is prepended to the strings. ; ; ==================================================================== pro so2ascread ; Ask for the number of data files to cycle over. configfile='' print,'' print,'>>> Give the name of configuration file. read,configfile ;configfile='so2ascread.in' configfile=strtrim(configfile,2) ; Some definitions needed later on lat=fltarr(5) lon=fltarr(5) line="" ; Check if the configuration file exists. fd_id=FILE_TEST(configfile,/READ) if fd_id ne 1 then begin print,' *** Error: config file does not exist' goto,endprog endif ; Read the configuration file. openr,1,configfile,ERROR=err if err ne 0 then begin print,"*** Error opening config file ",configfile," -- ",!ERROR_STATE.MSG goto,endprog endif datapath='' readf,1,datapath datapath=strtrim(datapath,2) nfiles=-1 readf,1,nfiles if nfiles le 0 then begin print,"*** No files to read" close,1 goto,endprog endif datafiles=strarr(nfiles) for nf=0,nfiles-1 do begin readf,1,line line=strtrim(line,2) datafiles(nf)=line endfor close,1 ; Use as colour table "Rainbow+white". ; The white is used for the background of the plot. ; The first colour of this table is black, which makes it not possible ; to use black as colour for any other line, such as the coast lines. ; Therefore it is better to start higher up in the colour table, e.g. ; at index 32 (some shade of blue) and go up to 254 (red). device, decomposed=0, retain=2 loadct, 39, /silent ; Settings for the map to plot: ; - centre longitute and latitude ; - size of the plot ; 1) For a 30x30 degree plot, corresponding to the ; "Ecuador / Colombia" region of the Volcanic ; SO2 service. POlat= 0.0 POlon=-77.5 POlim=[POlat-15.0,POlon-15.0,POlat+15.0,POlon+15.0] ; 2) For a world plot POlat=0.0 POlon=0.0 POlim=[POlat-90.0,POlon-180.0,POlat+90.0,POlon+180.0] ; No rotation Rot=0.0 ; The map. map_set, POlat, POlon, Rot, LIMIT=POlim, $ /cylindrical, /continents, con_color=2, $ e_horizon={fill:1, color:255} ; Loop over the files to process print,'' for nf=0,nfiles-1 do begin ; Read the name of the ASCII file to be read and remove ; leading and trailing spaces. print,'--- Reading file ',datafiles(nf) filename=datapath+datafiles(nf) ; Check whether the ASCII file exists. fd_id=FILE_TEST(filename,/READ) if fd_id ne 1 then begin print,' *** Error: ASCII file does not exist' goto,endprog endif ; Open the datafile openr,2,filename,ERROR=err if err ne 0 then begin print,"*** Error opening ",filename," -- ",!ERROR_STATE.MSG goto,endprog endif ; Loop over the lines in the file and read (part of) the data. ; Lines starting with '#' are comment lines and have to be skipped. ; The data lines could be read with the explicit data format given ; in the header of the data file. But in order to easily skip any ; comment line, the file is read below line by line into a string, ; and then the data is read from the string. while not eof(2) do begin ; Read a line from the data file as a string. readf,2,line ; Split the line read into its components, the data columns. DataArr=strsplit(line,count=ncol,/EXTRACT) ; Skip any comment line. if DataArr(0) eq "#" then goto,nextline ; There are forward scans and backscans, and this can ; be seen from the third entry on the data line. reads,DataArr(2),pixelid ; Backscan pixels should not be plotted and can thus ; be skipped while reading. if pixelid eq 3 then goto,nextline ; Read pixel corner latitudes. reads,DataArr(3),lat1 reads,DataArr(4),lat2 reads,DataArr(5),lat3 reads,DataArr(6),lat4 lat(0)=lat1 lat(1)=lat3 lat(2)=lat4 lat(3)=lat2 lat(4)=lat1 ; Read pixel corner longitudes. reads,DataArr(8),lon1 reads,DataArr(9),lon2 reads,DataArr(10),lon3 reads,DataArr(11),lon4 lon(0)=lon1 lon(1)=lon3 lon(2)=lon4 lon(3)=lon2 lon(4)=lon1 ; Read the SO2 data, which is given in DU. ; The VCD read is that of the 6-km plume height, ; i.e. the one used in the plots on the website. reads,DataArr(16),so2scdval reads,DataArr(17),so2scderr reads,DataArr(27),so2vcdval reads,DataArr(28),so2vcderr ; Read the cloud cover data. reads,DataArr(37),ccival reads,DataArr(38),ccfval reads,DataArr(39),ctpval ; Plot the SO2 slant column quant = so2scdval ; Do not plot "no data" pixels if quant gt -90. then begin ; Plot values from 0.5 (blue) to 3 (red) DU icol=fix((quant-0.5)*(254-32)/2.5)+32 if icol lt 32 then icol=32 if icol gt 255 then icol=254 polyfill,lon,lat,color=icol endif ; Read the next data line. nextline: endwhile ; Next data file. close,2 endfor map_continents, /coasts, color = 0, mlinethick=2 ; End of the program. endprog: print,'' end ; ====================================================================