datareq.pro

 InfoInfo
Search:    

;*******************************************************************************
function header

html = '<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"' $

return, html
end

;*******************************************************************************
function start_html, title=title, bgcolor=bgcolor

html = '<HEAD>' $

return, html
end

;*******************************************************************************
function end_html

html = '</BODY></HTML>' + string(10B)

return, html
end

;*******************************************************************************
function a, href=href, content=content

html = '<A HREF="' + href + '">' + content + '</A>' + string(10B)

return, html
end

;*******************************************************************************
function br

html = '<BR>' + string(10B)

return, html
end

;*******************************************************************************
function h1, align=align, content

if keyword_set(align) then begin

endif else begin

endelse

return, html
end

;*******************************************************************************
function h2, align=align, content

if keyword_set(align) then begin

endif else begin

endelse

return, html
end

;*******************************************************************************
function h3, align=align, content

if keyword_set(align) then begin

endif else begin

endelse

return, html
end

;*******************************************************************************
function h4, align=align, content

if keyword_set(align) then begin

endif else begin

endelse

return, html
end

;*******************************************************************************
function h5, align=align, content

if keyword_set(align) then begin

endif else begin

endelse

return, html
end

;*******************************************************************************
function table, align=align, border=border, cellpadding=cellpadding, width=width, $

html = '<TABLE ALIGN="' + align + '" BORDER="' + border + '" CELLPADDING="' + cellpadding + $

return, html
end

;*******************************************************************************
function tr, align=align, valign=valign, content=content

html = '<TR ALIGN="' + align + '" VALIGN="' + valign + '">' + string(10B) $

return, html
end

;*******************************************************************************
function th, content

html = '<TH>' + content + '</TH>' + string(10B)

return, html
end

;*******************************************************************************
function img, src=src, alt=alt, content=content, border=border

html = '<IMG SRC="' + src + '" ALT="' + alt + '" BORDER="' + border + '">' + string(10B)

return, html
end

;*******************************************************************************
pro doDataRequest, outputTypes, component, beginDateTime, endDateTime, stations, $

; globals
common HtmlFileCommons, gHtmlFile, gHtmlFileLUN
common OutputsCommons, gMakeBin, gMakeIAGA, gMakePNG, gMakePS, gBx, gBy, gBz

; ===== CONSTANTS

; kMaxInputFiles = total number of binary data files for all stations; if, for
; example, 9 stations are selected and two days of data are requested, then 18
; files would be needed...
kMaxInputFiles = 100

; ===== STRUCTURE DEFINITIONS

; StationPropStruct - station properties (name, IAGA code, etc.)
stationProp = {StationPropStruct, maccsName:, maccsCode:, iagaName:, iagaCode:, latitude:0.0D, longitude:0.0D}
; InputFileStruct - data file properties (name, path, associated station, Julian day number of file, etc)
inputFile = {InputFileStruct, exists:0, name:, pathPrefix:, urlPrefix:, jdn:0.0D, stationProp:stationProp}

; ===== OUTPUT TYPES

; convert comma-delimited outputTypes list to array
outputs = str_sep(outputTypes, ', ')

gMakeBin = 0
gMakeIAGA = 0
gMakePNG = 0
gMakePS = 0
foo = where(outputs, outputsCount)
if outputsCount ne 0 then begin

endif

; ===== COMPONENT

gBx = 0
gBy = 0
gBz = 0
case component of

endcase

; ===== STATIONS

; convert comma-delimited station code list to stations array
stations = str_sep(stations, ', ')
foo = where(stations, stationsCount)
if stationsCount ne 0 then begin

endif

; ===== DATE-TIMES

; prepare begining and ending julian day numbers
beginDT = str_sep(beginDateTime, ', ')
beginJDN = julday(beginDT[1], beginDT[2], beginDT[0], beginDT[3])
caldat, beginJDN, month, day, year, hour

endDT = str_sep(endDateTime, ', ')
endJDN = julday(endDT[1], endDT[2], endDT[0], endDT[3])
caldat, endJDN, month, day, year, hour

; ===== INPUT FILES

