No exact solution exists

There is no isometric map from the sphere to the plane. When you convert lat/lon coordinates from the sphere to x/y coordinates in the plane, you cannot hope that all lengths will be preserved by this operation. You have to accept some kind of deformation. Many different map projections do exist, which can achieve different compromises between preservations of lengths, angles and areas. For smallish parts of earth's surface, transverse Mercator is quite common. You might have heard about UTM. But there are many more.

The formulas you quote compute x/y/z, i.e. a point in 3D space. But even there you'd not get correct distances automatically. The shortest distance between two points on the surface of the sphere would go through that sphere, whereas distances on the earth are mostly geodesic lengths following the surface. So they will be longer.

Approximation for small areas

If the part of the surface of the earth which you want to draw is relatively small, then you can use a very simple approximation. You can simply use the horizontal axis x to denote longitude λ, the vertical axis y to denote latitude φ. The ratio between these should not be 1:1, though. Instead you should use cos(φ0) as the aspect ratio, where φ0 denotes a latitude close to the center of your map. Furthermore, to convert from angles (measured in radians) to lengths, you multiply by the radius of the earth (which in this model is assumed to be a sphere).

• x = r λ cos(φ0)
• y = r φ

This is simple equirectangular projection. In most cases, you'll be able to compute cos(φ0) only once, which makes subsequent computations of large numbers of points really cheap.

You need the bearing, which is the angle from point 1 point 2 on the surface of the globe. You can get that from turf.js, a powerful GIS calculation library.

Then you need a distance formula, to get the distance on the surface of the globe between the two points. Again, turf.js has a distance function.

Once you have the bearing and the distance, you've basically turned this into a 2D cartesian coordinates problem. Now you need to pick your units of x and y. What is a unit of x / y? 1 meter? 1 km? Lets say we that 1 unit of x or y is equal to 1 meter. If you find that point 2 is 35m away from point 1, with a bearing of 90 degrees, you know that the xy coordinates of point 2 would be (0,35), assuming `(0,0) for point 1. If point 2 is 1.4km away from point 1, with a bearing of 32 degrees, point 2's xy coords can be calculated with some simple trig: ( 1400 * cos(32) , 1400 * sin(32) ).

So a js function might look like this:

function getCoordsOfNextPoint( firstPoint, nextPoint, prevXY = {x: 0, y: 0} ) {

  const bearing = turf.bearing(firstPoint, nextPoint);
  const distance = turf.distance(firstPoint, nextPoint, {units: 'kilometers'});

  const xy = {
    x:  prevXY.x + distance * 1000 * Math.cos(bearing * Math.PI / 180),
    y:  prevXY.y + distance * 1000 * Math.sin(bearing * Math.PI / 180)
  }

  return xy

}

const myFirstPoint = [-75.343, 39.984]
const mySecondPoint = [-75.534, 39.123]
const secondPointXYCoords = getCoordsOfNextPoint(myFirstPoint, mySecondPoint)

You may want to tweak the units or parameters. Also watch out for the conventions on the units of what the bearing function returns, that can trip you up. (I know its tripped me up in the past). Turf is super powerful, definitely worth looking into.

Edit: stringing points together

To get from point 0 to point 1, you won't need the extra parameter prevXY. But every point after that will take the point returned by the function used for the 2 points before it. For example:

// Assuming x0y0 = { x: 0, y: 0}

const x1y1 = secondPointXYCoords = getCoordsOfNextPoint(myFirstPoint, mySecondPoint)
const x2y2 = secondPointXYCoords = getCoordsOfNextPoint(myFirstPoint, mySecondPoint, x1y1)
const x3y3 = vsecondPointXYCoords = getCoordsOfNextPoint(myFirstPoint, mySecondPoint, x2y2)

And you could write a loop to read your latlng points from an array and write to another array of xy points, which will basically translate your latlng points to xy points.

Hopefully that gets your started.

With a bit of maths (to follow) you can get the following:

Longitude

=ATAN(F2/E2)

where E2 and F2 are X and Y

Latitude

=ACOS(SQRT(E2^2+F2^2)/6371)

So below I have started with the latitude and longitude of a famous northern city, converted them to X and Y coordinates, then converted them back into latitude and longitude. Because of the square root, there is actually a mirror image in the southern hemisphere.

The original equations can be written as:

X=R cos(lat) cos(long) ...(1)

Y=R cos(lat) sin(long) ...(2)

where R is the radius of the Earth in miles.

Dividing (2) by (1) gives

Y/X = tan (long)

or

long = atan(Y/X) ....(3)

Substituting (3) into (1) gives

X = R cos(lat) cos(arctan(Y/X))

but

cos(arctan(Y/X)) => X/sqrt(X^2 + Y^2) 

(by considering a right-angled triangle with shorter sides X and Y)

so

