AI Experimentation Example – Digital Town Planning

I thought it would be fun to show you the kind of things that I am using AI for.

Here I take a visual basic script that created a set of functions required to convert eastings and northing into Longitude and Latitude. Myself and a colleague derived this several years ago (with the help of a Google , Ordnance Survey and time, and a lot of discussion) I wrote about this originally on cloudydatablog.net link

It took myself and another colleague 2 weeks to get working properly additionally I had been thinking about getting it working for a year before that!

Roll on chatGPT and I only need myself and 5 minutes!

Function PHId(North1, N0, aFo, PHI0, n, bFo)
PHI1 = ((North1 - N0) / aFo) + PHI0
M = Marc(bFo, n, PHI0, PHI1)
PHI2 = ((North1 - N0 - M) / aFo) + PHI1
Do While Abs(North1 - N0 - M) > 0.00001
PHI2 = ((North1 - N0 - M) / aFo) + PHI1
M = Marc(bFo, n, PHI0, PHI2)
PHI1 = PHI2
Loop
PHId = PHI2
End Function
 
Function Marc(bFo, n, P1, P2)
Marc = bFo * (((1 + n + ((5 / 4) * (n ^ 2)) + ((5 / 4) * (n ^ 3))) * (P2 - P1)) - (((3 * n) + (3 * (n ^ 2)) + ((21 / 8) * (n ^ 3))) * (Sin(P2 - P1)) * (Cos(P2 + P1))) + ((((15 / 8) * (n ^ 2)) + ((15 / 8) * (n ^ 3))) * (Sin(2 * (P2 - P1))) * (Cos(2 * (P2 + P1)))) - (((35 / 24) * (n ^ 3)) * (Sin(3 * (P2 - P1))) * (Cos(3 * (P2 + P1)))))
End Function
 
Function lon(East1, North1)
a = 6377563.396
b = 6356256.91
F0 = 0.9996012717
E0 = 400000
N0 = -100000
PHI0 = 0.855211333
LAM0 = -0.034906585
aFo = a * F0
bFo = b * F0
e2 = (aFo ^ 2 - bFo ^ 2) / aFo ^ 2
n = (aFo - bFo) / (aFo + bFo)
InitPHI = PHId(North1, N0, aFo, PHI0, n, bFo)
nuPL = aFo / ((1 - (e2 * (Sin(InitPHI)) ^ 2)) ^ 0.5)
rhoPL = (nuPL * (1 - e2)) / (1 - (e2 * (Sin(InitPHI)) ^ 2))
eta2PL = (nuPL / rhoPL) - 1
M = Marc(bFo, n, PHI0, InitPHI)
Et = East1 - E0
X = ((Cos(InitPHI)) ^ -1) / nuPL
XI = (((Cos(InitPHI)) ^ -1) / (6 * nuPL ^ 3)) * ((nuPL / rhoPL) + (2 * ((Tan(InitPHI)) ^ 2)))
XII = (((Cos(InitPHI)) ^ -1) / (120 * nuPL ^ 5)) * (5 + (28 * ((Tan(InitPHI)) ^ 2)) + (24 * ((Tan(InitPHI)) ^ 4)))
XIIA = (((Cos(InitPHI)) ^ -1) / (5040 * nuPL ^ 7)) * (61 + (662 * ((Tan(InitPHI)) ^ 2)) + (1320 * ((Tan(InitPHI)) ^ 4)) + (720 * ((Tan(InitPHI)) ^ 6)))
lon = (LAM0 + (Et * X) - ((Et ^ 3) * XI) + ((Et ^ 5) * XII) - ((Et ^ 7) * XIIA))
End Function
 
