Convert Eastings and Northings to Lat Long in Power Query

[macmy034]

Hey all;

Just thought I would drop this here in a hope of helping out people in the future. I have a number of projects that require me to use Easting and Northing coords; the issue is everyone wants it mapped onto a map - which is hard to do when there is no actual reference to the map!

So on my little journey of discovery; I have found several scripts along the way that helped me in a few ways, but not quite what I wanted. I finally found a function that does this in SQL but pushing things up to a SQL server and back down was slow and painful, and on occasion not return all the data!

So converting the SQL Function to Power Query seems the best option, its very fast - the only issue is you need to adjust the meridian line to the area you want to convert - If your doing this on a world wide scale, you would have to add another var and adjust it dynamically; alternatively you could have a lookup table with a bunch of if statements prior to the actual script running.

This original SQL Script was sourced from; http://vietmaximo.blogspot.com/2020/06/convert-australian-gda94.html and he should take 100% of the credit - I've simply switched it from SQL to a Power Query Function.

Hope it helps someone!

let LatLongConversion = (E, N, Conv) =>

    let  
        a = 6378137,
        f = 1/298.257222101,
        phizero = 0,
        lambdazero = 147,               //Central Meridian for Zone 55
        Nzero = 10000000,               //False Northing
        Ezero = 500000,                 //False Easting
        kzero = 0.9996,                 //Scale Factor

        b = a*(1-f),
        esq = 2*f-Number.Power(f,2),
        A0 = 1-esq/4-3*Number.Power(esq,2)/64-5*Number.Power(esq,3/256),
        A2 = 0.375*(esq+Number.Power(esq,2)/4+15*Number.Power(esq,3)/128),
        A4 = 15*(Number.Power(esq,2)+3*Number.Power(esq,3)/4)/256,
        A6 = 35*Number.Power(esq,3)/3072,

        Nprime     = N-Nzero,
        mprime     = Nprime/kzero,
        smn        = (a-b)/(a+b),
        G          = a*(1-smn)*(1-Number.Power(smn,2))*(1+9*Number.Power(smn,2)/4+225*Number.Power(smn,4)/64)*Number.PI/180.0,
        sigma      = mprime*Number.PI/(180*G),
        phiprime   = sigma+(3*smn/2-27*Number.Power(smn,3)/32)*Number.Sin(2*sigma)+(21*Number.Power(smn,2)/16-55*Number.Power(smn,4)/32)*Number.Sin(4*sigma)+(151*Number.Power(smn,3)/96)*Number.Sin(6*sigma)+(1097*Number.Power(smn,4)/512)*Number.Sin(8*sigma),
        rhoprime   = a*(1-esq)/Number.Power((1-esq*Number.Power((Number.Sin(phiprime)),2)),1.5),
        upsilonprime   =a/Number.Sqrt(1-esq*Number.Power((Number.Sin(phiprime)),2)),

        psiprime   = upsilonprime/rhoprime,
        tprime     = Number.Tan(phiprime),
        Eprime     = E-Ezero,
        chi        = Eprime/(kzero*upsilonprime),
        term_1     = tprime*Eprime*chi/(kzero*rhoprime*2),
        term_2     = term_1*Number.Power(chi,2)/12*(-4*Number.Power(psiprime,2)+9*psiprime*(1-Number.Power(tprime,2))+12*Number.Power(tprime,2)),
        term_3     = tprime*Eprime*Number.Power(chi,5)/(kzero*rhoprime*720)*(8*Number.Power(psiprime,4)*(11-24*Number.Power(tprime,2))-12*Number.Power(psiprime,3)*(21-71*Number.Power(tprime,2))+15*Number.Power(psiprime,2)*(15-98*Number.Power(tprime,2)+15*Number.Power(tprime,4))+180*psiprime*(5*Number.Power(tprime,2)-3*Number.Power(tprime,4))+360*Number.Power(tprime,4)),
        term_4     = tprime*Eprime*Number.Power(chi,7)/(kzero*rhoprime*40320)*(1385+3633*Number.Power(tprime,2)+4095*Number.Power(tprime,4)+1575*Number.Power(tprime,6)),
        term1      = chi*(1/Number.Cos(phiprime)),
        term2      = Number.Power(chi,3)*(1/Number.Cos(phiprime))/6*(psiprime+2*Number.Power(tprime,2)),
        term3      = Number.Power(chi,5)*(1/Number.Cos(phiprime))/120*(-4*Number.Power(psiprime,3)*(1-6*Number.Power(tprime,2))+Number.Power(psiprime,2)*(9-68*Number.Power(tprime,2))+72*psiprime*Number.Power(tprime,2)+24*Number.Power(tprime,4)),
        term4      = Number.Power(chi,7)*(1/Number.Cos(phiprime))/5040*(61+662*Number.Power(tprime,22)+1320*Number.Power(tprime,4)+720*Number.Power(tprime,6)),

        latitude   = (phiprime-term_1+term_2-term_3+term_4)*180/Number.PI,
        longitude  = lambdazero+180/Number.PI*(term1-term2+term3-term4),

        Conversion_Output = if Conv = "lat" then latitude else if Conv = "long" then longitude else null

    in Conversion_Output