; prepare array of InputFileStruct structures
inputFilesArrayTemp = replicate(inputFile, kMaxInputFiles)
; fill the 'inputFilesArrayTemp' array
fillInputFileArray, inputFilesArrayTemp, inputPathPrefix, inputURLPrefix, $

; get 'inputFilesArrayTemp' array size
inputsCount = 0
for i = 0, kMaxInputFiles - 1 do begin

endfor
; resize the 'inputFilesArray' array to simplify results of sort()
if inputsCount gt 0 then begin

endif

; ===== HTML OUTPUT

gHtmlFile = htmlFilePath
openw, gHtmlFileLUN, gHtmlFile, /get_lun
printf, gHtmlFileLUN, header()
printf, gHtmlFileLUN, start_html(title="MACCS Data Request Results", bgcolor="#ffffff")
printf, gHtmlFileLUN, h3("MACCS Data Request Results", align = "center")

if stationFilesCount eq 0 then begin

endif else begin

endelse

printf, gHtmlFileLUN, end_html()
close, gHtmlFileLUN

end

;*******************************************************************************
pro fillInputFileArray, inputFilesArray, inputPathPrefix, inputURLPrefix, $

; construct data file names; here's a sample to illustrate:
;

————-path prefi————->
; <-url prefix ->
; /var/apache/htdocs/space/webdata/cd/92/cd92223.bin

; data file extension
inputFileExtension = '.bin'

sizeArray = size(inputFilesArray)
numInputFiles = sizeArray[1]
if numInputFiles eq 0 then goto, TheEnd
sizeStationProps = size(stationProps)
numStationProps = sizeStationProps[1]
if numStationProps eq 0 then goto, TheEnd

fileIndex = 0
for stationIndex = 0, numStationProps - 1 do begin

; non-inclusive of end julian day number

; kDirDelimiter = '\'

endfor

TheEnd:

end

;*******************************************************************************
pro processInputFiles, inputFilesArray, beginJDN, endJDN, $

common HtmlFileCommons
common OutputsCommons

sizeArray = size(inputFilesArray)
numInputFiles = sizeArray[1]

; HTML strings for each of the table cells
dummyString =
kMaxStationDays = 100
stationHtml = replicate(dummyString, kMaxStationDays)
binHtml = replicate(dummyString, kMaxStationDays)
iagaHtml = replicate(dummyString, kMaxStationDays)

; loop over input files; group stations together
kNumStations = 9
requestedStations = strarr(kNumStations)
availableStations = strarr(kNumStations)
availableStationIndex = -1
stationCode = 'InitialNonEmptyStationCode'
stationIndex = -1

firstFileForThisStation = 1
lastFileForThisStation = 0
fileNumberForThisStation = 0

sizeArray = size(inputFilesArray)
numInputFiles = sizeArray[1]

firstFileIndex = 0
lastFileIndex = numInputFiles - 1

