Skip to content

15-matplotlib.org

Latest commit

 

History

History
245 lines (183 loc) · 6.26 KB

15-matplotlib.org

File metadata and controls

245 lines (183 loc) · 6.26 KB
Raw

Class 15: matplotlib 1 - graphing

Introduction

This class is about starting to plot data! Matplotlib is a super powerful plotting system. I am definitely not an expert. Check out the range of examples:

http://matplotlib.sourceforge.net/

Setup

I have pre-parsed the data for you.

mkdir -p ~/class/15
cd ~/class/15
wget http://vislab-ccom.unh.edu/~schwehr/rt/examples/2011-10-11.gga.dat.bz2

There is no need to uncompress the data!

bzcat 2011-10-11.gga.dat.bz2 | head
# x y z quality satellites hdop
-70.9395833333 43.1354166667 35.7 2 9 1.1
-70.9395766667 43.135415 36.1 2 9 1.1
-70.93957 43.1354133333 36.5 2 9 1.1
-70.9395666667 43.1354133333 37.0 2 9 1.1
-70.9395633333 43.1354133333 37.4 2 9 1.1
-70.9395633333 43.1354133333 37.8 2 9 1.1
-70.9395616667 43.1354133333 38.3 2 9 1.1
-70.9395616667 43.135415 38.7 2 9 1.1
-70.93956 43.1354133333 39.1 2 9 1.1

Making sure that you have numpy, scipy, and matplotlib installed

Loading the data

Start ipython with pylab loaded

ipython --pylab

We can load the data using numpy’s loadtxt function. It knows how to handle gzip and bzip2 compressed data!

loadtxt?
# q to quit out of the pager

data = loadtxt('2011-10-11.gga.dat.bz2')

data

# The results
array([[-70.93958333,  43.13541667,  35.7       ,   2.        ,
          9.        ,   1.1       ],
       [-70.93957667,  43.135415  ,  36.1       ,   2.        ,
          9.        ,   1.1       ],
       [-70.93957   ,  43.13541333,  36.5       ,   2.        ,
          9.        ,   1.1       ],
       ..., 
       [-70.93955167,  43.13545667,  43.        ,   2.        ,
          8.        ,   1.2       ],
       [-70.93955167,  43.13545667,  43.3       ,   2.        ,
          8.        ,   1.2       ],
       [-70.93955167,  43.13545833,  43.5       ,   2.        ,
          8.        ,   1.2       ]])

That’s not really what we want. We would like an array for each variable. There is a “unpack” option for loadtxt that will let us assign each of the 6 columns separately.

x,y,z,quality,satellites,hdop = loadtxt('2011-10-11.gga.dat.bz2', unpack=True)

average(x)
# -70.939601490675187

average(y)
#: 43.135434976022353

min(x)
# -70.939716666699994

max(x)
# -70.939486666700006

Simple plots

That’s nice, but can we plot those?

plot(x)

Ugly units, but it does give us a sense of what is going on. It looks pretty random. How about latitude?

plot(y)

Uh oh! Having both on the same plot is not what we want. Time to learn how to manage “figures”.

cla() # wipe what is in the "current" figure
figure(1)
plot(x)
figure(2) # make a new figure
plot(y)
figure(3)
plot(x,y)

Can we lable the hair ball?

You bet!

title('GPS wander for 1 day')
xlabel('Longitude')
ylabel('Latitude')

We can also put a dot at each point:

plot(x,y,'ro')  # Yuck!  too many points

Try something nicer - the lines with the average of x and y marked:

cla()
plot(x,y)
annotate('Center', xy = (average(x),average(y)))
plot(average(x),average(y), 'ro')

Using pyproj to calculate distances and directions

pyproj can reproject data, but it also has functions for calculating distances and bearings on great circles. We can use this to plot the direction and distance of the GPS wander during a day in the xy plane. Let’s give it a try with two points:

import pyproj

geod = pyproj.Geod(ellps='WGS84')

geod.inv(-70.93957666666667, 43.135415, -70.939585, 43.135403333333336)
# (-152.38532556888342, 27.614668733409243, 1.464169686122718)

The results are the direction from pt 1 to pt 2, the reverse direction, and the distance in meters.

Let’s make a script that will use pyproj to give us two lists: the direction and the distance. We will use the average of x and y as our reference point. Call is “wander.py”:

FIX: add explanation of import as

import pyproj
import numpy as np

def wander_list():

    geod = pyproj.Geod(ellps='WGS84')

    x,y,z,quality,satellites,hdop = np.loadtxt('2011-10-11.gga.dat.bz2', unpack=True)

    x_ave = np.average(x)
    y_ave = np.average(y)

    delta_m = []
    delta_dir = []
    for i in range(len(x)):
        d = geod.inv(x_ave,y_ave, x[i], y[i])
        delta_dir.append(d[0])
        delta_m.append(d[2])

    return delta_dir, delta_m

What????

ValueError: undefined inverse geodesic (may be an antipodal point)

If the values are too close together, it can’t calculate a direction. We have not talked about exceptions up to this point. We now need to learn a little about them. In python, when things go wrong, well designed code will often “raise” an exception. It is saying “hey! this is what is wrong”. The hope is that some of the code above that problem will know what to do. If no part of the program “catches” the exception and handles it, the program has to stop. Here is a little example to use in ipython:

try:
    raise Exception
    # Never get to here
except:
    # Do what is necessary to get things back on track
    print 'Opps!'

Modify the code in your for loop to look like this:

try:
    d = geod.inv(x_ave,y_ave, x[i],y[i])
    delta_dir.append(d[0])
    delta_m.append(d[2])
except:
    delta_dir.append(0) # just pick north
    delta_m.append(0) # No distance

Now it should work!