in LatLongConversion

[americanomuerto]

Thanks for this @macmy034. 

Outside of editing the central meridian, are there other fields to consider updating?

For instance, I believe my state coordinates are in feet. When I did  this conversion changing only the central meridian, all my data points were in Colombia!

[jtg80]

Thanks for providing this. It helped me start down my own little journey to convert state plane coordinates to lat/lon. Here's a couple of functions that work in Power Query - the first for lat, the second for lon.

= let

    get_latitude = (x as number, y as number) as number =>

    let

        pi = 3.14159265358979,

        fx = /* Enter false easting value */,

        fy = /* Enter false northing value */,

        loD = /* Enter latitude of origin value */,

        loR = loD * pi / 180,

        cmD = /* Enter central meridian value */,

        cmR = cmD * pi / 180,

        sp1D = /* Enter standard parallel 1 value */,

        sp1R = sp1D * pi / 180,

        sp2D = /* Enter standard parallel 2 value */,

        sp2R = sp2D * pi / 180,

        smjrM = /* Enter semimajor axis value in meters */,

        smjrF = smjrM * 3.2808333,

        smnrM = /* Enter semiminor axis value in meters */,

        smnrF = smnrM * 3.2808333,

        iflat = smjrM / (smjrM - smnrM),

        flat = 1 / iflat,

        e = Number.Sqrt(2 * flat - Number.Power(flat, 2)),

        m1 = Number.Cos(sp1R) / Number.Power((1 - Number.Power(e, 2) * Number.Power(Number.Sin(sp1R), 2)), 0.5),

        m2 = Number.Cos(sp2R) / Number.Power((1 - Number.Power(e, 2) * Number.Power(Number.Sin(sp2R), 2)), 0.5),

        t1 = Number.Tan((pi / 4) - (sp1R / 2)) / Number.Power((1 - (e * Number.Sin(sp1R))) / (1 + (e * Number.Sin(sp1R))), e / 2),

        t2 = Number.Tan((pi / 4) - (sp2R / 2)) / Number.Power((1 - (e * Number.Sin(sp2R))) / (1 + (e * Number.Sin(sp2R))), e / 2),

        TF = Number.Tan((pi / 4) - (loR / 2)) / Number.Power((1 - (e * Number.Sin(loR))) / (1 + (e * Number.Sin(loR))), e / 2),

        n = (Number.Ln(m1) - Number.Ln(m2)) / (Number.Ln(t1) - Number.Ln(t2)),

        F = m1 / (n * Number.Power(t1, n)),

        RF = smjrF * F * Number.Power(TF, n),

        RZ = Number.Power(Number.Power((x - fx), 2) + Number.Power((RF - (y - fy)), 2), 0.5),

        tz = Number.Power((RZ / (smjrF * F)), 1 / n),

        zz = Number.Atan((x - fx) / (RF - (y - fy))),

        latTR = (pi / 2) - (2 * Number.Atan(tz)),

        lat1R = (pi / 2) - (2 * Number.Atan((tz * Number.Power((1 - (e * Number.Sin(latTR))) / (1 + (e * Number.Sin(latTR))), e / 2)))),

        lat2R = (pi / 2) - (2 * Number.Atan((tz * Number.Power((1 - (e * Number.Sin(lat1R))) / (1 + (e * Number.Sin(lat1R))), e / 2)))),

        lat3R = (pi / 2) - (2 * Number.Atan((tz * Number.Power((1 - (e * Number.Sin(lat2R))) / (1 + (e * Number.Sin(lat2R))), e / 2)))),

        lat4R = (pi / 2) - (2 * Number.Atan((tz * Number.Power((1 - (e * Number.Sin(lat3R))) / (1 + (e * Number.Sin(lat3R))), e / 2)))),

        latFR = (pi / 2) - (2 * Number.Atan((tz * Number.Power((1 - (e * Number.Sin(lat4R))) / (1 + (e * Number.Sin(lat4R))), e / 2)))),

        latTD = (latTR * 180) / pi,

        lat1D = (lat1R * 180) / pi,

        lat2D = (lat2R * 180) / pi,

        lat3D = (lat3R * 180) / pi,

        lat4D = (lat4R * 180) / pi,

        latFD = (latFR * 180) / pi

    in

        latFD

