在 Python 中将英国国家网格坐标转换为 WGS64 纬度和经度

分类：编程问答

如何解决在 Python 中将英国国家网格坐标转换为 WGS64 纬度和经度

我无法将东/北坐标表从英国国家网格系统转换为 WGS84 系统中的纬度/经度坐标。我已经用来自军械测量局的 Helmert 转换电子表格中的常量编写了以下函数集：https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/coordinates

import math
import pandas as pd

North = 666099.000 # Northing to be transformed
East = 253697.000 # Easting to be transformed

a = 6377563.396 # Semi-major axis for OGSB36
#a = 6378137.0000 # Semi-major axis for WGS84

b = 6356256.909 # Semi-minor axis for OGSB36
#b = 6356752.3142 # Semi-minor axis for WGS84

f0 = 0.9996012717 # Central Meridan Scale

e0 = 400000 # True origin Easting 
n0 = -100000 # True origin Northing

PHI0 = 0.855211333 # True origin latitude (Radians) i.e. N 49 0' 0''
DecimalPHI0 = 49.00000000 # True origin latitude (Degrees)

LAM0 = -0.034906585 # True origin longitude (Radians) i.e. W 2 0' 0''
DecimalLAM0 = -2.00000000 # True origin longitude (Degrees)

def InitialLat(North,n0,af0,PHI0,n,bf0):
    """
    Compute initial value for Latitude (PHI) IN RADIANS.
    Input:
     - northing of point (North) and northing of false origin (n0) in meters;
     - semi major axis multiplied by central meridian scale factor (af0) in meters;
     - latitude of false origin (PHI0) IN RADIANS;
     - n (computed from a,b and f0) and
     - ellipsoid semi major axis multiplied by central meridian scale factor (bf0) in meters.
    """
    #First PHI value (PHI1)
    PHI1 = ((North - n0) / af0) + PHI0

    def Marc(bf0,PHI1):
        """
        Compute meridional arc.
        Input:
         - ellipsoid semi major axis multiplied by central meridian scale factor (bf0) in meters;
         - n (computed from a,b and f0);
         - lat of false origin (PHI0) and initial or final latitude of point (PHI) IN RADIANS.
        """
        Marc = bf0 * (((1 + n + ((5 / 4) * (n ** 2)) + ((5 / 4) * (n ** 3))) * (PHI1 - PHI0))
        - (((3 * n) + (3 * (n ** 2)) + ((21 / 8) * (n ** 3))) * (math.sin(PHI1 - PHI0)) * (math.cos(PHI1 + PHI0)))
        + ((((15 / 8) * (n ** 2)) + ((15 / 8) * (n ** 3))) * (math.sin(2 * (PHI1 - PHI0))) * (math.cos(2 * (PHI1 + PHI0))))
        - (((35 / 24) * (n ** 3)) * (math.sin(3 * (PHI1 - PHI0))) * (math.cos(3 * (PHI1 + PHI0)))))
        return Marc
    
    # Calculate M
    M = Marc(bf0,PHI1)
    
    #Calculate new PHI value (PHI2)
    PHI2 = ((North - n0 - M) / af0) + PHI1
    
    #Iterate to get final value for InitialLat
    while abs(North - n0 - M) > 0.00001:
        PHI2 = ((North - n0 - M) / af0) + PHI1
        M = Marc(bf0,PHI2)
        PHI1 = PHI2
    
    InitialLat = PHI2
    return InitialLat