Function lat(East1, North1)
a = 6377563.396
b = 6356256.91
F0 = 0.9996012717
E0 = 400000
N0 = -100000
PHI0 = 0.855211333
LAM0 = -0.034906585
aFo = a * F0
bFo = b * F0
e2 = (aFo ^ 2 - bFo ^ 2) / aFo ^ 2
n = (aFo - bFo) / (aFo + bFo)
InitPHI = PHId(North1, N0, aFo, PHI0, n, bFo)
nuPL = aFo / ((1 - (e2 * (Sin(InitPHI)) ^ 2)) ^ 0.5)
rhoPL = (nuPL * (1 - e2)) / (1 - (e2 * (Sin(InitPHI)) ^ 2))
eta2PL = (nuPL / rhoPL) - 1
M = Marc(bFo, n, PHI0, InitPHI)
Et = East1 - E0
VII = (Tan(InitPHI)) / (2 * nuPL * rhoPL)
VIII = ((Tan(InitPHI)) / (24 * rhoPL * nuPL ^ 3)) * (5 + (3 * ((Tan(InitPHI)) ^ 2)) + eta2PL - (9 * ((Tan(InitPHI)) ^ 2) * eta2PL))
IX = ((Tan(InitPHI)) / (720 * rhoPL * nuPL ^ 5)) * (61 + (90 * ((Tan(InitPHI)) ^ 2)) + (45 * ((Tan(InitPHI)) ^ 4)))
lat = (InitPHI - ((Et ^ 2) * VII) + ((Et ^ 4) * VIII) - ((Et ^ 6) * IX))
End Function
 
Function degrees(radians)
degrees = 180 * radians / 3.14159265358979
End Function
 
Function trunc(value)
If value > 0 Then
trunc = Int(value)
Else
trunc = Int(value + 1)
End If
End Function

Now I ask chatGPT to convert this function from Visual Basic for applications into Postgres SQL without resorting to any of the built in PostGIS extension functions. This is what it came up with. All not tested but just wow being able to get that back in basically the time it takes to copy and paste the starting function in is phenomenal. Now from a quick review of the code everything looks like a straightforward translation of mathematical operations but that alone would have taken me hours to get right. Its not that I wouldn’t have done it because I couldn’t have done it its a case of I wouldn’t have done it because I don’t have the time… and there for me is the difference AI is making for me – its significantly increasing the amount of work I can do.

SQL