in

    get_latitude

= let

    get_longitude = (x as number, y as number) as number =>

    let

        pi = 3.14159265358979,

        fx = /* Enter false easting value */,

        fy = /* Enter false northing value */,

        loD = /* Enter latitude of origin value */,

        loR = loD * pi / 180,

        cmD = /* Enter central meridian value */,

        cmR = cmD * pi / 180,

        sp1D = /* Enter standard parallel 1 value */,

        sp1R = sp1D * pi / 180,

        sp2D = /* Enter standard parallel 2 value */,

        sp2R = sp2D * pi / 180,

        smjrM = /* Enter semimajor axis value in meters */,

        smjrF = smjrM * 3.2808333,

        smnrM = /* Enter semiminor axis value in meters */,

        smnrF = smnrM * 3.2808333,

        iflat = smjrM / (smjrM - smnrM),

        flat = 1 / iflat,

        e = Sqrt(2 * flat - Number.Power(flat, 2)),

        m1 = Number.Cos(sp1R) / Number.Power((1 - Number.Power(e, 2) * Number.Power(Number.Sin(sp1R), 2)), 0.5),

        m2 = Number.Cos(sp2R) / Number.Power((1 - Number.Power(e, 2) * Number.Power(Number.Sin(sp2R), 2)), 0.5),

        t1 = Number.Tan((pi / 4) - (sp1R / 2)) / Number.Power((1 - (e * Number.Sin(sp1R))) / (1 + (e * Number.Sin(sp1R))), e / 2),

        t2 = Number.Tan((pi / 4) - (sp2R / 2)) / Number.Power((1 - (e * Number.Sin(sp2R))) / (1 + (e * Number.Sin(sp2R))), e / 2),

        TF = Number.Tan((pi / 4) - (loR / 2)) / Number.Power((1 - (e * Number.Sin(loR))) / (1 + (e * Number.Sin(loR))), e / 2),

        n = (Number.Ln(m1) - Number.Ln(m2)) / (Number.Ln(t1) - Number.Ln(t2)),

        F = m1 / (n * Number.Power(t1, n)),

        RF = smjrF * F * Number.Power(TF, n),

        RZ = Number.Power(Number.Power((x - fx), 2) + Number.Power((RF - (y - fy)), 2), 0.5),

        tz = Number.Power((RZ / (smjrF * F)), 1 / n),

        zz = Number.Atan((x - fx) / (RF - (y - fy))),

        lonR = ((zz / n) + cmR),

        lonD = (lonR * 180) / pi

    in

        lonD

in

    get_longitude

[americanomuerto]

@jtg80  Why is the conversion from meter to feet not 0.304800609601219?

[jtg80]

smnrF is something that I put in there but never used. Honestly, I probably should have removed it, But to answer your questions, if you are going from meters to feet the conversion is 3.2808333. If you are going the other way (feet to meters) then it is 0.304800609601219