def E_N_to_Lat(East,North,a,b,e0,f0,LAM0):
    """
    Un-project Transverse Mercator eastings and northings back to latitude.
    Input:
     - eastings (East) and northings (North) in meters; _
     - ellipsoid axis dimensions (a & b) in meters; _
     - eastings (e0) and northings (n0) of false origin in meters; _
     - central meridian scale factor (f0) and _
     - latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
    """

    #Convert angle measures to radians
    Pi = math.pi
    RadPHI0 = PHI0 * (Pi / 180)
    RadLAM0 = LAM0 * (Pi / 180)
    
    # Compute af0,bf0,e squared (e2),n and Et
    af0 = a * f0
    bf0 = b * f0
    e2 = ((af0 ** 2) - (bf0 ** 2)) / (af0 ** 2)
    n = (af0 - bf0) / (af0 + bf0)
    Et = East - e0

    # Compute initial value for latitude (PHI) in radians
    PHId = InitialLat(North,RadPHI0,bf0)
    
    # Compute nu,rho and eta2 using value for PHId
    nu = af0 / (math.sqrt(1 - (e2 * ((math.sin(PHId)) ** 2))))
    rho = (nu * (1 - e2)) / (1 - (e2 * (math.sin(PHId)) ** 2))
    eta2 = (nu / rho) - 1
    
    # Compute Latitude
    VII = (math.tan(PHId)) / (2 * rho * nu)
    VIII = ((math.tan(PHId)) / (24 * rho * (nu ** 3))) * (5 + (3 * ((math.tan(PHId)) ** 2)) + eta2 - (9 * eta2 * ((math.tan(PHId)) ** 2)))
    IX = ((math.tan(PHId)) / (720 * rho * (nu ** 5))) * (61 + (90 * ((math.tan(PHId)) ** 2)) + (45 * ((math.tan(PHId)) ** 4)))
    
    E_N_Lat = (180 / Pi) * (PHId - ((Et ** 2) * VII) + ((Et ** 4) * VIII) - ((Et ** 6) * IX))
    return(E_N_Lat)