X=R cos(lat) X / sqrt(X^2+Y^2)

cos(lat) = sqrt(X^2+Y^2) / X

lat = acos(sqrt(x^2+Y^2)/X) ....(4)

Equations (3) and (4) are used to determine the latitude and longitude from X and Y.

You have so many parentheses in your latitude formula that it’s hard to see what goes with what.

Let $\phi = \mathrm{LatY} \times \frac\pi{180},$ that is, if LatY is the latitude in degrees then $\phi$ is the latitude in radians. Let $h$ be mapHeight and let $w$ be mapWidth. Then your formula for $y$ becomes this in mathematical notation: $$ y = \frac h2 - \frac{w \ln\left(\tan\left(\frac\pi4 + \frac\phi2\right)\right) }{2 \pi} $$

This is similar to the formula found at http://mathworld.wolfram.com/MercatorProjection.html except for the scaling and translation factors (which you want in order to fit the output on the display).

Solving the equation for $\phi$ (still in radians), $$ \phi = 2 \arctan\left(\exp\left(\frac{2\pi}{w} \left(\frac h2 - y \right) \right)\right) - \frac\pi2. $$ Multiply by $\frac{180}{\pi}$ to get the answer in degrees.

Again it’s hard to be sure due to the profusion of parentheses, but the attempted formula seems to be equivalent to the mathematical equation

$$ \mathrm{lat} = \frac{\exp\left(-\frac{\left(Y - \frac h2 \right)}{w} \times 2 \pi \right) - \tan\left(\frac\pi4\right) \times 2}{\pi / 180}, $$

which is clearly quite different. The fact that $\tan\left(\frac\pi4\right)$ (which is just equal to $1$) occurs in there should be a red flag indicating that something was done in the wrong order.

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...

As has already been noted in comments, the real solution to your problem is to research map projections, choose an appropriate projection, and use that projection to convert latitude and longitude to points on your display. However, for small areas of the earth's surface, a simple "flat earth" approximation might be sufficient.

The major problem with your code is that the length of a degree of latitude is not the same as the length of a degree of longitude, except at the equator. If you move one degree north, you move the same distance regardless of your location: 1/360th of the length of a great circle. But if you move one degree west, you have moved a shorter distance depending on your latitude, because you are moving on a "small circle" instead of a great circle. If you are near the north pole then the radius of this small circle is much less than the radius of a great circle. The only exception is if you are at the equator, and then the "small circle" is actually a great circle. So how to compensate?

If your latitude is $\phi$, where $\phi = 0^\circ$ at the equator and $\phi = 90^\circ$ at the north pole, then the radius of the small circle through your location is $R \cos \phi$, where $R$ is the radius of the earth. The radius of a great circle is $R$. If you move one degree north then you have moved a distance of $2 \pi R / 360$, whereas if you move one degree west you have moved a distance of $(2 \pi R \cos \phi)/360$. In your code, a pixel corresponds to some measure of distance, so "degrees per pixel" is inversely proportional to distance; but the distance per degree is different in the X and Y directions, as explained above. So to compensate, you should divide the degrees per pixel in the X direction by $\cos \phi$. Note that the cosine function in most programming languages requires an angle in radians, so you may need to convert degrees to radians before computing the cosine.

Update

If utm coordinates are suitable for your purposes, then note that in lieu of pyproj, the utm package is now available and is simpler to use. Just be aware of lng, lat ordering vs. easting, northing ordering for the respective from_latlon and to_latlon functions, per the docs.

The below conversion example would become:

import utm
x, y = utm.from_latlon(input_lat, input_lon)

Original

Use pyproj for converting your lng, lat pairs to a projected coordinate system. In other words you need to convert from your geographic coordinate system (most likely EPSG code 4326) to a local projected coordinate system, e.g. a local UTM zone or regional system such as the British National Grid (EPSG code 27700).

import pyproj as proj

# setup your projections
crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic
crs_bng = proj.Proj(init='epsg:27700') # use a locally appropriate projected CRS

# then cast your geographic coordinate pair to the projected system
x, y = proj.transform(crs_wgs, crs_bng, input_lon, input_lat)

Note that pyproj.transform() also works on numpy arrays, so you can therefore transform your lon, lat arrays to x, y arrays. You can then use numpy's built-in array methods to normalise your values.

For example:

import numpy as np
x = (x - x.min()) / (x.max() - x.min())
y = (y - y.min()) / (y.max() - y.min())

However, if you are using sklearn already, then you may as well use sklearn.preprocessing.normalize.

Here's what worked for me, without so much bs.

int x =  (int) ((MAP_WIDTH/360.0) * (180 + lon));
int y =  (int) ((MAP_HEIGHT/180.0) * (90 - lat));

The lat,lon coordinates were given to me by Android devices. So they should be in the same standard used by all Google Earth/Map products.