[LukeAB]

hi macmy034,

sorry for being daft, but by 'converting to Power Query' do you mean you have this in Excel and somehow then link the PQ to PowerApps? Not sure I follow how that is done , can you flesh that out? 

tx 

LAB

[macmy034]

Uh no...

It was a SQL Function and I converted it to a Power Query function - It can then be used within Excel or PowerBi to convert eastings and northings to lat / long.

Nothing to do with PowerApps? 

What are you attempting to do?

[LukeAB]

ok tx. I got the link to here via PowerApps .

on PowerApps it only will capture the LatLong of the device locn. I want the UTM coords. I can do it via a call back to SQL server USP and do that way if I have to. Another way would be bodging it up in PowerApps but a bit of work to get the variable declared and then some of the functions (like Power and SqrRoot) to see if available 

[macmy034]

Hi Luke,

unfortunately I'm unfamiliar with powerapps and how to convert data within it - you could store the data in a location and run this function against it, either in powerautomate with power query or in sql at a higher level.

I have the original sql query if you need it (it should be in the above link I initially added to the post)

[LukeAB]

all good. I have in sql as well, just trying to avoid going up and back from it. But is pretty fast, so may not really be an issue, as storing the data on SQL

tx again

[Alice_Drummond]

Thanks so much!!! This is so so helpful

[whisp1986]

If you are wondering how to convert this excellent bit of PowerQuery code from AUS to BNG you have to provide the parameters from this website into the variables at the start: OSGB36 / British National Grid - United Kingdom Ordnance Survey - EPSG:27700, The LAT value will still be an issue and to get around this you need to add 49 to the output LAT value at the bottom of the code.

NB - The original SQL code states that this is approx. 100m out

I hope this helps.

let LatLongConversion = (E, N, Conv) =>

    let  
        a = 6377563.396,
        f = 1/299.3249646,
        phizero = 0,
        lambdazero = -2,               //Central Meridian for Zone 55
        Nzero = -100000,               //False Northing
        Ezero = 400000,                 //False Easting
        kzero = 0.9996,                 //Scale Factor

        b = a*(1-f),
        esq = 2*f-Number.Power(f,2),
        A0 = 1-esq/4-3*Number.Power(esq,2)/64-5*Number.Power(esq,3/256),
        A2 = 0.375*(esq+Number.Power(esq,2)/4+15*Number.Power(esq,3)/128),
        A4 = 15*(Number.Power(esq,2)+3*Number.Power(esq,3)/4)/256,
        A6 = 35*Number.Power(esq,3)/3072,

        Nprime     = N-Nzero,
        mprime     = Nprime/kzero,
        smn        = (a-b)/(a+b),
        G          = a*(1-smn)*(1-Number.Power(smn,2))*(1+9*Number.Power(smn,2)/4+225*Number.Power(smn,4)/64)*Number.PI/180.0,
        sigma      = mprime*Number.PI/(180*G),
        phiprime   = sigma+(3*smn/2-27*Number.Power(smn,3)/32)*Number.Sin(2*sigma)+(21*Number.Power(smn,2)/16-55*Number.Power(smn,4)/32)*Number.Sin(4*sigma)+(151*Number.Power(smn,3)/96)*Number.Sin(6*sigma)+(1097*Number.Power(smn,4)/512)*Number.Sin(8*sigma),
        rhoprime   = a*(1-esq)/Number.Power((1-esq*Number.Power((Number.Sin(phiprime)),2)),1.5),
        upsilonprime   =a/Number.Sqrt(1-esq*Number.Power((Number.Sin(phiprime)),2)),

        psiprime   = upsilonprime/rhoprime,
        tprime     = Number.Tan(phiprime),
        Eprime     = E-Ezero,
        chi        = Eprime/(kzero*upsilonprime),
        term_1     = tprime*Eprime*chi/(kzero*rhoprime*2),
        term_2     = term_1*Number.Power(chi,2)/12*(-4*Number.Power(psiprime,2)+9*psiprime*(1-Number.Power(tprime,2))+12*Number.Power(tprime,2)),
        term_3     = tprime*Eprime*Number.Power(chi,5)/(kzero*rhoprime*720)*(8*Number.Power(psiprime,4)*(11-24*Number.Power(tprime,2))-12*Number.Power(psiprime,3)*(21-71*Number.Power(tprime,2))+15*Number.Power(psiprime,2)*(15-98*Number.Power(tprime,2)+15*Number.Power(tprime,4))+180*psiprime*(5*Number.Power(tprime,2)-3*Number.Power(tprime,4))+360*Number.Power(tprime,4)),
        term_4     = tprime*Eprime*Number.Power(chi,7)/(kzero*rhoprime*40320)*(1385+3633*Number.Power(tprime,2)+4095*Number.Power(tprime,4)+1575*Number.Power(tprime,6)),
        term1      = chi*(1/Number.Cos(phiprime)),
        term2      = Number.Power(chi,3)*(1/Number.Cos(phiprime))/6*(psiprime+2*Number.Power(tprime,2)),
        term3      = Number.Power(chi,5)*(1/Number.Cos(phiprime))/120*(-4*Number.Power(psiprime,3)*(1-6*Number.Power(tprime,2))+Number.Power(psiprime,2)*(9-68*Number.Power(tprime,2))+72*psiprime*Number.Power(tprime,2)+24*Number.Power(tprime,4)),
        term4      = Number.Power(chi,7)*(1/Number.Cos(phiprime))/5040*(61+662*Number.Power(tprime,22)+1320*Number.Power(tprime,4)+720*Number.Power(tprime,6)),

        latitude   = 49+(phiprime-term_1+term_2-term_3+term_4)*180/Number.PI,
        longitude  = lambdazero+180/Number.PI*(term1-term2+term3-term4),

        Conversion_Output = if Conv = "lat" then latitude else if Conv = "long" then longitude else null

    in Conversion_Output

