Commit 63f64722 authored by Ian Boyd's avatar Ian Boyd
Browse files

Upload New File

parent cb5de09f
PRO calc_var_si_conversion
;Procedure that can be called externally to perform checks on the VAR_UNITS value in the metadata and,
;based on the VAR_UNITS input, calculate and save to a file the VAR_SI_CONVERSION value.
; ----------
;Written by Ian Boyd - iboyd@astro.umass.edu
;
; History:
; 20061012: Introduced to IDLCR8HDF - Version 2.0
; (Previously these checks were performed in the Extract_and_Test routine)
; 20080302: var_unit_arr and unit_pre_arr added which hold, respectively, the valid VAR_UNITS and
; corresponding VAR_SI_CONVERSION values, and the set of UNIT_PREFIXs (previously
; values from the AVDC TAV file have been used). This has been done so that the
; routine can be stand-alone (i.e. not dependent on also having to read in a TAV file),
; and also because the original Envisat table.dat file does not contain the
; VAR_SI_CONVERSION values. The input parameters have been changed to reflect this: bu
; and up arrays (previously containing the VAR_UNIT and UNIT_PREFIX values from the TAV
; file), and nbu and nup (the number of elements in the bu and up arrays) are no longer
; used. The integer flag tab_type has been added to account for the different handling of
; some of the VAR_UNIT and VAR_SI_CONVERSION values between AVDC and original Envisat.
; The STOP_WITH_ERROR routine is no longer called in the event of an error, but an Error
; State string is returned instead. Bug fixed when testing for a UNIT_PREFIX - previously
; only the first character of the VAR_UNIT value was checked for a possible UNIT_PREFIX,
; thus 'deka' (da) was excluded. Bug Fixed when the calculated Power Value of the last
; base SI unit shown in VAR_SI_CONVERSION is '1', this is now set so that the power value
; does not appear - Version 3.0
; 20090611: Galileo added to var_unit_arr (AVDC); Last VAR_SI_CONVERSION value for ppmv, ppbv,
; pptv, and ppv changed to DIMENSIONLESS (was ppv) in var_unit_arr (AVDC);
; VAR_SI_CONVERSION values for molec changed to 0;1.66054E-24;mol (was 0;1;molec) in
; var_unit_arr (AVDC); VAR_SI_CONVERSION values for DU changed to 0;4.4615E-4;mol m-2
; (was 0;2.6867E20;molec m-2) in var_unit_arr; Stop power units being added to
; DIMENSIONLESS (e.g. for VAR_UNITS=ppmv2); Change EVDC VAR_SI_CONVERSION values in
; var_unit_arr to match Envisat Metadata Guidelines values - Version 3.01
; 20091117: EVDC VAR_SI_CONVERSION values set to the same as AVDC (only one var_unit_arr set);
; Fix bug which, in some cases, does not account for repeated units in VAR_UNITS when
; determining the VAR_SI_CONVERSION (e.g. VAR_UNITS=W m-2 sr-1 m-1); Put VAR_SI_CONVERSION
; units in power value order; Correctly scale multiple units by the power
; e.g. W2 = (m2 kg s-3)^2 (previously assumed only a single unit was being scaled);
; Generate error if third VAR_SI_CONVERSION value is 'DIMENSIONLESS' or 'NONE' but also
; includes extra values - Version 4.0b1
; 20101001: Set up for GEOMS compliance e.g. DIMENSIONLESS changed to 1; MJD2000 changed to MJD2K;
; NONE removed - Version 4.0b2
; 20101221: Bug fix - need to correctly account for when units cancel each other out. Was still
; assigned the value DIMENSIONLESS - Version 4.0b3
; 20110303: Bug fix - when a dimensionless unit includes a power value (e.g. ppmv2), the base unit
; in VAR_SI_CONVERSION stays as '1'; Add 'dB' to var_unit_arr - Version 4.0b4
; 20110719: Add 'pH' - Version 4.0b5
; 20200210: Create a version that can be called with arguments; Remove requirement for tab_type;
; This version renamed calc_var_si_conversion.pro/.sav/
;
; Inputs: vuvalue - a string containing the Metadata VAR_UNITS value to be checked (everything
; after the '=')
; rd - integer flag to indicate if VAR_DATA_TYPE is real/double (-1) or not (1). Used to
; make VAR_SI_CONVERSION values of the same type
;
; Output: var_si_conv - a string containing the VAR_SI_CONVERSION value determined by the program
; (saved to the file 'vsc_calc.txt')
; errstate - code returned to the program that called calc_var_si_conversion as follows:
; 0 = VAR_SI_CONVERSION calculated OK and written to 'vsc_calc.txt'
; 2 = Input arguments not valid (vsc_calc.txt not created)
;
var_unit_arr=['%;0;0.01;1','A;0;1;A','C;0;1;s A','cd;0;1;cd','d;0;86400;s','deg;0;1.74533E-2;rad',$
'degC;273.15;1;K','1;0;1;1','h;0;3600;s','Hz;0;1;s-1','J;0;1;m2 kg s-2','K;0;1;K',$
'l;0;1E-3;m3','lm;0;1;cd sr','lx;0;1;cd sr m-2','m;0;1;m','min;0;60;s',$
'MJD2K;0;86400;s','mol;0;1;mol','molec;0;1.66054E-24;mol','Np;0;1;1',$
'N;0;1E3;m kg s-2','Pa;0;1;m-1 kg s-2','pH;0;1E-12;m2 kg s-2 A-2','photons;0;1;photons',$
'ppbv;0;1E-9;1','ppmv;0;1E-6;1','pptv;0;1E-12;1','ppv;0;1;1','psu;0;1;psu',$
'rad;0;1;rad','s;0;1;s','sr;0;1;sr','V;0;1;m2 kg s-3 A-1','W;0;1;m2 kg s-3','kg;0;1;kg',$
'DU;0;4.4614E-4;mol m-2','Gal;0;1E-2;m s-2','dB;0;1;1']
unit_pre_arr=['Y;yotta;1E24','Z;zetta;1E21','E;exa;1E18','P;peta;1E15','T;tera;1E12','G;giga;1E9',$
'M;mega;1E6','k;kilo;1E3','h;hecto;1E2','da;deka;1E1','d;deci;1E-1','c;centi;1E-2',$
'm;milli;1E-3','u;micro;1E-6','n;nano;1E-9','p;pico;1E-12','f;femto;1E-15',$
'a;atto;1E-18','z;zepto;1E-21','y;yocto;1E-24']
;Check for valid input arguments
args = command_line_args(COUNT = iNrArgs)
var_si_conv='' & errstate=2 ;initialize outputs
if ( iNrArgs ne 2 ) then exit, status=errstate $
else if (STRTRIM(args(1),2) ne '-1') and (STRTRIM(args(1),2) ne '1') then exit, status=errstate $
else begin
vuvalue = STRTRIM(args(0), 2)
rd = FIX(args(1))
errstate=0
endelse
CD,CURRENT=curdir
curdir=curdir+path_sep()
nta=N_ELEMENTS(var_unit_arr)
;extract var_unit_arr/unit_pre_arr sub-values into vu/up arrays
vuhold=STRSPLIT(var_unit_arr[0],';',/Extract) ;test for number of sub-values
nvu=N_ELEMENTS(vuhold) & vu=STRARR(nvu,nta)
FOR j=0,nta-1 DO BEGIN
vuhold=STRSPLIT(var_unit_arr[j],';',/Extract)
vu[*,j]=vuhold
ENDFOR
nta=N_ELEMENTS(unit_pre_arr)
vuhold=STRSPLIT(unit_pre_arr[0],';',/Extract) ;test for number of sub-values
nup=N_ELEMENTS(vuhold) & up=STRARR(nup,nta)
FOR j=0,nta-1 DO BEGIN
vuhold=STRSPLIT(unit_pre_arr[j],';',/Extract)
up[*,j]=vuhold
ENDFOR
;separate out metadata sub-values into component parts, and set-up holding arrays
vp=STRSPLIT(STRTRIM(vuvalue,2),' ',/Extract) & vpn=N_ELEMENTS(vp)
vpx=STRARR(2,vpn) ;0 holds VAR_UNIT value, 1 holds POWER component
bpx=STRARR(vpn) ;holding string array for base unit
ex=INTARR(vpn)+1 ;holding integer array for power value (defaults to 1)
tm=DBLARR(vpn) ;holding array for scale factor
j=0
WHILE (j LE vpn-1) AND (errstate EQ '') DO BEGIN
;test to see if the sub-value is a base unit in the TAV file
ti=WHERE(vp[j] EQ vu[0,*],tcnt)
IF tcnt NE 0 THEN vpx[0,j]=vp[j] $ it is a base unit
ELSE BEGIN ;separate out into unit and power values as required
vpx[0,j]=STRMID(vp[j],0,1) ;save first character of vp(j)
stopchk=0
IF (STRLEN(vp[j]) GE 2) THEN BEGIN
FOR k=1,STRLEN(vp[j])-1 DO BEGIN
ah=STRMID(vp[j],k,1)
test1=(BYTE(ah) GE 65) AND (BYTE(ah) LE 90) ;A-Z
test2=(BYTE(ah) GE 97) AND (BYTE(ah) LE 122) ;a-z
IF (test1[0]) OR (test2[0]) THEN BEGIN
IF stopchk EQ 0 THEN vpx[0,j]=vpx[0,j]+ah ELSE stopchk=2
;if stopchk EQ 2 THEN this means that VAR_UNITS is not legal
ENDIF ELSE IF stopchk EQ 0 THEN BEGIN
stopchk=1 ;first non-alpha character so check for numeric or '-' character
test1=(BYTE(ah) EQ 45) OR ((BYTE(ah) GE 49) AND (BYTE(ah) LE 57)) ;- or 1-9
IF NOT test1[0] THEN stopchk=2 ELSE vpx[1,j]=ah
ENDIF ELSE IF stopchk EQ 1 THEN BEGIN
;need to check for a numeric character only
test1=(BYTE(ah) GE 48) AND (BYTE(ah) LE 57) ;0-9
IF NOT test1[0] THEN stopchk=2 ELSE vpx[1,j]=vpx[1,j]+ah
ENDIF
ENDFOR
IF vpx[1,j] NE '' THEN ex[j]=FIX(vpx[1,j])
ENDIF
IF stopchk EQ 2 THEN vpx[0,j]=vp[j] ;in the event that the value is not valid, so will create error
;Do TAV check on the VAR_UNIT value
ti=WHERE(vpx[0,j] EQ vu[0,*],tcnt)
ENDELSE
IF tcnt NE 0 THEN BEGIN ;VAR_UNIT is a BASE_VALUE
bemult=(DOUBLE(vu[nvu-2,ti[0]]))^ex[j] & tm[j]=bemult
bpx[j]=vu[nvu-3,ti[0]]+';'+STRTRIM(STRING(format='(E8.1)',bemult),2)+';'+vu[nvu-1,ti[0]]
ENDIF ELSE IF vpx[0,j] EQ 'g' THEN BEGIN ;check for VAR_UNIT EQ g (gram) for AVDC style TAV file
bemult=0.001d^ex[j] & tm[j]=bemult
bpx[j]='0;'+STRTRIM(STRING(format='(E8.1)',bemult),2)+';kg'
ENDIF ELSE BEGIN ;separate out vpx value into prefix and base-value and test
;check for valid prefix - first check for 'deka' (da)
pref=STRMID(vpx[0,j],0,2) & bas=STRMID(vpx[0,j],2)
pi=WHERE(pref EQ up[0,*],pcnt)
IF pcnt EQ 0 THEN BEGIN ;test for the remaining prefixes
pref=STRMID(vpx[0,j],0,1) & bas=STRMID(vpx[0,j],1)
pi=WHERE(pref EQ up[0,*],pcnt)
ENDIF
IF pcnt NE 0 THEN BEGIN
pmult=DOUBLE(up[nup-1,pi[0]])
;check for valid base
ti=WHERE(bas EQ vu[0,*],tcnt)
IF tcnt NE 0 THEN BEGIN
bpmult=(DOUBLE(vu[nvu-2,ti[0]])*pmult)^ex[j] & tm[j]=bpmult
bpx[j]=vu[nvu-3,ti[0]]+';'+STRTRIM(STRING(format='(E8.1)',bpmult),2)+';'+vu[nvu-1,ti[0]]
ENDIF ELSE IF bas EQ 'g' THEN BEGIN ;check for VAR_UNIT EQ g (gram) for AVDC style TAV file
bpmult=(pmult*0.001D)^ex[j] & tm[j]=bpmult
bpx[j]='0;'+STRTRIM(STRING(format='(E8.1)',bpmult),2)+';kg'
ENDIF ELSE errstate=2
ENDIF ELSE errstate=2
ENDELSE
j=j+1
ENDWHILE
IF errstate EQ 0 THEN BEGIN ;No errors detected so continue
;Create VAR_SI_CONVERSION value
tmult=1.0D
FOR j=0,vpn-1 DO tmult=tmult*tm[j]
;reformat the multiplier e.g. 1.000E+003 becomes 1E3
;convert to Exponential form if necessary
IF (tmult EQ 273.15D) OR (tmult MOD 60.D EQ 0.D) OR ((tmult GE 0.01D) AND (tmult LT 1.D2) $
AND (tmult*1.D4 MOD 1.D2 EQ 0.D)) THEN tmults=STRTRIM(STRUPCASE(tmult),2) $ ;keep the same format
ELSE tmults=STRTRIM(STRING(format='(E14.6)',tmult),2)
epos=STRPOS(tmults,'E') & ppos=STRPOS(tmults,'.')
IF epos NE -1 THEN BEGIN ;remove unnecessary characters after 'E'
ep=FIX(STRMID(tmults,epos+1)) & epx=STRTRIM(ep,2)
tmults=STRMID(tmults,0,epos+1)+epx
ENDIF
IF ppos NE -1 THEN BEGIN ;remove any trailing zeroes after the decimal place
IF epos NE -1 THEN ep=STRMID(tmults,ppos+1,epos-(ppos+1)) ELSE ep=STRMID(tmults,ppos+1)
WHILE STRMID(ep,STRLEN(ep)-1,1) EQ '0' DO ep=STRMID(ep,0,STRLEN(ep)-1)
IF ep NE '' THEN tmults=STRMID(tmults,0,ppos+1)+ep ELSE tmults=STRMID(tmults,0,ppos)
IF epos NE -1 THEN tmults=tmults+'E'+epx
ENDIF
;Scale the units by the power e.g. W2 = (m2 kg s-3)^2
FOR j=0,vpn-1 DO BEGIN
vsc=STRSPLIT(bpx[j],';',/EXTRACT)
vspl=STRSPLIT(vsc[2],' ',/EXTRACT,COUNT=vscnt)
sichk=STRARR(vscnt) & pwchk=sichk
IF (vpx[1,j] NE '') AND (vpx[1,j] NE '1') AND (vsc[2] NE '1') THEN BEGIN
bpx[j]=vsc[0]+';'+vsc[1]+';'
FOR k=0,vscnt-1 DO BEGIN
;separate out SI units and power values
sires=STRSPLIT(vspl[k],'-0123456789',/Extract)
sichk[k]=sires[0] ;SI Unit
IF STRLEN(sichk[k]) NE STRLEN(vspl[k]) THEN $
pwchk[k]=STRMID(vspl[k],STRLEN(sichk[k])) $
ELSE pwchk[k]='1' ;Power value
pwm=FIX(pwchk[k])*FIX(vpx[1,j])
IF k EQ 0 THEN sp='' ELSE sp=' '
bpx[j]=bpx[j]+sp+sichk[k]+STRTRIM(pwm,2)
ENDFOR
ENDIF
ENDFOR
;Put together VAR_SI_CONVERSION
vsc=STRSPLIT(bpx[0],';',/EXTRACT)
;IF vsc[2] EQ '1' THEN vpx[1,0]=''
var_si_conv=vsc[0]+';'+tmults+';'+vsc[2]
;add remaining base units to VAR_SI_CONVERSION
IF vpn GT 1 THEN $
FOR j=1,vpn-1 DO BEGIN
vsc=STRSPLIT(bpx[j],';',/Extract)
var_si_conv=var_si_conv+' '+vsc[2]
ENDFOR
;check VAR_SI_CONVERSION for repeated SI units e.g. m m-3 becomes m-2
schk=STRSPLIT(var_si_conv,' ;',/Extract) & scnt=N_ELEMENTS(schk)
IF scnt GT 3 THEN BEGIN ;more than one SI unit in VAR_SI_CONVERSION
sichk=STRARR(scnt-2) & pwchk=sichk
FOR j=0,scnt-3 DO BEGIN ;separate out SI units and power values
sires=STRSPLIT(schk[j+2],'-0123456789',/Extract)
sichk[j]=sires[0] ;SI Unit
IF STRLEN(sichk[j]) NE STRLEN(schk[j+2]) THEN $
pwchk[j]=STRMID(schk[j+2],STRLEN(sichk[j])) $
ELSE pwchk[j]='1' ;Power value
ENDFOR
j=0 & pwval=0
WHILE j LT scnt-3 DO BEGIN
si=WHERE((sichk[j] NE '') AND (sichk[j] EQ sichk[j+1:scnt-3]),sicnt)
IF sicnt NE 0 THEN BEGIN
FOR k=0,sicnt-1 DO BEGIN
si[k]=si[k]+j+1
pwval=FIX(pwchk[j])+FIX(pwchk[si[k]])
IF (pwval[0] EQ 0) AND (k EQ sicnt-1) THEN BEGIN
sichk[j]='' & pwchk[j]=''
ENDIF ELSE IF (pwval[0] EQ 1) and (k EQ sicnt-1) THEN pwchk[j]='' $
ELSE pwchk[j]=STRTRIM(pwval[0],2)
;make null all the other SI values
sichk[si[k]]='' & pwchk[si[k]]=''
ENDFOR
ENDIF ELSE IF pwchk[j] EQ '1' THEN pwchk[j]=''
j=j+1
ENDWHILE
;Also do check on the power value of the last SI Unit
IF pwchk[scnt-3] EQ '1' THEN pwchk[scnt-3]=''
;Put units in power value order
pwhold=pwchk
oi=WHERE(pwhold EQ '',ocnt)
IF ocnt NE 0 THEN pwhold[oi]='1'
pws=SORT(FIX(pwhold))
sichk=sichk[REVERSE(pws)] & pwchk=pwchk[REVERSE(pws)]
var_si_conv=schk[0]+';'+schk[1]+';'+sichk[0]+pwchk[0]
si=WHERE(sichk NE '',sicnt)
IF sicnt EQ 0 THEN var_si_conv=var_si_conv+'1' $ ;i.e. values cancelled out
ELSE BEGIN
FOR j=1,scnt-3 DO BEGIN
si=WHERE(sichk[0:j-1] NE '',sicnt)
IF (sichk[j] EQ '') OR (sicnt EQ 0) THEN var_si_conv=var_si_conv+sichk[j]+pwchk[j] $
ELSE var_si_conv=var_si_conv+' '+sichk[j]+pwchk[j]
ENDFOR
ENDELSE
ENDIF
IF rd LT 0 THEN BEGIN
;VAR_DATA_TYPE is Real or Double so make VAR_SI_CONVERSION values floating point
tchkh=STRSPLIT(var_si_conv,';',/Extract) & tup=STRUPCASE(tchkh)
FOR j=0,1 DO BEGIN
IF STRPOS(tup[j],'.') EQ -1 THEN BEGIN
IF STRPOS(tup[j],'E') EQ -1 THEN tchkh[j]=tchkh[j]+'.0' $
ELSE tchkh[j]=STRMID(tup[j],0,STRPOS(tup[j],'E'))+'.0'+STRMID(tup[j],STRPOS(tup[j],'E'))
ENDIF
ENDFOR
var_si_conv=tchkh[0]+';'+tchkh[1]+';'+tchkh[2]
ENDIF
;Check for invalid VAR_UNITS - third VAR_SI_CONVERSION is 1 plus extra values
vsc=STRSPLIT(var_si_conv,';',/EXTRACT)
vsc[2]=STRTRIM(vsc[2],2)
IF (errstate EQ 0) AND (STRMID(vsc[2],0,1) EQ '1') AND (STRLEN(vsc[2]) GT 1) THEN BEGIN
errstate=2 & var_si_conv=''
ENDIF
ENDIF
IF errstate EQ 0 THEN BEGIN
OPENW, ou, curdir+'vsc_calc.txt', /GET_LUN
PRINTF, ou, var_si_conv
FREE_LUN, ou
ENDIF
exit, status = errstate
END ;procedure calc_var_si_conversion
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment