Python创建旋转地理坐标网格的问题排查与优化方案问询

问题背景 #

需要为项目创建等间距地理坐标网格，后续导入QGIS使用。最初基于4个参考坐标（定义左上角及X、Y方向间距）编写Python代码生成扩展网格，但出现两个问题：

• X、Y方向点排列不整齐，呈现抖动状
• 点间距离不符合要求

推测问题1源于外推方式导致的精度丢失，问题2源于坐标集定义不当。因未找到基于坐标定义旋转网格的方法才采用当前方案，现希望基于「左上角起始坐标、水平/垂直间距、行列数」这些输入，寻求解决办法或理想实现方案。

原实现代码 #

from string import ascii_lowercase
from geojson import Point, Feature, FeatureCollection, dump
import json  # 补充原代码缺失的导入

def getExtraPolatedPoint(p1,p2, ratio):
    'Creates a line extrapoled in p1->p2 direction'
    b = (p1[0]+ratio*(p2[0]-p1[0]), p1[1]+ratio*(p2[1]-p1[1]) )
    return b

L = [letter1+letter2 for letter1 in ascii_lowercase for letter2 in ascii_lowercase]

with open(r"4points.geojson") as json_file:
    res = json.load(json_file)
    
features = []
feature_collection = FeatureCollection(features)

tl = None

for a in res["features"]:
    b = a["geometry"]["coordinates"]
    if tl is not None:
        if tl[1]<b[1]:tl=b
        if tr[0]<b[0]:tr=b
        if br[1]>b[1]:br=b
        if bl[0]>b[0]:bl=b
    else:
        tl = tr = br = bl = b

refP1 = tl
refP2x = tr
refP2y = bl
refP2xy = br

for i in range(16):
    newP1  = getExtraPolatedPoint(refP1,refP2x,i)
    newP2xy = getExtraPolatedPoint(refP2y,refP2xy,i)
    
    for j in range(75):
        gridPoint =  getExtraPolatedPoint(newP1,newP2xy,j)
        point = Point(gridPoint)
        features.append(Feature(geometry=point, properties={"ID":L[i]+'{:02}'.format(j+1) }))

问题分析 #

• 原方案直接对经纬度进行线性外推，而经纬度属于球面坐标，线性计算会忽略球面曲率，导致精度丢失，出现点排列抖动的问题
• 手动通过四个参考点判断网格方向的逻辑存在误差，容易因坐标顺序或判断条件错误，导致点间间距不符合要求

解决方案 #

使用pyproj库处理球面坐标的精确计算，核心思路是基于方位角和球面距离生成网格，步骤如下：

1. 基于左上角起始点和参考点计算水平、垂直方向的方位角
2. 利用fwd函数根据方位角、间距、行列数，计算每行的起始点和终点
3. 通过npts函数在起始点与终点之间生成等间距的点，确保网格整齐且间距准确

核心实现代码 #

from pyproj import Geod
from geojson import Point, Feature, FeatureCollection, dump
from string import ascii_lowercase

# 初始化地理坐标系（WGS84）
geod = Geod(ellps="WGS84")

# 输入参数
start_lon, start_lat = tl[0], tl[1]  # 左上角起始坐标
x_spacing = 1000  # 水平间距（单位：米）
y_spacing = 1000  # 垂直间距（单位：米）
rows = 16  # 行数
cols = 75  # 列数

# 计算水平、垂直方向的方位角
_, az_x, _ = geod.inv(start_lon, start_lat, refP2x[0], refP2x[1])
_, az_y, _ = geod.inv(start_lon, start_lat, refP2y[0], refP2y[1])

# 生成网格点集合
features = []
id_prefixes = [letter1+letter2 for letter1 in ascii_lowercase for letter2 in ascii_lowercase]

for row_idx in range(rows):
    # 计算当前行的起始点（沿垂直方向移动row_idx个间距）
    row_start_lon, row_start_lat, _ = geod.fwd(start_lon, start_lat, az_y, y_spacing * row_idx)
    # 计算当前行的终点（沿水平方向移动cols个间距）
    row_end_lon, row_end_lat, _ = geod.fwd(row_start_lon, row_start_lat, az_x, x_spacing * cols)
    # 生成当前行的所有等间距点
    row_points = geod.npts(row_start_lon, row_start_lat, row_end_lon, row_end_lat, cols-1)
    # 插入行起始点到列表头部
    row_points.insert(0, (row_start_lon, row_start_lat))
    
    # 为每个点添加属性并构造Feature
    for col_idx, (lon, lat) in enumerate(row_points):
        point = Point((lon, lat))
        features.append(Feature(
            geometry=point,
            properties={"ID": f"{id_prefixes[row_idx]}{col_idx+1:02}"}
        ))

# 保存为GeoJSON文件
with open("rotated_grid.geojson", "w") as f:
    dump(FeatureCollection(features), f)

更新说明 #

• 采用经纬度地理坐标系，旋转需求源于需贴合地图上的特定对象
• 最终借助pyproj库，利用方位角（azimuth）、fwd函数生成最大坐标，再通过npts函数生成指定间距的点，解决了问题

内容的提问来源于stack exchange，提问作者significant