in LatLongConversion

[spandy34]

Hi @macmy034 

Thank you for this.  THis is what I am wanting but due to my knowledge of PBI being limited I am unsure of the stages I would do.  Do you have any instructions or notes as to how I would do this.  I have two columns in a table which hold the eastings and northings and I want two columns for long / lat .  I know how to create a power query but I just need help linking referencing it to my table called Service Requests.  

I live in the UK and will only be looking at the data in the UK.  Alternativley could I use the SQL code you used and create two new columns usin SQL.  Your thoughts would be appreciated.

Thanks for your time.

[macmy034]

Hey Sandy;

The SQL Code is available in the Link I initially posted (I'm on a work Pc right now, so can't access it)

alternatively if you want to do it in power query, you can follow the video I made here:

https://www.youtube.com/watch?v=N8AiLx1mLG0

this demonstrates how to calculate Latitude, do the same for Longitude but type in "long" in the "Conv" field and it will return longitude.

Small note - if your not familiar with PowerQuery functions, this is going to pass every value in the column to the function, one by one - In the Demo I'm just doing this in a table with 1 row, but if you have 200,000 - it will send all 200,000 for conversion.

Your central meridian line for london is -3. If its using the British Geological Survey measurements youll have to use the python script below.

Reach out if you get stuck 🙂

[AndrewSV]

This is incredible, thank you so much!

As I understand it, this take 3 input parameters for E, N, and Conv. Is there any way to apply the query to whole columns? For example, if I have a table with an Easting Column and a Northing Column.

[macmy034]

Yeah totally; you just need to copy and paste the query into its own "blank query" so it becomes a function, then use "Invoke custom function" under the "add column" section of Power Query.

Make sure you check the correct zone is being converted (I highlight it in the below video), you can find yours here;

https://mangomap.com/robertyoung/maps/69585/what-utm-zone-am-i-in-#

click the area your converting, and take the "Central meridian" number. In my example im in zone 55 - so its 147 (Queensland, Australia)

also youll have to run the function twice, once to calc lat, once to calc long.

vid:

https://www.youtube.com/watch?v=N8AiLx1mLG0

Hope this helps!

[AndrewSV]

Hey again,

I've come across another set of data I'm using which only lists Easting & Northing in the British National Grid system (epsg:27700). Is there a function that can convert this to WSG84 Longitude/Latitude?

