utm_coor

import math

PI = math.pi
FOURTHPI = PI / 4;
deg2rad = PI / 180;
rad2deg = 180.0 / PI;

Ellipsoid = []
#id, Ellipsoid name, Equatorial Radius, square of eccentricity
Ellipsoid.append({'id':-1, 'name':'Placeholder','equatorial_radius':0,'eccentricity_squared':0})
Ellipsoid.append({'id': 1, 'name':'Airy','equatorial_radius':6377563,'eccentricity_squared':0.00667054})
Ellipsoid.append({'id': 2, 'name':'Australian National','equatorial_radius':6378160,'eccentricity_squared':0.006694542})
Ellipsoid.append({'id': 3, 'name':'Bessel 1841','equatorial_radius':6377397,'eccentricity_squared':0.006674372})
Ellipsoid.append({'id': 4, 'name':'Bessel 1841 (Nambia)' ,'equatorial_radius':6377484,'eccentricity_squared':0.006674372})
Ellipsoid.append({'id': 5, 'name':'Clarke 1866','equatorial_radius':6378206,'eccentricity_squared':0.006768658})
Ellipsoid.append({'id': 6, 'name':'Clarke 1880','equatorial_radius':6378249,'eccentricity_squared':0.006803511})
Ellipsoid.append({'id': 7, 'name':'Everest','equatorial_radius':6377276,'eccentricity_squared':0.006637847})
Ellipsoid.append({'id': 8, 'name':'Fischer 1960 (Mercury)' ,'equatorial_radius':6378166,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 9, 'name':'Fischer 1968','equatorial_radius':6378150,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 10, 'name':'GRS 1967','equatorial_radius':6378160,'eccentricity_squared':0.006694605})
Ellipsoid.append({'id': 11, 'name':'GRS 1980','equatorial_radius':6378137,'eccentricity_squared':0.00669438})
Ellipsoid.append({'id': 12, 'name':'Helmert 1906','equatorial_radius':6378200,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 13, 'name':'Hough','equatorial_radius':6378270,'eccentricity_squared':0.00672267})
Ellipsoid.append({'id': 14, 'name':'International','equatorial_radius':6378388,'eccentricity_squared':0.00672267})
Ellipsoid.append({'id': 15, 'name':'Krassovsky','equatorial_radius':6378245,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 16, 'name':'Modified Airy','equatorial_radius':6377340,'eccentricity_squared':0.00667054})
Ellipsoid.append({'id': 17, 'name':'Modified Everest','equatorial_radius':6377304,'eccentricity_squared':0.006637847})
Ellipsoid.append({'id': 18, 'name':'Modified Fischer 1960','equatorial_radius':6378155,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 19, 'name':'South American 1969','equatorial_radius':6378160,'eccentricity_squared':0.006694542})
Ellipsoid.append({'id': 20, 'name':'WGS 60','equatorial_radius':6378165,'eccentricity_squared':0.006693422})
Ellipsoid.append({'id': 21, 'name':'WGS 66','equatorial_radius':6378145,'eccentricity_squared':0.006694542})
Ellipsoid.append({'id': 22, 'name':'WGS-72','equatorial_radius':6378135,'eccentricity_squared':0.006694318})
#zone-paraguay
Ellipsoid.append({'id': 23, 'name':'WGS-84','equatorial_radius':6378137,'eccentricity_squared':0.00669438})

def UTMtoLL(ellipsoid, UTMEasting, UTMNorthing, UTMZone):
        '''
        converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
        East Longitudes are positive, West longitudes are negative.
        North latitudes are positive, South latitudes are negative
        Lat and Long are in decimal degrees.
        Written by Chuck Gantz- chuck.gantz@globalstar.com
        '''
        k0 = 0.9996
        a = Ellipsoid[ellipsoid]['equatorial_radius']
        eccSquared = Ellipsoid[ellipsoid]['eccentricity_squared']
        eccPrimeSquared = 0.0
        e1 = (1-math.sqrt(1-eccSquared))/(1+math.sqrt(1-eccSquared))
        N1 = 0.0
        T1 = 0.0
        C1 = 0.0
        R1 = 0.0
        D = 0.0
        M = 0.0
        LongOrigin = 0.0
        mu = 0.0
        phi1 = 0.0
        phi1Rad = 0.0
        x = 0.0
        y = 0.0
        ZoneNumber = 0
        NorthernHemisphere = 1 # 1 for northern hemispher, 0 for southern
        x = UTMEasting - 500000.0 #remove 500,000 meter offset for longitude
        y = UTMNorthing

        # LA zona UTM viene dada por un número y una letra. El número indica el uso y la letra la zona
        ZoneNumber=int(UTMZone[0:(len(UTMZone)-1)])
        ZoneLetter=UTMZone[(len(UTMZone)-1):]
        if (ord(ZoneLetter.upper()) - ord('N')) >= 0:
                NorthernHemisphere = 1 # point is in northern hemisphere
        else:
                NorthernHemisphere = 0 # point is in southern hemisphere
                y -= 10000000.0 # remove 10,000,000 meter offset used for southern hemisphere
        # print ZoneNumber, ZoneLetter, NorthernHemisphere
        LongOrigin = (ZoneNumber - 1)*6 - 180 + 3#;  //+3 puts origin in middle of zone
        eccPrimeSquared = (eccSquared)/(1-eccSquared)
        M = y / k0
        mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))
        phi1Rad = mu+(3*e1/2-27*math.pow(e1,3)/32)*math.sin(2*mu)+(21*math.pow(e1,2)/16-55*math.pow(e1,4)/32)*math.sin(4*mu)+(151*math.pow(e1,3)/96)*math.sin(6*mu)
        phi1 = phi1Rad*rad2deg
        N1 = a/math.sqrt(1-eccSquared*math.sin(phi1Rad)*math.sin(phi1Rad))
        T1 = math.tan(phi1Rad)*math.tan(phi1Rad)
        C1 = eccPrimeSquared*math.cos(phi1Rad)*math.cos(phi1Rad)
        R1 = a*(1-eccSquared)/math.pow(1-eccSquared*math.sin(phi1Rad)*math.sin(phi1Rad), 1.5)
        D = x/(N1*k0);
        Lat = phi1Rad - (N1*math.tan(phi1Rad)/R1)*(math.pow(D,2)/2-(5+3*T1+10*C1-4*math.pow(C1,2)-9*eccPrimeSquared)*math.pow(D,4)/24+(61+90*T1+298*C1+45*math.pow(T1,2)-252*eccPrimeSquared-3*math.pow(C1,2))*math.pow(D,6)/720)
        Lat = Lat * rad2deg
        Long = (D-(1+2*T1+C1)*math.pow(D,3)/6+(5-2*C1+28*T1-3*math.pow(C1,2)+8*eccPrimeSquared+24*T1*T1)*math.pow(D,5)/120)/math.cos(phi1Rad)
        Long = LongOrigin + Long * rad2deg
        return Lat,Long

if __name__ == '__main__':
    coor = open('ccr.csv', 'r')
    lines  = coor.readlines()
    nc = open('ccrn.csv', 'a+')
    for a in lines:
        lon = int(a.split(';')[0])
        lat = int(a.split(';')[1])
        result = UTMtoLL(23, lon,lat,'21J')
        nc.write('%s;%s;%s;%s\n' % (lon,lat,result[0], result[1]))
    nc.close()
    coor.close()