Mississippi River Valley

A while back I wrote a program to calculate which cities are within 300 miles of the Mississppi River. I used a file called place2k.txt
from http://www.census.gov/geo/www/gazetteer/places2k.html . I only took 85 data points for the river, so the distances might have as much as 15 miles in error. I included all points within 315 miles of the river, just to be safe.

Here are my data points for the Mississppi River.

Here is the output.

from __future__ import with_statement

#infile is the output of usefuldata
def doit(infile, outfile, v=False):
    river = dmsfile("mississippi river.txt")
    n = file(outfile,'w')
    with file(infile) as f:
        header = True
        for line in f:
            line = line.rsplit("t")
            if header:
                header = False
                string = ""
                for x in line:
                    x = x.strip()
                    string += '%st' % x
                string += '%st' % "distance"
                if v:
                    print string
                n.write(string+"n")
                continue
            lat = float(line[9])
            lon = float(line[10])
            distance, point = dist2river(lon, lat, river)
            if distance < 315:
                d = (point[0] < lon) and "w" or "e" #east or west
                string = ""
                for x in line:
                    x = x.strip()
                    string += '%st' % x
                string += '%st' % str(distance)
                if v:
                    print string
                n.write(string+"n")
    n.close()

def dist2river(lon,lat,river):
    'returns the distance and closest point on the river to the point'
    smallest = distance(river[0][0],river[0][1],lon,lat) #check against initial point
    for point in river:
        d = distance(point[0],point[1],lon,lat)
        if d <= smallest: #see if this point is closer
            thepoint=point
            smallest=d
    return smallest, thepoint


def count(filename):
    'counts the number of lines in a file'
    i=0
    with file(filename) as f:
        for line in f:
            i=i+1
    print i

def dms(d,m,s):
    return d+m/60.+s/3600.

def dmsfile(filen):
    """
    Imports a coordinates from tsv file and returns it as a list
    """
    a=[]
    with file(filen) as f:
        for line in f:
            t=line.rsplit("t")
            lon=dms(int(t[0]),int(t[1]),int(t[2]))
            lat=dms(int(t[3]),int(t[4]),int(t[5]))
            a.append((-lon,lat))
    return a

def distance(Lon1,Lat1,Lon2,Lat2):
    """
    Calculates the distance between two coordinates
    Based on calcDistance in http://nmviewogc.cr.usgs.gov/javascript/aimsMap.js
    """
    import math as Math
    Lon1 = Lon1 * Math.pi / 180;
    Lon2 = Lon2 * Math.pi / 180;
    Lat1 = Lat1 * Math.pi / 180;
    Lat2 = Lat2 * Math.pi / 180;
    LonDist = Lon1-Lon2;
    LatDist = Lat1-Lat2;
    A = Math.pow(Math.sin(LatDist / 2),2) + Math.cos(Lat1) * Math.cos(Lat2) * Math.pow(Math.sin(LonDist /2),2);
    C = 2 * Math.asin(Math.sqrt(A));
    D = (3963.1676 - 13 * Math.sin((Lat1 + Lat2) / 2)) * C;
    return D

if __name__ == '__main__':
    doit("places2kmod.txt","places2kdist.txt")
Advertisements