for fileIndex = firstFileIndex, lastFileIndex do begin

    thisFile = inputFilesArray[fileIndex]

    thisMaccsCode = thisFile.stationProp.maccsCode

    thisMaccsName = thisFile.stationProp.maccsName

    fileNumberForThisStation = fileNumberForThisStation + 1

    ;; see if this file exists

    if thisFile.exists then begin

      thisFileExists = 1

    endif else begin

      thisFileExists = 0

    endelse

    ;; see if this is the last file for this station

    if fileIndex lt lastFileIndex then begin

      nextFile = inputFilesArray[fileIndex + 1]

      nextMaccsCode = nextFile.stationProp.maccsCode

      nextMaccsName = nextFile.stationProp.maccsName

      if nextMaccsCode ne thisMaccsCode then begin

        lastFileForThisStation = 1

      endif

    endif else begin

      lastFileForThisStation = 1

    endelse

    ;; see if this is the first file for this station

    if firstFileForThisStation then begin

      stationIndex = stationIndex + 1

      requestedStations[stationIndex] = thisMaccsCode

      stationHtml[stationIndex] = stationHtml[stationIndex] + h5(thisMaccsName)

      firstFileForThisStation = 0

    endif

    if thisFile.exists then begin

      thisStationIsAvailableAndItsNotInTheListOfAvailableStations = 1

      for index = 0, kNumStations - 1 do begin

        if availableStations[index] eq thisMaccsCode then begin

          thisStationIsAvailableAndItsNotInTheListOfAvailableStations = 0

        endif

      endfor

      if thisStationIsAvailableAndItsNotInTheListOfAvailableStations then begin

        availableStationIndex = availableStationIndex + 1

        availableStations[availableStationIndex] = thisMaccsCode

      endif

    endif

    jdn1 = thisFile.jdn

    caldat, jdn1, month1, day1, year1

    jdn2 = jdn1 + 1

    caldat, jdn2, month2, day2, year2

    ;; new file (same station or new station)

    if gMakeBin then begin

      fileName = thisFile.name

      ;; ammend html

      href = thisFile.urlPrefix + fileName

      if thisFile.exists then begin

        binHtml[stationIndex] = binHtml[stationIndex] + a(href=href, content=fileName)

        binHtml[stationIndex] = binHtml[stationIndex] + br()

      endif else begin

        binHtml[stationIndex] = binHtml[stationIndex] + fileName + ' (unavailable)'

        binHtml[stationIndex] = binHtml[stationIndex] + br()

      endelse

    endif

    if gMakeIAGA then begin

      ;; build IAGA-format file name

      iagaCode = thisFile.stationProp.iagaCode

      iagaFileName = string(format='(i4.4,i2.2,i2.2,i2.2)', year1, month1, day1, day2) + 'SC.' + iagaCode

      iagaFilePath = outputPathPrefix + iagaFileName

      ;; ammend html

      href = outputURLPrefix + iagaFileName

      if thisFile.exists then begin

        iagaHtml[stationIndex] = iagaHtml[stationIndex] + a(href=href, content=iagaFileName)

        iagaHtml[stationIndex] = iagaHtml[stationIndex] + br()

        ;; process data

        writeIAGAFile, iagaFilePath, thisFile

      endif else begin

        iagaHtml[stationIndex] = iagaHtml[stationIndex] + iagaFileName + ' (unavailable)'

        iagaHtml[stationIndex] = iagaHtml[stationIndex] + br()

      endelse

    endif

    if lastFileForThisStation then begin

      ;; ...when i get to the bottom i go back to the top of the slide...

      firstFileForThisStation = 1

      lastFileForThisStation = 0

    endif

endfor

caldat, beginJDN, month1, day1, year1, hour1
caldat, endJDN, month2, day2, year2, hour2

if gMakePNG or gMakePS then begin

endif

if gMakePNG then begin

    pngFileName = plotFileName + '.png'

    pngFilePath = outputPathPrefix + pngFileName

    pngFileURL = outputURLPrefix + pngFileName

    content = img(src=pngFileURL, alt=pngFileName, border='0')

    pngHtml = '<CENTER>' + a(href=pngFileURL, content=content) + '</CENTER>'

    ;; process data

    writePlotFile, 'PNG', pngFilePath, inputFilesArray, beginJDN, endJDN, availableStations

endif

if gMakePS then begin

    psFileName = plotFileName + '.ps'

    psFilePath = outputPathPrefix + psFileName

    psFileURL = outputURLPrefix + psFileName

    psHtml = a(href=psFileURL, content=psFileName)

    ;; process data

    writePlotFile, 'PS', psFilePath, inputFilesArray, beginJDN, endJDN, availableStations

endif

if gMakeBin or gMakeIAGA then begin

endif

if gMakePNG then begin

endif

if gMakePS then begin

endif

end

;*******************************************************************************
pro writeIAGAFile, iagaFilePath, inFileStruct

;; open IAGA output file
openw, iagaFileLUN, iagaFilePath, /get_lun

stationProp = inFileStruct.stationProp
iagaName = stationProp.iagaName
iagaCode = stationProp.iagaCode
latitude = string(format='(f5.1)', stationProp.latitude)
longitude = string(format='(f5.1)', stationProp.longitude)

; data column headings
Hx = " " + iagaCode + "X"
Hy = " " + iagaCode + "Y"
Hz = " " + iagaCode + "Z"

time = systime(0)

