Saturday, October 22, 2011

This is a small script to visualize Medtronic / Minimed CareLink blood glucose measurement data that has been exported to a CSV file to produce a nice PDF graph with gnuplot, like this: Plot
#!/usr/bin/env python

import sys
import csv
import datetime
import subprocess

# number of days to plot, going back from today
num_days = 7

if len(sys.argv) != 2:
    print 'syntax: %s input.csv' % sys.argv[0]
    sys.exit(1)

start_date = datetime.date.today() - datetime.timedelta(num_days)

out_files = {}
averages = {}

def filename(y, m, d):
    return '%02d-%02d-%02d.dat' % (y, m, d)

with open(sys.argv[1]) as f:
    # skip 11 header lines
    for i in range(11):
        f.next()

    for row in csv.DictReader(f):
        (d, m, y) = (int(x) for x in row['Date'].split('/'))
        date = datetime.date(y + 2000, m, d)
        if date >= start_date:
            # Taking the raw value is more precise
            #sg = row['Sensor Glucose (mmol/L)']
            if row['Raw-Type'] == 'GlucoseSensorData':
                raw = row['Raw-Values']
                assert raw.startswith('AMOUNT=')
                sg = int(raw[raw.find('=') + 1:raw.find(',')]) / 18.

                if (y, m, d) in out_files:
                    out = out_files[(y, m, d)]
                else:
                    out = out_files[(y, m, d)] = open(filename(y, m, d), 'wt')

                time = row['Time']
                (hours, minutes) = int(time[:2]), int(time[3:5])
                minutes -= minutes % 5 # to make windows overlap
                if (hours, minutes) in averages:
                    averages[(hours, minutes)][0] += sg
                    averages[(hours, minutes)][1] += 1
                else:
                    averages[(hours, minutes)] = [sg, 1]

                print >> out, row['Time'], sg

# close output files
for f in out_files.values():
    f.close()

# write averages
average_file = 'average.dat'
avg, avg_n = 0., 0
with open(average_file, 'wt') as f:
    for ((hours, minutes), (sigma, n)) in sorted(averages.items()):
        val = float(sigma) / n
        avg += val
        avg_n += 1
        print >> f, "%02d:%02d:00 %f" % (hours, minutes, val)
avg /= avg_n

# create gnuplot script
plot_file = 'plot.gnuplot'
eps_file = 'plot.eps'
with open(plot_file, 'wt') as f:
    plots = ', '.join('"%s" using 1:2 linewidth 2 title "%02d/%02d/%02d"' % (filename(y, m, d), d, m, y) for (y, m, d) in sorted(out_files.keys()))
    print >> f, '''
set terminal postscript eps color solid
set output "%s"
set style data lines
set title "Sensor glucose (mmol/L)"
set xdata time
set timefmt "%%H:%%M:%%S"
set format x "%%H:%%M"
set yrange [2:12]
set ytics 1
plot 4 notitle linewidth 3 linecolor rgb "black", 9 notitle linewidth 3 linecolor rgb "black", %f notitle linewidth 5 linecolor rgb "red", "%s" using 1:2 title "average" linewidth 5 linecolor rgb "red",''' % (eps_file, avg, average_file), plots

# Run gnuplot, and convert eps to pdf
subprocess.check_call('gnuplot ' + plot_file, shell=True)
subprocess.call('epstopdf ' + eps_file, shell=True)