python - 将本地坐标投影到带参考点的全球 GPS

Question

我在 geopandas GeoDataFrame 中有一堆形状（例如shapely LineStrings 或Polygons）。形状指定本地 200x200 米网格中的坐标，即所有坐标都在 (0, 0) 和 (200, 200) 之间。

我现在想在全球范围内“放置”这些行。为此，我想指定一个 GPS 点（具有给定的纬度/经度）作为参考。

我的第一个（天真的）方法是使用 geolib，获取所有形状的坐标（在本地 X/Y 中）并应用以下转换并“重新创建”形状：

# Convert coordinates to GPS location
from shapely.geometry import LineString
from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84  # the base geodesic (i.e. the world)

origin = (48.853772345870176, 2.350983211585546)  # this is somewhere in Paris, for example

def local_to_latlong(x, y, orientation=0, scale=1):
    """ Two step process.
    - First walk x meters to east from origin.
    - Then, from that point, walk y meters north from origin.
    
    Optional: 
    - orientation allows to "spin" the coordinates
    - scale allows to grow/shrink the distances
    """
    go_X = geod.Direct(*origin, orientation + 90, x * scale)  # x is East-coordinate
    go_Y = geod.Direct(go_X["lat2"], go_X["lon2"], orientation + 0, y * scale)  # y is North-coordinate
    return go_Y["lat2"], go_Y["lon2"]


original_line = LineString([(0,0), (100,100), (200,100)])
global_line = LineString([local_to_latlong(x, y) for y, x in original_line.coords])

但是，我希望这不是最聪明的方法，并且有更聪明的方法......

我想将这种转换应用于 GeoDataFrame 中的任何形状。理想情况下，它可以使用“to_crs”，但我不确定如何转换形状，以便它们“参考原点”以及使用哪个 crs。

Answer 1

• 鉴于您的来源是 EPSG:4326，您可以估计 UTM 区域
• 有了这个，你可以得到原点的 UTM 区域坐标
• 将您的自定义 200x200 米区域转换为 UTM 区域的坐标
• 最后to_crs()用来转换成EPSG:4326

import shapely.geometry
import geopandas as gpd
import pandas as pd
import numpy as np

# generate some polygons (squares), where grid is 200*200
gdf = gpd.GeoDataFrame(
    geometry=pd.DataFrame(
        np.repeat(np.sort(np.random.randint(0, 200, [20, 2]), axis=1), 2, axis=1)
    ).apply(lambda d: shapely.geometry.box(*d), axis=1)
)
# chage to linestrings, clearer when we plot
gdf["geometry"] = gdf["geometry"].exterior

origin = (2.350983211585546, 48.853772345870176)  # this is somewhere in Paris, for example

# work out utm crs of point.  utm is in metres
gdf_o = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(origin)], crs="EPSG:4326")
crs = gdf_o.estimate_utm_crs()
# where is origin in utm zone
xo,yo = gdf_o.to_crs(crs).loc[0,"geometry"].xy

# translate custom zone to co-ordinates of utm zone
# assume point is center of 200x200 grid (hence subtract 100)
gdf_gps = gdf["geometry"].translate(xoff=xo[0]-100, yoff=yo[0]-100).set_crs(crs).to_crs("epsg:4326")

# plot on map to show it has worked...
m = gdf_gps.explore()
m = gdf_o.explore(m=m, color="red", marker_kwds={"radius":20})
m