; pad with spaces
totalChars = 36
padString, iagaName, totalChars
padString, iagaCode, totalChars
padString, latitude, totalChars
padString, longitude, totalChars
totalChars = 40
padString, time, totalChars

printf, iagaFileLUN, " Format | IAGA2000 |"
printf, iagaFileLUN, " Station | " + iagaName + "|"
printf, iagaFileLUN, " IAGA Code | " + iagaCode + "|"
printf, iagaFileLUN, " Geodetic Latitude | " + latitude + "|"
printf, iagaFileLUN, " Geodetic Longitude | " + longitude + "|"
printf, iagaFileLUN, " Elevation | |"
printf, iagaFileLUN, " Observed | XYZ |"
printf, iagaFileLUN, " Digital Sampling | 0.5 seconds |"
printf, iagaFileLUN, " Data Interval | Average 5-second centered |"
printf, iagaFileLUN, " Data type | Variation |"
printf, iagaFileLUN, " Variables | 9 |"
printf, iagaFileLUN, " # This file was written by the MACCS Data Request Web |"
printf, iagaFileLUN, " # Application on " + time + "|"
printf, iagaFileLUN, " # For more information visit [WWW]http://space.augsburg.edu |"
; " Year| MT| DY| HR| MN| SC| ???X| ???Y| ???Z|"
; " 2000| 1| 1| 0| 0| 0| 15262.9| 974.9| 48404.2|"
printf, iagaFileLUN, " Year| MT| DY| HR| MN| SC|" + Hx + "|" + Hy + "|" + Hz + "|"

fillMagneticFieldArrays, bx, by, bz, hour, minute, second, inFileStruct
sizeArray = size(bx)
count = sizeArray[1]

month = indgen(count)
day = indgen(count)
year = indgen(count)

caldat, inFileStruct.jdn, m, d, y

month[*] = m
day[*] = d
year[*] = y

table = findgen(9, count)
table[0, *] = year
table[1, *] = month
table[2, *] = day
table[3, *] = hour
table[4, *] = minute
table[5, *] = second
table[6, *] = bx
table[7, *] = by
table[8, *] = bz

printf, format='(((x, i4, "|"), 5(x, i2, "|"), 3(x, f9.2, "|")))', iagaFileLUN, $

printf, iagaFileLUN, "EOF"

free_lun, iagaFileLUN

end

;*******************************************************************************
; plotType is 'PNG' or 'PS'
;
pro writePlotFile, plotType, plotFilePath, inputFilesArray, beginJDN, endJDN, availableStations

common OutputsCommons

kMissing = 32767.0
kRecordCount = 17280

if plotType eq 'PNG' then begin

endif
if plotType eq 'PS' then begin

endif

caldat, beginJDN, month1, day1, year1, hour1, minute1
caldat, endJDN, month2, day2, year2, hour2, minute2
hours = round((endJDN - beginJDN)*24.0)

case 1 of

endcase

step = round(1/tph)
for i=0, majors do begin

endfor
!x.range = [hour1, hour1 + hours]
!x.ticks=majors
!x.minor=minors
!x.style=1

sizeArray = size(availableStations)
numAvailableStations = sizeArray[1]
stationCount = 0
for index = 0, (numAvailableStations - 1) do begin

endfor

yticks = stationCount + 1
yticknames = strarr(yticks + 1)
yticknames[0] = ' '
yticknames[yticks] = ' '
orderedTicknames = ['CHB', 'CDR', 'PGG', 'RBY', 'PEB', 'GJO', 'IGL', 'CRV', 'NAN']
yticknameIndex = 1
for orderIndex = 0, 8 do begin

    for stationIndex = 0, (numAvailableStations - 1) do begin

      station = availableStations[stationIndex]

      if station eq orderedTicknames[orderIndex] then begin

        yticknames[yticknameIndex] = station

        yticknameIndex = yticknameIndex + 1

      endif

    endfor

endfor

yscale = 200 ; todo calculate yscale?
yrange = yscale*yticks
!y.ticks=yticks
!y.tickname=yticknames
!y.range=[0, yrange]
!y.minor=-1
!y.style=1

plotInitialized = 0
firstFileForThisStation = 1
lastFileForThisStation = 0
fileNumberForThisStation = 0