-- Function to compute PHId
CREATE OR REPLACE FUNCTION compute_phid(
    north1 DOUBLE PRECISION, 
    n0 DOUBLE PRECISION, 
    afo DOUBLE PRECISION, 
    phi0 DOUBLE PRECISION, 
    n DOUBLE PRECISION, 
    bfo DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
DECLARE
    phi1 DOUBLE PRECISION := ((north1 - n0) / afo) + phi0;
    m DOUBLE PRECISION := compute_marc(bfo, n, phi0, phi1);
    phi2 DOUBLE PRECISION := ((north1 - n0 - m) / afo) + phi1;
BEGIN
    LOOP
        IF ABS(north1 - n0 - m) <= 0.00001 THEN
            RETURN phi2;
        END IF;
        phi2 := ((north1 - n0 - m) / afo) + phi1;
        m := compute_marc(bfo, n, phi0, phi2);
        phi1 := phi2;
    END LOOP;
END;
$$ LANGUAGE plpgsql;

-- Function to compute Marc
CREATE OR REPLACE FUNCTION compute_marc(
    bfo DOUBLE PRECISION, 
    n DOUBLE PRECISION, 
    p1 DOUBLE PRECISION, 
    p2 DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
BEGIN
    RETURN bfo * (
        ((1 + n + ((5.0 / 4.0) * (n ^ 2)) + ((5.0 / 4.0) * (n ^ 3))) * (p2 - p1)) -
        (((3.0 * n) + (3.0 * (n ^ 2)) + ((21.0 / 8.0) * (n ^ 3))) * SIN(p2 - p1) * COS(p2 + p1)) +
        ((((15.0 / 8.0) * (n ^ 2)) + ((15.0 / 8.0) * (n ^ 3))) * SIN(2.0 * (p2 - p1)) * COS(2.0 * (p2 + p1))) -
        (((35.0 / 24.0) * (n ^ 3)) * SIN(3.0 * (p2 - p1)) * COS(3.0 * (p2 + p1)))
    );
END;
$$ LANGUAGE plpgsql;

-- Function to compute longitude
CREATE OR REPLACE FUNCTION compute_longitude(
    east1 DOUBLE PRECISION, 
    north1 DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
DECLARE
    a DOUBLE PRECISION := 6377563.396;
    b DOUBLE PRECISION := 6356256.91;
    f0 DOUBLE PRECISION := 0.9996012717;
    e0 DOUBLE PRECISION := 400000;
    n0 DOUBLE PRECISION := -100000;
    phi0 DOUBLE PRECISION := 0.855211333;
    lam0 DOUBLE PRECISION := -0.034906585;
    afo DOUBLE PRECISION := a * f0;
    bfo DOUBLE PRECISION := b * f0;
    e2 DOUBLE PRECISION := (afo ^ 2 - bfo ^ 2) / afo ^ 2;
    n DOUBLE PRECISION := (afo - bfo) / (afo + bfo);
    init_phi DOUBLE PRECISION := compute_phid(north1, n0, afo, phi0, n, bfo);
    nu_pl DOUBLE PRECISION := afo / SQRT(1 - (e2 * SIN(init_phi) ^ 2));
    rho_pl DOUBLE PRECISION := (nu_pl * (1 - e2)) / (1 - (e2 * SIN(init_phi) ^ 2));
    eta2_pl DOUBLE PRECISION := (nu_pl / rho_pl) - 1;
    et DOUBLE PRECISION := east1 - e0;
    x DOUBLE PRECISION := (1 / COS(init_phi)) / nu_pl;
    xi DOUBLE PRECISION := ((1 / COS(init_phi)) / (6 * nu_pl ^ 3)) * ((nu_pl / rho_pl) + (2 * TAN(init_phi) ^ 2));
    xii DOUBLE PRECISION := ((1 / COS(init_phi)) / (120 * nu_pl ^ 5)) * (5 + (28 * TAN(init_phi) ^ 2) + (24 * TAN(init_phi) ^ 4));
    xiia DOUBLE PRECISION := ((1 / COS(init_phi)) / (5040 * nu_pl ^ 7)) * (61 + (662 * TAN(init_phi) ^ 2) + (1320 * TAN(init_phi) ^ 4) + (720 * TAN(init_phi) ^ 6));
BEGIN
    RETURN lam0 + (et * x) - ((et ^ 3) * xi) + ((et ^ 5) * xii) - ((et ^ 7) * xiia);
END;
$$ LANGUAGE plpgsql;

-- Function to compute latitude
CREATE OR REPLACE FUNCTION compute_latitude(
    east1 DOUBLE PRECISION, 
    north1 DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
DECLARE
    a DOUBLE PRECISION := 6377563.396;
    b DOUBLE PRECISION := 6356256.91;
    f0 DOUBLE PRECISION := 0.9996012717;
    e0 DOUBLE PRECISION := 400000;
    n0 DOUBLE PRECISION := -100000;
    phi0 DOUBLE PRECISION := 0.855211333;
    afo DOUBLE PRECISION := a * f0;
    bfo DOUBLE PRECISION := b * f0;
    e2 DOUBLE PRECISION := (afo ^ 2 - bfo ^ 2) / afo ^ 2;
    n DOUBLE PRECISION := (afo - bfo) / (afo + bfo);
    init_phi DOUBLE PRECISION := compute_phid(north1, n0, afo, phi0, n, bfo);
    nu_pl DOUBLE PRECISION := afo / SQRT(1 - (e2 * SIN(init_phi) ^ 2));
    rho_pl DOUBLE PRECISION := (nu_pl * (1 - e2)) / (1 - (e2 * SIN(init_phi) ^ 2));
    eta2_pl DOUBLE PRECISION := (nu_pl / rho_pl) - 1;
    et DOUBLE PRECISION := east1 - e0;
    vii DOUBLE PRECISION := TAN(init_phi) / (2 * nu_pl * rho_pl);
    viii DOUBLE PRECISION := (TAN(init_phi) / (24 * rho_pl * nu_pl ^ 3)) * (5 + (3 * TAN(init_phi) ^ 2) + eta2_pl - (9 * TAN(init_phi) ^ 2 * eta2_pl));
    ix DOUBLE PRECISION := (TAN(init_phi) / (720 * rho_pl * nu_pl ^ 5)) * (61 + (90 * TAN(init_phi) ^ 2) + (45 * TAN(init_phi) ^ 4));
BEGIN
    RETURN init_phi - ((et ^ 2) * vii) + ((et ^ 4) * viii) - ((et ^ 6) * ix);
END;
$$ LANGUAGE plpgsql;

-- Function to convert radians to degrees
CREATE OR REPLACE FUNCTION radians_to_degrees(
    radians DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
BEGIN
    RETURN radians * 180.0 / PI();
END;
$$ LANGUAGE plpgsql;

-- Function to truncate value
CREATE OR REPLACE FUNCTION trunc_value(
    value DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
BEGIN
    IF value > 0 THEN
        RETURN FLOOR(value);
    ELSE
        RETURN CEIL(value - 1);
    END IF;
END;
$$ LANGUAGE plpgsql;