convert latitude and longitude to x and y grid system using python

if you were to use a 3D system, these functions will do:

def arc_to_deg(arc):
    """convert spherical arc length [m] to great circle distance [deg]"""
    return float(arc)/6371/1000 * 180/math.pi

def deg_to_arc(deg):
    """convert great circle distance [deg] to spherical arc length [m]"""
    return float(deg)*6371*1000 * math.pi/180

def latlon_to_xyz(lat,lon):
    """Convert angluar to cartesian coordiantes

    latitude is the 90deg - zenith angle in range [-90;90]
    lonitude is the azimuthal angle in range [-180;180] 
    """
    r = 6371 # https://en.wikipedia.org/wiki/Earth_radius
    theta = math.pi/2 - math.radians(lat) 
    phi = math.radians(lon)
    x = r * math.sin(theta) * math.cos(phi) # bronstein (3.381a)
    y = r * math.sin(theta) * math.sin(phi)
    z = r * math.cos(theta)
    return [x,y,z]

def xyz_to_latlon (x,y,z):
    """Convert cartesian to angular lat/lon coordiantes"""
    r = math.sqrt(x**2 + y**2 + z**2)
    theta = math.asin(z/r) # https://stackoverflow.com/a/1185413/4933053
    phi = math.atan2(y,x)
    lat = math.degrees(theta)
    lon = math.degrees(phi)
    return [lat,lon]

UTM projections are in meters. So you could use something like the utm lib at this link:

https://pypi.python.org/pypi/utm

Googling python lat lon to UTM will point to several options.

UTM zones are 6 degrees of longitude wide and start from 0 at the prime meridian. The origin of each UTM zone is on the equator (x-axis) with the y-axis at the western most degree of longitude. This makes the grid positive to the north and east. You could calculate your distance from these results. Values are most accurate in the middle of the UTM zone.

You also should know what datum your original lat lon values are based on and use the same datum in your conversion.

The following gets you pretty close (answer in km). If you need to be better than this, you have to work harder at the math - for example by following some of the links given.

import math
dx = (lon1-lon2)*40000*math.cos((lat1+lat2)*math.pi/360)/360
dy = (lat1-lat2)*40000/360

Variable names should be pretty obvious. This gives you

dx = 66.299 km (your link gives 66.577)
dy = 2.222 km (link gives 2.225)

Once you pick coordinates (for example, lon1, lat1) as your origin, it should be easy to see how to compute all the other XY coordinates.

Note - the factor 40,000 is the circumference of the earth in km (measured across the poles). This gets you close. If you look at the source of the link you provided (you have to dig around a bit to find the javascript which is in a separate file) you find that they use a more complex equation:

function METERS_DEGLON(x)
{  
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return((111415.13 * cos(d2r))- (94.55 * cos(3.0*d2r)) + (0.12 * cos(5.0*d2r)));
   }
}

function METERS_DEGLAT(x)
{
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return(111132.09 - (566.05 * cos(2.0*d2r))+ (1.20 * cos(4.0*d2r)) - (0.002 * cos(6.0*d2r)));
   }
}

It looks to me like they are actually taking account of the fact that the earth is not exactly a sphere... but even so when you are making the assumption you can treat a bit of the earth as a plane you are going to have some errors. I'm sure with their formulas the errors are smaller...