|
 |
Linuxskolen del 19 (Linuxmagasinet 3/2005)
I forrige artikkel viste vi hvordan du kan lage en enkel plotterutine, hvordan du kunne utvide denne til å oppfatte museklikk, samt hvordan denne rutinen kunne kjøres på kommandolinjen. Her ser vi på mer avansert plotting med IDL
Vi skal ta for oss en datafil fra ALOMARs værstasjon. Denne gjør vi først noen små forandringer i, som å konvertere datoformatet i hver linje til Juliansk kalender, samt sjekke dataene for uønskede datapunkter. Årsaken til at vi konverterer datoformatet er at man da blir stående med et langt flyttall (long float) som man kan regne med, og ikke en lang tekst streng (string) som er noe uhåndterlig i matematisk sammenheng. Filtreringen jeg nevnte er nødvendig for å fjerne datapunkter i datasettet da vindsensoren som benyttes av og til leverer transienter hvor vindhastigheten er uforholdsmessig stor i forhold til faktisk vind. Siden jeg benytter dynamisk skalering av plottene, vil slike transienter føre til at vindplottet får en y-skala som strekker seg opp til maksverdien for en slik transient, mens den faktiske vindhastigheten omtrent drukner i den nederste delen av skalaen.
multiplott
Det vi skal ende med er et plott som viser vind, temperatur, luftfuktighet og lufttrykk. Med andre ord skal vi produsere noe som i IDL-verdenen kalles for "Multiplot". Dette plottet skal også kjøres i bestemte intervaller for å vise utviklingen utover dagen.
La oss først se litt på selve datafilen. Eksemplet under viser bare en liten del av en slik fil som i realiteten inneholder én linje data for hvert sekund I løpet av 24 timer. Dette skulle bli 86400 linjer!
Eksempellinjen under viser hvordan en slik linje ser ut før jeg konverterer formatet på dato og klokkeslett:
dato | vind | snittvind | vindretning | snitt vindretning | temperatur | snitt-temp | trykk | snitt-trykk | sikt | snitt-sikt | solstråling | snitt solstråling | fuktighet | snitt-fuktighet | vindkjøleeffekt | snitt vindkjøleeffekt
2005-05-08 00:00:00 | 7.1 | 7.1 | 132.0 | 133.23 | 1.9 | 2.06 | 959.1 | 959.1 | 2997.9 | 2997.7| 28.7 | 31.3 | 78.2 | 77.89 | -2.82 | -8.78377578000098
Til å konvertere datofeltet og fjerne eventuelle feil fra vindsensoren benytter jeg to små programsnutter. Det ene er en Perl-rutine, laget av min gode venn og Perl-guru Martin Langteigen (langteigen@serenityofspace.com) som konverterer dato-delen til Juliansk dagnummer. Den andre lille snutten er et lite C-program, skrevet av en kollega av meg I Australia. Det sjekker filene for vindhastigheter over det som er rimelig å forvente i forhold til datasettet forøvrig. De verdiene som er over den oppgitte grensen erstattes med verdien til neste linje som forventes å være ok. Dette har fungert utmerket. Både konvertering og filter kjøres som en integrert del av selve IDL-programmet som jeg vil vise om litt. Etter en slik konvertering ser samme data-linje slik ut:
2453498.50000 | 7.1 | 7.1 | 132.0 | 133.23 | 1.9 | 2.06 | 959.1 | 959.1 | 2997.9 | 2997.7 | 28.7 | 31.3 | 78.2 | 77.89 | -2.82 | -8.78377578000098
Som du ser er nå første del av linjen blitt til et langt flyttall. Dette er samme dato som tidligere, men nå regnet ut som et Juliansk dagnummer fra år 0. For å lære mer om dette kan du starte her: http://no.wikipedia.org/wiki/Juliansk_kalender.
La oss starte med selve programmet nå. Det første vi gjør er å sette opp selve rammen rundt det:
Pro autoweather_houm_wtph_sub
; 1105-2005: Værplott I LinMag.
; Av Kolbjorn Bekkelund
device, true=24, decomposed=0, retain=2
loadct, 5
end
Jeg forklarte forrige gang hva device-linjen gjør, så det vil jeg ikke ta opp denne gangen. Det neste vi må gjøre er å fortelle IDL hvor de ønskede datafilene kan hentes fra, definere et par variabler, konvertere dato og filtrere bort uønskede vinder, samt hente systemklokken. Da ser det slik ut:
Pro autoweather_houm_wtph_sub
; 1105-2005: Værplott I LinMag.
; Av Kolbjorn Bekkelund
device, true=24, decomposed=0, retain=2
loadct, 5
lastfile = ''
fil_dir = ''
spawn, '/data/weather/convert'
fil_dir = '/data/weather/converted/*'
filnam = FILE_SEARCH(fil_dir,count=iff)
dato= systime()
end
lastfile og fildir defineres og settes til "ingen" verdi. Videre kjører jeg en kommando som heter "spawn". Denne er en utrolig hendig sak som lar meg kjøre Linux-kommandoer direkte I IDL! Dette benytter jeg meg av ofte siden det lar meg kombinere IDL og Linux på en meget smidig måte. "convert" er et bash-script som jeg har laget for å kombinere de to nevnte filtrene, samt lage meg en ny datafil når disse er kjørt. "convert" ser slik ut:
#!/bin/bash
cd /data/weather/online/
rm -rf /data/weather/converted/latest.log
rm -rf /data/weather/converted/substitute.dat
a=`ls -1rt 2005*.log | tail -1`
echo $a
# remove error winds
/data/weather/substitute -t 40 $a /data/weather/converted/substitute.dat
# convert date and time to julian day
/usr/bin/perl /data/weather/cnvrt.pl /data/weather/converted/substitute.dat > /data/weather/converted/$a
rm -rf /data/weather/converted/substitute.dat
echo $a
exit
Etter dette går jeg inn med IDL-rutinen FILE_SEARCH og henter den siste konverterte datafilen i angitte katalog. Det er denne vi skal plotte fra. For å kunne gjøre noe som helst med datasettet, må vi lese det inn i IDL. Dette gjøres ved at vi lager en to-dimensjonal array. Det er flere måter å gjøre dette på. Jeg bruker en metode hvor jeg utnytter en av IDLs egne rutiner for å lese inn ASCII-baserte data. Denne kalles "ASCII_TEMPLATE" og må først "læres opp" til å håndtere mitt datasett. I dette ligger det at jeg kan fortelle ASCII_TEMPLATE hvordan jeg vil at den skal lese inn dataene. Om den skal ta vare på eventuelle "headere", hvordan den skal sortere data (kommaseparerte, tab, kolon etc), samt hvor mange linjer den skal ta med. Dette vil jeg ikke demonstrere her, men bare si at det jeg ender ut med er en mal for hvordan mitt datasett skal leses inn, og denne malen henviser jeg til i koden som "weather_template_new.dat". Programmet vårt ser nå slik ut:
Pro autoweather_houm_wtph_sub
; 1105-2005: Værplott I LinMag.
; Av Kolbjorn Bekkelund
device, true=24, decomposed=0, retain=2
loadct, 5
lastfile = ''
fil_dir = ''
spawn, '/data/weather/convert'
fil_dir = '/data/weather/converted/*'
filnam = FILE_SEARCH(fil_dir,count=iff)
dato= systime()
FileTable=filnam[iff-1] ; get the last file in the online-dir
shortdate = strmid((Filetable), 24, 8)
Restore, '/data/weather/eval/weather_template_new.dat'
data=read_ascii(FileTable, template = weather_template_new)
data=data.field01
end
Vi har nå en array på 86400 x 14 celler. Med andre ord,en ganske stor sak!
Det neste som må gjøres er å sette opp selve grunnplottet. Det gjøres slik:
!Y.OMARGIN=[4, 1]
!X.OMARGIN = [2,7]
!p.multi = [0,1,4]
!p.font = -1
;!p.charsize=2.4
window, 0, xsize=800, ysize=800, title = 'ALOMAR Observatory - Wind, temperature, humidity and pressure since midnight'
Her forteller jeg IDL at jeg ønsker å h litt "luft" på sidene i det ferdige plottet. Det gjøres med !Y.OMARGIN etc. Videre sier jeg at jeg ønsker å plotte 4 grafer i samme plott ved hjelp av funksjonen !p.multi. Jeg sier da at jeg ønsker alt i en kolonne med fire rader. Grafene blir m.a.o. stående over hverandre. Jeg setter også opp fontstørrelse og type, slik at det ser ok ut. Tilslutt definerer jeg selve vinduet som blir til det bildet vi ender ut med når alt er ferdig. Vinduet settes til 800 x 800 pixler i størrelse. Da slipper man scrolling for å se hele plottet selv på relativt små skjermer. Tittelen øverst på plottet kommer fra parameret "title". Programmet vokser og ser nå slik ut:
Pro autoweather_houm_wtph_sub
; 1105-2005: Værplott I LinMag.
; Av Kolbjorn Bekkelund
device, true=24, decomposed=0, retain=2
loadct, 5
lastfile = ''
fil_dir = ''
spawn, '/data/weather/convert'
fil_dir = '/data/weather/converted/*'
filnam = FILE_SEARCH(fil_dir,count=iff)
dato= systime()
FileTable=filnam[iff-1] ; get the last file in the online-dir
shortdate = strmid((Filetable), 24, 8)
Restore, '/data/weather/eval/weather_template_new.dat'
data=read_ascii(FileTable, template = weather_template_new)
data=data.field01
!Y.OMARGIN=[4, 1]
!X.OMARGIN = [2,7]
!p.multi = [0,1,4]
!p.font = -1
;!p.charsize=2.4
window, 0, xsize=800, ysize=800, title = 'ALOMAR Observatory - Wind, temperature, humidity and pressure since midnight'
end
Det neste jeg vil gjøre er å kalkulere antall målinger (antall linjer) i datafilen. Dette kan jeg blant annet nytte til å finne ut tidspunktet for siste måling. Slik kan jeg følge med under plottingen om siste måling og plottetidspunkt skiller for mye. Siden jeg plotter fra datafilen mens den faktisk akkumuleres på serveren, kan jeg derfor oppdage på et tidlig stadium om værstasjonen feiler siden tiden for siste data blir for gammel I forhold til forventet. Dette gjør jeg slik:
number_samples = N_Elements(data(0,*))
data_set_time = strmid(strjoin(jul2cal(data(0,(number_samples-1)), /TO_STRING, /MDY)), 10, 8)
Jeg går her inn i den første kolonnen i datafilen (data(0) , og henter ut innføringen på den siste linjen med "number_samples-1". Denne verdien formatterer jeg så til normal dato og klokkeslett med "jul2kal", komprimerer bort eventuelle mellomrom med "strjoin" og til slutt henter jeg ut bare tidspunktet med "strmid" fra posisjon 10 og 8 karakterer utover. Lett som bare det! Det neste jeg skal vise et eksempel på er hvordan jeg henter ut en del tilleggsinformasjon til plottene mine. Det er jo ikke bare ønskelig å vise noen fine kurver, men det bør jo også legges inn en del tekstinfo som for eksempel hvor høy maks vindhastighet har vært, maks temperatur, minimum temperatur og så videre. Dette kalkulerer jeg for alle grafene før jeg starter selve plottingen. Slik ser det ut for vind-delen:
maxgust = max(data(1,*), maxindex)
maxgust_time = where(data(1,*) EQ maxgust); find the time of max gust in julian
number_maxgust_time = n_elements(maxgust_time); find number of same value
for i = 0, number_maxgust_time-1 do begin
maxgust_time_max = maxgust_time(number_maxgust_time-1)
endfor
maxgust_time_first = jul2cal(data(0,maxgust_time(0)), /TO_STRING, /MDY)
maxgust_time_latest = strmid(strjoin(jul2cal(data(0,maxgust_time_max), /TO_STRING, /MDY)), 10, 8); calculate date/time at max windgust
maxmean = max(data(2,*), maxindex)
maxmean_time = where(data(2,*) EQ maxmean); find the time of max gust in julian
number_maxmean_time = n_elements(maxmean_time); find number of same value
for i = 0, number_maxmean_time-1 do begin
maxmean_time_max = maxmean_time(number_maxmean_time-1)
endfor
maxmean_time_first = jul2cal(data(0,maxmean_time(0)), /TO_STRING, /MDY)
maxmean_time_latest = strmid(strjoin(jul2cal(data(0,maxmean_time_max), /TO_STRING, /MDY)), 10, 8); calculate date/time at max mean wind
Jeg velger her ut de verdiene jeg er på jakt etter ved hjelp av standardrutiner I IDL, og dette gjør jeg fra kolonne 1 og 2 I arrayen. Med andre ord fra faktisk kolonne 2 og 3 siden den første heter 0 og inneholder juliansk dagnummer som forklart tidligere. Etter dette står jeg med all nødvendig info om vind. Det samme gjør jeg for tempereatur, luftfuktighet og lufttrykk, bare at her kalkulrerer jeg også minimum-verdiene I tillegg. Minimum vind er ikke interessant siden den ikke kan bli lavere enn 0 m/s.
Det neste som må regnes ut og settes opp er tidsskalen for grafene. Dette gjøres slik:
trange=jd2js([MIN(data[0,*]),MAX(data[0,*])])
jstime=reform(jd2js(data[0,*]))
setup_time_axis, trange, 0, $
xticks=xticks, xminor=xminor, xtickv=xtickv,xrange=xrange, $
xtickname=xtickname, $
XSTYLE = xstyle, TIMEFORMAT = 'HH:MM'
Her foretar jeg en del krumspring. Det første jeg gjør er å hente ut max og min fra kolonnen med julianske verdier. Disse tar jeg så og konverterer til julianske sekunder før dette går inn i oppsettet for akser og "tickmarks". Jeg setter også opp at tidsskalaen skal benytte timer og minutter (HH:MM). Nå er vi faktisk klare til å plotte en graf. Jeg tar her for meg den delen av programmet som omhandler plotting av vinddata. Det ser slik ut:
yrange_wind_gust_max = Max(data(1, *))
PLOT, jstime, reform(data[1, *]), $
xticks=xticks,xminor=xminor,xtickv=xtickv,$
xrange=xrange,xtickname=xtickname,$
ystyle=8, background=255, color=0, yrange=yrange_wind_gust, charsize=2.4, $
YTITLE = 'Wind Speed (m/s)', yticklen=1, ygridstyle=1,/nodata, $
XTITLE = 'Local temp, humidity, pressure and wind @ '+ data_set_time +', plotted: '+ dato +''
OPLOT, jstime, reform(data[1, *]), color=105
OPLOT, jstime, reform(data[2, *]), color=37
axis, yaxis=1, color=0, charsize=2.4
Hva gjør jeg så her da? Vel, jeg starter med å definere y-skalaen for grafen. Denne settes til den høyeste verdien av alle innføringene i kolonne 1. Videre benytter jeg rutinen PLOT til å utføre selve plottingen. Her plotter jeg vinden fra kolonne 1 i arrayen som en funksjon av jstime som setter opp tidsskalaen i henhold til det vi fant ut for litt siden. Jeg setter også opp type tickmark, fontstørrelse, bakgrunns- og forgrunnsfarge, rutemønster og tittel på denne grafen. Dere ser også at jeg på slutten legger til OPLOT et par ganger. Den første plotter selve grafen for maksvind på ny for å få den i ønsket farge, mens OPLOT nummer to legger en kurve for middelvind på samme graf. Den siste innføringen - "axis", setter opp en ekstra y-akse på høyre side av grafen slik at avlesningen blir enklere. Nå begynner det å ligne på noe! Imidlertid mangler det ennå noe informativ tekst i plottene som maks vind, middelvind etc. Dette legger vi i samme plott slik:
xyouts, 100, 191, 'Max gust : ', color=105, charsize=1.2, /device
xyouts, 170, 191, color=0, decimals(maxgust, 2), charsize=1.2, /device
xyouts, 250, 191, color=0, maxgust_time_latest, charsize=1.2, /device
xyouts, 365, 191, 'Max mean: ', color=37, charsize=1.2, /device
xyouts, 443, 191, color=0, decimals(maxmean, 2), charsize=1.2, /device
xyouts, 520, 191, color=0, maxmean_time_latest, charsize=1.2, /device
Tallene etter xyouts forteller i antall pixler (y,x) hvor skriften skal plasseres i det totale plottet (ikke bare innenfor denne ene grafen!), mens color og charsize skulle være selvforklarende. Det som imidlertid ikke er så selvforklarende er "decimals". Dette er en IDL-rutine som formatterer verdien av "maxgust" (verdien av maks vind fra array) til et fornuftig kommatall. Uten "decimals" hadde jeg fått mange desimaler etter komma i de enkelte leddene og det ønsker jeg ikke. "/device" forteller IDL at teksten skal legges til plottet og ikke sendes til noen andre enheter som for eksempel skrivere etc. Uten /device dukker mao ikke teksten din opp i plottet!
Da gjenstår det bare å sette opp den delen av programmet som faktisk lar deg lage et bilde av plottet ditt. Dette kan være en png-fil eller andre formater din versjon av IDL støtter eller er lisensiert for. Jeg benytter meg av PNG siden jeg skal bruke bildene på web. Kommandoene for dette ser slik ut:
write_png, '/data/weather/plots/old-plots/'+ shortdate +'-wtph.png', tvrd(/true)
write_png, '/data/weather/plots/last-multiplot_wtph.png', tvrd(/true); save last plot for browsing on web
Her har jeg to linjer hvor den øverste lager et bilde som går inn i en samling av slike. Ett for hvert døgn. Den benytter variabelen "shortdate" som vi lagde helt i starten av programmet for å skille de enkelte dagene fra hverandre. Bildet plottes hver hver kjøring av programmet, men får KUN nytt navn en gang i døgnet. Linje to lager det bildet som legges ut på web. Dette bildet har alltid samme navn og blir dermed overskrevet hver gang.
Programmet ser nå slik ut i en meget forkortet versjon:
Pro autoweather_houm_wtph_sub
; 1105-2005: Værplott I LinMag.
; Av Kolbjorn Bekkelund
device, true=24, decomposed=0, retain=2
loadct, 5
lastfile = ''
fil_dir = ''
spawn, '/data/weather/convert'
fil_dir = '/data/weather/converted/*'
filnam = FILE_SEARCH(fil_dir,count=iff)
dato= systime()
FileTable=filnam[iff-1] ; get the last file in the online-dir
shortdate = strmid((Filetable), 24, 8)
Restore, '/data/weather/eval/weather_template_new.dat'
data=read_ascii(FileTable, template = weather_template_new)
data=data.field01
!Y.OMARGIN=[4, 1]
!X.OMARGIN = [2,7]
!p.multi = [0,1,4]
!p.font = -1
;!p.charsize=2.4
window, 0, xsize=800, ysize=800, title = 'ALOMAR Observatory - Wind, temperature, humidity and pressure since midnight'
number_samples = N_Elements(data(0,*))
data_set_time = strmid(strjoin(jul2cal(data(0,(number_samples-1)), /TO_STRING, /MDY)), 10, 8)
maxgust = max(data(1,*), maxindex)
maxgust_time = where(data(1,*) EQ maxgust); find the time of max gust in julian
number_maxgust_time = n_elements(maxgust_time); find number of same value
for i = 0, number_maxgust_time-1 do begin
maxgust_time_max = maxgust_time(number_maxgust_time-1)
endfor
maxgust_time_first = jul2cal(data(0,maxgust_time(0)), /TO_STRING, /MDY)
maxgust_time_latest = strmid(strjoin(jul2cal(data(0,maxgust_time_max), /TO_STRING, /MDY)), 10, 8) ; calculate date/time at max windgust
maxmean = max(data(2,*), maxindex)
maxmean_time = where(data(2,*) EQ maxmean); find the time of max gust in julian
number_maxmean_time = n_elements(maxmean_time); find number of same value
for i = 0, number_maxmean_time-1 do begin
maxmean_time_max = maxmean_time(number_maxmean_time-1)
endfor
maxmean_time_first = jul2cal(data(0,maxmean_time(0)), /TO_STRING, /MDY)
maxmean_time_latest = strmid(strjoin(jul2cal(data(0,maxmean_time_max), /TO_STRING, /MDY)), 10, 8) ; calculate date/time at max mean wind
trange=jd2js([MIN(data[0,*]),MAX(data[0,*])])
jstime=reform(jd2js(data[0,*]))
setup_time_axis, trange, 0, $
xticks=xticks, xminor=xminor, xtickv=xtickv,xrange=xrange, $
xtickname=xtickname, $
XSTYLE = xstyle, TIMEFORMAT = 'HH:MM'
yrange_wind_gust_max = Max(data(1, *))
PLOT, jstime, reform(data[1, *]), $
xticks=xticks,xminor=xminor,xtickv=xtickv,$
xrange=xrange,xtickname=xtickname,$
ystyle=8, background=255, color=0, yrange=yrange_wind_gust, charsize=2.4, $
YTITLE = 'Wind Speed (m/s)', yticklen=1, ygridstyle=1,/nodata, $
XTITLE = 'Local temp, humidity, pressure and wind @ '+ data_set_time +', plotted: '+ dato +''
OPLOT, jstime, reform(data[1, *]), color=105
OPLOT, jstime, reform(data[2, *]), color=37
axis, yaxis=1, color=0, charsize=2.4
xyouts, 100, 191, 'Max gust : ', color=105, charsize=1.2, /device
xyouts, 170, 191, color=0, decimals(maxgust, 2), charsize=1.2, /device
xyouts, 250, 191, color=0, maxgust_time_latest, charsize=1.2, /device
xyouts, 365, 191, 'Max mean: ', color=37, charsize=1.2, /device
xyouts, 443, 191, color=0, decimals(maxmean, 2), charsize=1.2, /device
xyouts, 520, 191, color=0, maxmean_time_latest, charsize=1.2, /device
write_png, '/data/weather/plots/old-plots/'+ shortdate +'-wtph.png', tvrd(/true)
write_png, '/data/weather/plots/last-multiplot_wtph.png', tvrd(/true); save last plot for browsing on web
end
Det siste vi skal se på er hvordan dette blir automatisert. For å gjøre dette kreves det tre steg:
1. Lage en hovedrutine I IDL som så kjører den sub-rutinen vi har laget.
2. Lage et lite bash-script som så kaller opp IDL med hovedrutinen som parameter. Dette bash-scriptet kjøres av cron.
3. Forandre window, 0, xsize=800, ysize=800 til window, 0, /pixmap, xsize=800, ysize=800
Dette hindrer at det hver gang IDL kjører, dukker opp et plott på skjermen som så forsvinner igjen straks plottingen er ferdig. Om noen bruker den maskinen dette kjøres på vil det virke forferdelig irriternede! /pixmap skjuler hele prosessen!
Hovedrutinen ser slik ut:
Pro autoweather_accumulated
autoweather_houm_wtph_sub
autoweather_houm_sv_sub ; denne er for andre deler av værstasjonen ikke vist her
exit, status=45
end
Bash-scriptet som så kaller på denne rutinen ser slik ut:
#!/bin/sh
#
# Run IDL programs to produce online weather plots on the ozone pc
# By Kolbjorn Bekkelund
#
source /home/ozon/.bashrc
cd /data/weather/eval/
/usr/local/bin/idl autoweather_accumulated.pro
exit
Nå mangler bare crontab-innføringen for dette:
*/5 * * * * /usr/local/bin/autorun_idl_weather
Vips! Været plottes av IDL hvert 5. minutt uten tilsyn!
Dette runder av IDL for denne gang. Det er mulig at jeg kommer tilbake med mindre tips og triks, men jeg har nå fått vist dere litt av det IDL og Linux kan benyttes til. De som ønsker å se hele programmet, kan sende meg en email. Da legger jeg også ved de IDL-rutinene som ikke følger med selve IDL, men som jeg har benyttet for å lage plottene.
|
 |
|