sizeArray = size(inputFilesArray)
numInputFiles = sizeArray[1]

firstFileIndex = 0
lastFileIndex = numInputFiles - 1

for fileIndex = firstFileIndex, lastFileIndex do begin

    thisFile = inputFilesArray[fileIndex]

    thisMaccsCode = thisFile.stationProp.maccsCode

    fileNumberForThisStation = fileNumberForThisStation + 1

    ;; see if this file exists

    if thisFile.exists then begin

      thisFileExists = 1

    endif else begin

      thisFileExists = 0

    endelse

    ;; see if this is the last file for this station

    if fileIndex lt lastFileIndex then begin

      nextFile = inputFilesArray[fileIndex + 1]

      nextMaccsCode = nextFile.stationProp.maccsCode

      if nextMaccsCode ne thisMaccsCode then begin

        lastFileForThisStation = 1

      endif

    endif else begin

      lastFileForThisStation = 1

    endelse

    if thisFileExists then begin

      ;; get data

      fillMagneticFieldArrays, bx, by, bz, hour, minute, second, thisFile

      if gBx eq 1 then begin

        b = bx

      endif

      if gBy eq 1 then begin

        b = by

      endif

      if gBz eq 1 then begin

        b = bz

      endif

      sizeArray = size(b)

      bCount = sizeArray[1]

      x = findgen(bCount)

      day = fileNumberForThisStation - 1

      x = day*24.0 + hour + minute/(60.0) + second/(60.0 * 60.0)

    endif else begin

      b = findgen(kRecordCount)

      b[*] = kMissing

      x = findgen(kRecordCount)

      day = fileNumberForThisStation - 1

      x = day*24.0 + x*24.0/kRecordCount

    endelse

    ;; figure out which tick mark this data corresponds to

    ytick = 1

    for i=1, yticks - 1 do begin

      if yticknames[i] eq thisMaccsCode then begin

        ytick = i

      endif

    endfor

    if firstFileForThisStation then begin

      ;; assuming b and x arrays have same size

      samplingRate = 1.0/5.0

      lag = beginJDN - thisFile.jdn

      firstIndex = round(lag * 24.0 * 60.0 * 60.0 * samplingRate)

      bAll = b[firstIndex:*]

      xAll = x[firstIndex:*]

      firstFileForThisStation = 0

    endif else begin

      bAll = [bAll, b]

      xAll = [xAll, x]

    endelse

    if lastFileForThisStation then begin

      bGood = where(bAll ne kMissing, bGoodCount)

      if (bGoodCount ne 0) then begin

        bAvg = total(bAll(bGood))/bGoodCount

      endif else begin

        bAvg = !values.d_nan

      endelse

      ;; filter out missing data by assigning NaN

      bBad = where(bAll eq kMissing, bBadCount)

      if (bBadCount ne 0) then begin

        bAll(bBad) = !values.d_nan

      endif

      kMaxVal = 10000.0

      kMinVal = -10000.0

      if plotInitialized then begin

        oplot, xAll, (bAll - bAvg + ytick*yscale), max_value=kMaxVal, min_value=kMinVal

      endif else begin

        if plotType eq 'PS' then begin

          ;; We need to specify 6.5" time axis width, but there's no

          ;; apparently easy way to do this. The 'xmargin' keyword of

          ;; the 'axis' procedure allows the width to be specified in

          ;; character units only. The 'position' keyword in 'plot'

          ;; procedure doesn't specify actual axis length since it

          ;; includes the label areas. We'll just have to fudge.

          ;; Using a ruler, the label area was measured to be 5/8 of

          ;; an inch. We'll just excise it from the right margin...

          fudge = 5.0/8.0

          pos = [1.0*psxdpi, 1.0*psydpi, (7.5 + fudge)*psxdpi, 10.0*psydpi]

          plot, xAll, (bAll - bAvg + ytick*yscale), position=pos, /device, max_value=kMaxVal, min_value=kMinVal

        endif else begin

          plot, xAll, (bAll - bAvg + ytick*yscale), max_value=kMaxVal, min_value=kMinVal

        endelse

        plotInitialized = 1

      endelse

    endif ;;lastFileForThisStation

    if lastFileForThisStation then begin

      ;; ...when i get to the bottom i go back to the top of the slide...

      firstFileForThisStation = 1

      lastFileForThisStation = 0

      fileNumberForThisStation = 0

    endif

