pyproj is a Python wrapper around PROJ.4. Here's a quick walkthrough.

Initialize a geodetic converter:

>>> from pyproj import Geod
>>> g = Geod(ellps='clrk66')

where ellps='clrk66' selects Clarke's 1866 reference ellipsoid. help(Geod.__new__) gives a list of possible ellipsoids.

Calculate the distance between two points, as well as the local heading, try

>>> lat1,lng1 = (40.7143528, -74.0059731)  # New York, NY
>>> lat2,lng2 = (49.261226, -123.1139268)   # Vancouver, Canada
>>> az12,az21,dist = g.inv(lng1,lat1,lng2,lat2)
>>> az12,az21,dist
(-59.10918706123901, 84.99453463527395, 3914198.2912370963)

which gives forward and back azimuths as well as the geodesic distance in meters. Not that longitude comes before latitude in the these pyproj argument lists.

Calculate the terminus of a geodesic from an initial point, azimuth, and distance with:

>>> lng3,lat3,az3 = g.fwd(lng1,lat1,az12, dist)
>>> lat3,lng3,az3
(49.26122600306212, -123.11392684861474, 84.99453467574762)

Plan your trip with:

>>> pts = g.npts(lng1,lat1,lng2,lat2,npts=5)
>>> pts.insert(0, (lng1, lat1))
>>> pts.append((lng2, lat2))
>>> import numpy
>>> npts = numpy.array(pts)
>>> npts
array([[ -74.0059731 ,   40.7143528 ],
       [ -80.93566289,   43.52686057],
       [ -88.48167748,   45.87969433],
       [ -96.61187851,   47.6930911 ],
       [-105.22271807,   48.89347605],
       [-114.13503215,   49.42510006],
       [-123.1139268 ,   49.261226  ]])

To plot the above New York to Vancouver route on a flat map, we need a Proj instance:

>>> from pyproj import Proj
>>> awips221 = Proj(proj='lcc', R=6371200, lat_1=50, lat_2=50,
...     lon_0=-107, ellps='clrk66')
>>> awips218 = Proj(proj='lcc', R=6371200, lat_1=25, lat_2=25,
...     lon_0=-95, ellps='clrk66')  #x_0=-llcrnrx,y_0=-llcrnry)

#llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower
#    left hand corner of projection region.

where proj='lcc' selects the Lambert conformal conic projection for the x/y points, and ellps='clrk66' selects the reference ellipsoid for the lat/lng coordinates. The other coordinates are LCC parameters that select the AWIPS 221 and AWIPS 226 projections respectively (lat_1 corresponds to Latin1, lat_2 corresponds to Latin2, and lon_0 corresponds to Lov; see this description of the two-standard-parallel LCC and its PROJ.4 parameters).

Convert our lat/lng pairs into grid points:

>>> awips221(lng1, lat1)
(2725283.842678774, 5823260.730665273)
>>> x221,y221 = awips221(npts[:,0], npts[:,1])
>>> # xy221 = numpy.concatenate((a1, a2, ...), axis=0)  # numpy-2.0
>>> xy221 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype)
>>> xy221[:,0] = x221
>>> xy221[:,1] = y221
>>> xy221
array([[ 2725283.84267877,  5823260.73066527],
       [ 2071820.3526011 ,  5892518.49630526],
       [ 1422529.71760395,  5967565.49899035],
       [  775650.03731228,  6046475.43928965],
       [  129946.46495299,  6127609.80532071],
       [ -515306.57275941,  6209785.69230076],
       [-1160447.80254759,  6292455.41884832]])

Finally, you can convert points from one projection to another.

>>> from pyproj import transform
>>> x218,y218 = transform(awips221, awips218, x221, y221)
>>> xy218 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype)
>>> xy218[:,0] = x218
>>> xy218[:,1] = y218
>>> xy218
array([[ 1834251.59591561,  4780900.70184736],
       [ 1197541.13209718,  5028862.9881648 ],
       [  542391.04388716,  5258740.71523961],
       [ -131577.34942316,  5464828.45934687],
       [ -822685.42269932,  5641393.59760613],
       [-1527077.85176048,  5783597.16169582],
       [-2239159.34620498,  5888495.91009021]])

Another useful coordinate system is the Universal Transverse Mercator projection which slices the world into zones.

>>> p = Proj(proj='utm', zone=10, ellps='clrk66')

Putting everything together, here's a route map based on digital lat/lng pairs stored in a text file:

>>> from numpy import array
>>> from pylab import plot, show
>>> from pyproj import Geod, Proj
>>> latlng = array([[float(x) for x in ln.split()]
...                for ln in open('coords', 'r')
...                if not ln.startswith('#')])
>>> g = Geod(ellps='WGS84')
>>> az12s,az21s,dists = g.inv(latlng[:-1,1], latlng[:-1,0],
...                           latlng[1:,1], latlng[1:,0])
>>> print('total distance: %g m' % dists.sum())
total distance: 2078.93 m
>>> mlng = latlng[:,1].mean()
>>> zone = int(round((mlng + 180) / 6.))
>>> p = Proj(proj='utm', zone=zone, ellps='WGS84')
>>> xs,ys = p(latlng[:,1], latlng[:,0])
>>> lines = plot(xs, ys, 'r.-')
>>> show()

I've written up a simple script using this approach: maproute.py. I've also written up a simple script to draw a map with labeled points: maplabel.py.

Note that you can easily get lat/lng pairs using geopy (ebuild in my Gentoo overlay):

>>> import geopy
>>> g = geopy.geocoders.Google()
>>> place1,(lat1,lng1) = g.geocode("New York, NY")
>>> place2,(lat2,lng2) = g.geocode("Vancouver, Canada")
>>> place1,(lat1,lng1)
(u'New York, NY, USA', (40.7143528, -74.0059731))
>>> place2,(lat2,lng2)
(u'Vancouver, BC, Canada', (49.261226, -123.1139268))

If you're looking for a more compact C++ package for geographic conversions, GeographicLib looks promising.