def E_N_to_Long(East,LAM0):
    """
    Un-project Transverse Mercator eastings and northings back to longitude.
    Input:
     - eastings (East) and northings (North) in meters;
     - ellipsoid axis dimensions (a & b) in meters;
     - eastings (e0) and northings (n0) of false origin in meters;
     - central meridian scale factor (f0) and
     - latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
    """
    # Convert angle measures to radians
    Pi = 3.14159265358979
    RadPHI0 = PHI0 * (Pi / 180)
    RadLAM0 = LAM0 * (Pi / 180)
    
    # Compute af0,n and Et
    af0 = a * f0
    bf0 = b * f0
    e2 = ((af0 ** 2) - (bf0 ** 2)) / (af0 ** 2)
    n = (af0 - bf0) / (af0 + bf0)
    Et = East - e0
    
    # Compute initial value for latitude (PHI) in radians
    PHId = InitialLat(North,rho and eta2 using value for PHId
    nu = af0 / (math.sqrt(1 - (e2 * ((math.sin(PHId)) ** 2))))
    rho = (nu * (1 - e2)) / (1 - (e2 * (math.sin(PHId)) ** 2))
    eta2 = (nu / rho) - 1
    
    # Compute Longitude
    X = ((math.cos(PHId)) ** -1) / nu
    XI = (((math.cos(PHId)) ** -1) / (6 * (nu ** 3))) * ((nu / rho) + (2 * ((math.tan(PHId)) ** 2)))
    XII = (((math.cos(PHId)) ** -1) / (120 * (nu ** 5))) * (5 + (28 * ((math.tan(PHId)) ** 2)) + (24 * ((math.tan(PHId)) ** 4)))
    XIIA = (((math.cos(PHId)) ** -1) / (5040 * (nu ** 7))) * (61 + (662 * ((math.tan(PHId)) ** 2)) + (1320 * ((math.tan(PHId)) ** 4)) + (720 * ((math.tan(PHId)) ** 6)))

    E_N_Long = (180 / Pi) * (RadLAM0 + (Et * X) - ((Et ** 3) * XI) + ((Et ** 5) * XII) - ((Et ** 7) * XIIA))
    return E_N_Long


def E_N_to_Lat_Long(North,East):
    Lat = E_N_to_Lat(East,DecimalPHI0,DecimalLAM0)
    Long = E_N_to_Long(East,DecimalLAM0)
    return [Lat,Long]

这一切运行良好，但它只有 +/- 3m 左右的转换精度，这对于我需要它的应用程序来说不够好。使用 Python 转换数千个坐标是否有更好的选择？

解决方法

这里有几个用于从 WGS84 转换为 BNG 的 Python 库，它们都利用 OSTN15 进行准确的转换： https://github.com/urschrei/convertbng

https://grid-banger.readthedocs.io/en/latest/README.html

可以在此处找到 OS 的更多信息：https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

相关推荐

导入项目后报错问题
依赖报错 idea导入项目后依赖报错，解决方案：https://blog.csdn.net/weixin_42420249/article/details/81191861 依赖版本报错：更换其他版本 无法下载依赖可参考：https://blog.csdn.net/weixin_42628809/a
idea不能识别yaml文件
使用mybatis plus常见错误
错误1：代码生成器依赖和mybatis依赖冲突 启动项目时报错如下 2021-12-03 13:33:33.927 ERROR 7228 [ main] o.s.b.d.LoggingFailureAnalysisReporter : *************************** APPL
gradle常见问题与错误
错误1：gradle项目控制台输出为乱码 # 解决方案：https://blog.csdn.net/weixin_43501566/article/details/112482302 # 在gradle-wrapper.properties 添加以下内容 org.gradle.jvmargs=-Df
Mybatis Plus传入参数0不起作用
错误还原：在查询的过程中，传入的workType为0时，该条件不起作用 <select id="xxx"> SELECT di.id, di.name, di.work_type, di.updated... <where> <if test=&qu
linux中make编译源码包失败
报错如下，gcc版本太低 ^ server.c:5346:31: 错误：‘struct redisServer’没有名为‘server_cpulist’的成员 redisSetCpuAffinity(server.server_cpulist); ^ server.c: 在函数‘hasActiveC
微服务启动错误：Command line is too long. Shorten command line for
解决方案1 1、改项目中.idea/workspace.xml配置文件，增加dynamic.classpath参数 2、搜索PropertiesComponent，添加如下 <property name="dynamic.classpath" value="tru
vue常见错误
删除根组件app.vue中的默认代码后报错：Module Error (from ./node_modules/eslint-loader/index.js): 解决方案：关闭ESlint代码检测，在项目根目录创建vue.config.js，在文件中添加 module.exports = { lin
spark常见错误
查看spark默认的python版本 [root@master day27]# pyspark /home/software/spark-2.3.4-bin-hadoop2.7/conf/spark-env.sh: line 2: /usr/local/hadoop/bin/hadoop: No s
matplotlib报错：AttributeError: module 'backend_interagg' has no attribute 'FigureCanvas'. Did you mean: 'FigureCanvasAgg'?
使用本地python环境可以成功执行 import pandas as pd import matplotlib.pyplot as plt # 设置字体 plt.rcParams['font.sans-serif'] = ['SimHei'] # 能正确显示负号 p
gitlab登录失败，报错：This challenge page was accidentally cached by an intermediary and is no longer available.
设置时间 控制面板
后端开发常见错误
错误1：Request method ‘DELETE‘ not supported 错误还原：controller层有一个接口，访问该接口时报错：Request method ‘DELETE‘ not supported 错误原因：没有接收到前端传入的参数，修改为如下 参考 错误2：cannot r
docker常见错误
错误1：启动docker镜像时报错：Error response from daemon: driver failed programming external connectivity on endpoint quirky_allen 解决方法：重启docker -> systemctl r
idea常见错误
错误1：private field ‘xxx‘ is never assigned 按Altʾnter快捷键，选择第2项 参考：https://blog.csdn.net/shi_hong_fei_hei/article/details/88814070 错误2：启动时报错，不能找到主启动类 #
pip安装依赖失败
报错如下，通过源不能下载，最后警告pip需升级版本 Requirement already satisfied: pip in c:\users\ychen\appdata\local\programs\python\python310\lib\site-packages (22.0.4) Coll
maven常见错误
错误1：maven打包报错 错误还原：使用maven打包项目时报错如下 [ERROR] Failed to execute goal org.apache.maven.plugins:maven-resources-plugin:3.2.0:resources (default-resources)
cloud项目中常见错误
错误1：服务调用时报错 服务消费者模块assess通过openFeign调用服务提供者模块hires 如下为服务提供者模块hires的控制层接口 @RestController @RequestMapping("/hires") public class FeignControl
springboot常见错误
错误1：运行项目后报如下错误 解决方案 报错2：Failed to execute goal org.apache.maven.plugins:maven-compiler-plugin:3.8.1:compile (default-compile) on project sb 解决方案：在pom.
springmvc拦截器中使用redisTemplate报空指针异常
参考 错误原因 过滤器或拦截器在生效时，redisTemplate还没有注入 解决方案：在注入容器时就生效 @Component //项目运行时就注入Spring容器 public class RedisBean { @Resource private RedisTemplate<String
vite报错：npm ERR! command C:WINDOWSsystem32cmd.exe /d /s /c C:UsersychenAppDataLocalTempnpx-dde1c848.cmd
使用vite构建项目报错 C:\Users\ychen\work>npm init @vitejs/app @vitejs/create-app is deprecated, use npm init vite instead C:\Users\ychen\AppData\Local\npm-
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot