Boreal triaxial projection

import math

pi=3.14159265358979
degrees=pi/180

############################ ENTER DATA ############################
mapin=open('c:\KPV\Maps\Wikipedia blank maps\BlankMap-Equirectangular.svg', 'rb')
mapout=open('c:\KPV\Maps\Wikipedia blank maps\BlankMap-BorealTriaxial.svg', 'wb')
####################################################################

l4c="1234" # the last four characters
readcoordinates=False
watchit=False
numberstatus=0  # 0: not within a coordinate pair
                # 1: read the first number
                # 2: between the two coordinates
                # 3: read the second number
num1=''
num2=''
numchars=['0','1','2','3','4','5','6','7','8','9','.','-','+','e']

mainlon=[-70*degrees,+20*degrees,+110*degrees]
midlon=[-25*degrees,+65*degrees,-160*degrees]

def roundit (sn=3.14159265358979):
    rnum=round(sn,5) # the number stands for the number of digits after the point
    position=str(sn).find('.')
    snum=str(sn)[:position+6]
    return snum

def rounddist(rda,rdb):
    rdd=rda-rdb
    if rdd>pi: rdd=rdd-2*pi
    if rdd<-pi: rdd=rdd+2*pi
    return rdd

def calculatecoordinates (x_i, y_i):
    R = pi/2 - y_i
    s0 = 0
    s1 = 0
    for i in range(3):
        ximi0=rounddist(x_i,mainlon[i])
        xims0=rounddist(x_i,mainlon[s0])
        if abs(ximi0)<=abs(xims0): s0 = i
        #print i, ximi0/degrees, xims0/degrees
    xims0=rounddist(x_i,mainlon[s0])
    for i in range(3):
        ximi1=(rounddist(x_i,midlon[i]))
        xims1=(rounddist(x_i,midlon[s1]))
        if (ximi1*xims0<=0) and (abs(ximi1)<=abs(xims1)): s1 = i
        #print i, ximi1/degrees, xims1/degrees
    sh=rounddist(x_i,mainlon[s0])/rounddist(midlon[s1],mainlon[s0])
    sign = 1
    if sh < 0:
        sign = -1
        sh = -sh
    if sh > 1: sh = 1
    sh = 1 - (1-sh) ** (1-R/pi)
    alpha = sign * sh * (rounddist(midlon[s1],mainlon[s0]))+mainlon[s0]-mainlon[1]
    x = R * math.sin(alpha)
    y = -R * math.cos(alpha)
    mapout.write(roundit(x/degrees)+","+roundit(y/degrees))

while True:
    try:
        ch=mapin.read(1)
        if ch=='': break
        if ch=='$': watchit=True
        l4c=l4c[1:]+str(ch)
        if readcoordinates and (ch=='"'):
            readcoordinates=False
            if (l4c[2] in numchars) and (numberstatus==3):
                numberstatus=0
                calculatecoordinates(eval(num1)*pi/180,(eval(num2))*pi/180)
                num1=""
                num2=""
        if l4c==' d="': readcoordinates=True
        if readcoordinates:
            if (numberstatus==0) and (ch in numchars): numberstatus=1
            if (numberstatus==1) and (not (ch in numchars)): numberstatus=2
            if (numberstatus==2) and (ch in numchars): numberstatus=3
            if (numberstatus==3) and (not (ch in numchars)):
                numberstatus=0
                calculatecoordinates((eval(num1))*pi/180,(eval(num2))*pi/180)
                num1=""
                num2=""
            if numberstatus==1: num1=num1+ch
            if numberstatus==3: num2=num2+ch
            if numberstatus==0: mapout.write(ch)
        else:
            mapout.write(ch)
    except KeyboardInterrupt: break

mapout.close()
mapin.close()