Many thanks

[macmy034]

Have you tried the central median of -3?

Can you provide some sample data points so I can have a bit of a play around.

[AndrewSV]

Thanks for the speedy response! Below is a snapshot of the data with lat/long already populated in order for you to sense-check your results against the actuals:

Identifier	GridType	Easting	Northing	Longitude	Latitude
3700001	UKOS	493160	181046	-0.658677	51.520658
3700002	UKOS	493309	181271	-0.656471	51.522656
3700004	UKOS	493400	181264	-0.655161	51.522578
3700006	UKOS	493749	180917	-0.650225	51.519401
3700008	UKOS	493767	180688	-0.650026	51.51734
3700009	UKOS	493782	180655	-0.649819	51.517041
3700010	UKOS	493816	180304	-0.649422	51.51388
3700011	UKOS	493831	180248	-0.649221	51.513374
3700012	UKOS	493950	180225	-0.647513	51.513148

[macmy034]

Sorry for not getting back to you Andrew; Christmas took over and I forgot about this!

I can't seem to get it working within PowerBi with our original function, unfortunately, so I wrote a quick Python script to do it via the BGS Web apps:

import requests
import json
import pandas as pd

easting_values = [493160, 493309, 493400, 493749, 493767, 493782, 493816, 493831, 493950]
northing_values = [181046, 181271, 181264, 180917, 180688, 180655, 180304, 180248, 180225]

data = []

for easting,northing in zip(easting_values,northing_values):
    response = requests.get(f'https://webapps.bgs.ac.uk/data/webservices/CoordConvert_LL_BNG.cfc?method=BNGtoLatLng&easting={easting}&northing={northing}', verify=False)
    if response.status_code == 200:
        data_point = json.loads(response.text)
        data.append(data_point)
    else:
        print(f'Request failed with status code {response.status_code} for easting: {easting} and northing: {northing}')

df = pd.DataFrame(data)

This code returns a pandas df. You can use Power Query to call this python script and it will go off and get the values for you. Be warned, 2 things are occuring when using this tool - Your passing values to the internet and I am ignoring the verify statements (SSL junk). Your also going to find its much slower than a native query, so depending on how many values your trying to change, it could take a while.

If you want to run this directly from Power Query you can use the following M Code:

let
    Source = Python.Execute("import requests#(lf)import json#(lf)import pandas as pd#(lf)#(lf)easting_values = [493160, 493309, 493400, 493749, 493767, 493782, 493816, 493831, 493950]#(lf)northing_values = [181046, 181271, 181264, 180917, 180688, 180655, 180304, 180248, 180225]#(lf)#(lf)data = []#(lf)#(lf)for easting,northing in zip(easting_values,northing_values):#(lf)    response = requests.get(f'https://webapps.bgs.ac.uk/data/webservices/CoordConvert_LL_BNG.cfc?method=BNGtoLatLng&easting={easting}&northing={northing}', verify=False)#(lf)    if response.status_code == 200:#(lf)        data_point = json.loads(response.text)#(lf)        data.append(data_point)#(lf)    else:#(lf)        print(f'Request failed with status code {response.status_code} for easting: {easting} and northing: {northing}')#(lf)#(lf)df = pd.DataFrame(data)"),
    df1 = Source{[Name="df"]}[Value],
    #"Changed Type" = Table.TransformColumnTypes(df1,{{"DEGMINSECLNG", type text}, {"EASTING", Int64.Type}, {"LONGITUDE", type number}, {"NORTHING", Int64.Type}, {"DEGMINSECLAT", type text}, {"LATITUDE", type number}})
in
    #"Changed Type"

 This is a quick and dirty way of doing it, im sure you could clean it up a bit 🙂

plotted it looks like this in London:

If your going to run this in a production enviroment (PowerBi Workspace) youll need to have Python installed on the gateway machine as per:

https://community.powerbi.com/t5/Service/Python-personal-gateway-credentials-A-problem-occurred-whil...

Hope this helps!

[v-jayw-msft · Accepted solution]

Hi @macmy034 ,

Awesome! Thanks for sharing!

Best Regards,

Jay