endfor

TheEnd:

if gBx eq 1 then begin

    compName = 'Bx'

endif
if gBy eq 1 then begin

    compName = 'By'

endif
if gBz eq 1 then begin

    compName = 'Bz'

endif

title = 'MACCS ' + compName + ' ' $

    + string(format='(C(CYI4.4, "/", CMOI2.2, "/", CDI2.2, X, CHI2.2))', beginJDN) $

    + ' to ' $

    + string(format='(C(CYI4.4, "/", CMOI2.2, "/", CDI2.2, X, CHI2.2))', endJDN)

xyouts, 0.50, 0.98, '!N!3'+title, align=0.5, /normal
xyouts, 0.06, 0.50, '!N!3' + 'Station Code', orientation=90, align=0.5, size=1.0, /normal
xyouts, 0.50, 0.02, '!N!3' + 'UT', align=0.5, size=1.0, /normal
xyouts, 0.04, 0.00, '!N!3'+'SCALE: '+strcompress(string(yscale))+'nT/TickMark', align=0.0, size=0.8, /normal
xyouts, 0.98, 0.02, '!N!3BOSTON UNIVERSITY', align=1.0, size=0.8, /normal
xyouts, 0.98, 0.00, '!N!3AUGSBURG COLLEGE', align=1.0, size=0.8, /normal

if plotType eq 'PNG' then begin

    ;; toggle top-to-bottom transfer order

; !order = not !order ; bug fixed in IDL 5.4

    theImage=tvrd()

    ;; retore top-to-bottom transfer order

; !order = not !order ; bug fixed in IDL 5.4

    write_png, plotFilePath, theImage

    free_lun, pngFileLUN

    device, /close

endif
if plotType eq 'PS' then begin

    device, /close_file

endif

end

;*******************************************************************************
; returns bx, by, and bz arrays for the specified input file
pro fillMagneticFieldArrays, bx, by, bz, hour, minute, second, inFileStruct

kRecordSize=15
kRecordCount=17280
; to speed up for testing...
; kRecordCount=kRecordCount/12

inFileName = inFileStruct.pathPrefix + inFileStruct.name
openr, inFileLUN, inFileName, /get_lun
data = bytarr(kRecordSize, kRecordCount)
readu, inFileLUN, data
free_lun, inFileLUN

; extract time value arrays from the file data array
hour = data(0,*)
minute = data(1,*)
second = data(2,*)
; extract field value arrays from the file data array
bx = long(data( 3: 6, *), 0, 1, kRecordCount)
by = long(data( 7:10, *), 0, 1, kRecordCount)
bz = long(data(11:14, *), 0, 1, kRecordCount)
; swap bytes if big-endian (necessary for solaris on sparc since data is little
;
; non-student version
;
byteorder, bx, /lswap, /swap_if_big_endian
byteorder, by, /lswap, /swap_if_big_endian
byteorder, bz, /lswap, /swap_if_big_endian
; divide field values by 1000. to get nanotesla
bx = bx/1000.
by = by/1000.
bz = bz/1000.

; flatten arrays into one dimension
hour = reform(hour, kRecordCount, /overwrite)
minute = reform(minute, kRecordCount, /overwrite)
second = reform(second, kRecordCount, /overwrite)
bx = reform(bx, kRecordCount, /overwrite)
by = reform(by, kRecordCount, /overwrite)
bz = reform(bz, kRecordCount, /overwrite)

end

;*******************************************************************************
; pads the specified string with spaces
pro padString, str, chars

if strlen(str) gt chars then begin

    ;; truncate

    str = strmid(str, 0, chars)

    goto, TheEnd

endif else begin

    ;; pad

    while strlen(str) lt chars do begin

      str = str + ' '

    endwhile

endelse

TheEnd:
theString = str

end

This is a Wiki Spot wiki. Wiki Spot is a non-profit organization that helps communities collaborate via wikis.