Python 地理空间坐标格式转换
Python Geo-Spatial Coordinate Format Conversion
我有一个包含 6 列坐标对的数据框：度|分|秒（纬度和经度）。这称为 NAD83 格式。我想将它们转换成一个只有 2 列的十进制格式的新数据框，称为 NAD27。
我常用的库geopy几乎支持所有格式，所以实际上并没有专门的转换函数。我浏览了此处的文档以确保：
https://geopy.readthedocs.io/en/1.10.0/
python是否有任何其他方法可以转换为 NAD27？
感谢阅读
假设您的 DataFrame df 包含列 lonD、lonM、lonS、latD、latM 和 latS.
然后下面应该工作，在内部使用 geopandas、shapely 和 pyproj。
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
def dms_to_dec(d, m, s):
sign = 1 - 2 * np.signbit(d)
return d + sign * m / 60 + sign * s / 3600
points = df.apply(lambda row: Point(dms_to_dec(*row[['lonD', 'lonM', 'lonS']]),
dms_to_dec(*row[['latD', 'latM', 'latS']])),
axis=1)
gdf_nad83 = gpd.GeoDataFrame(df, geometry=points, crs={'init': 'EPSG:4269'})
gdf_nad27 = gdf_nad83.to_crs({'init': 'EPSG:4267'})
因为我 运行 也喜欢这个，并且发现 df.apply() 方法太慢，所以我转而使用 MultiPoint() 对象并使用向量化操作，然后将该单个对象转换为Point()s 与 list().
此外，我们需要考虑 DMS 列可能仅在 D 列上包含 - 符号这一事实。如果是这种情况，而且你很幸运，DataFrame 是使用 numpy 浮点数创建的，那么 "-0.0" 可能已存储为 numpy.NZERO（负零），在这种情况下，我们仍然可以使用 [=57 恢复符号=].否则，标志可能会丢失，赤道以南或 zero-th 子午线以西的点将显示为北或东。
明确一点：D、M、S 坐标表示法只是一种不同的纬度和经度坐标表示法，其中 D、M 和 S 代表 度、（弧）分钟和（弧）秒。小数是另一种，它将度值与弧分和弧秒组合成一个数字；弧分是 1/60 度，弧秒是 1/3600 度，因此您可以做一些数学运算将这些值加在一起（保留度数的符号）。 GeoPy 希望使用十进制值，因此您需要将 arc-seconds 和 arc-minutes 折叠成度值。
另一方面，NAD83 和 NAD27 不是 geodetic datums or geodetic systems，并且此类系统与符号无关。它们只是一种标准化的方式，用于指定要使用的坐标系以及坐标系锚定到的参考点。
也就是说，geopandas 可以用于在不同的大地基准之间形成运行。该项目分别接受 CRS strings to define what coordinate system to use when interpreting points (of which the geodetic datum is a component); using a coordinate system database such as https://spatialreference.org/ to find the EPSG codes for NAD83 and NAD27, gives us EPSG:4269 and EPSG:4267。请注意，您不必在此处创建数据框，如果您只需要 conversion.
，则 GeoSeries 就足够了
因此，鉴于您有度、分和秒，您需要将这些值转换为十进制坐标以提供给 geopandas。而且您希望快速高效。您可以通过使用向量化计算来做到这一点（其中 numpy 直接在数据的机器表示上使用非常快速的算术运算将计算应用于所有行，而不是 Python 表示）。
我在这里坚持相同的约定，输入是一个 Pandas DataFrame df，其中包含列 lonD、lonM、lonS， latD、latM 和 latS。使用 geopandas、numpy 和 shapely:
import geopandas as gpd
import numpy as np
from shapely.geometry import asMultiPoint
def vec_dms_to_dec(d, m, s):
"""convert d, m, s coordinates to decimals
Can be used as a vectorised operation on whole numpy arrays,
each array must have the same shape.
Handles signs only present on the D column, transparently.
Note that for -0d Mm Ss inputs, the sign might be have been lost!
However, if it was preserved as np.NZERO, this function will
recover it with np.signbit().
"""
assert d.shape == m.shape == s.shape
# account for signs only present on d
if (m >= 0).all() and (s >= 0).all():
# all s and m values are without signs
# so only d carries this info. Use the sign *bit* so negative
# and positive zero are distinguished correctly.
sign = np.where(np.signbit(d), np.ones_like(d) * -1.0, np.ones_like(d))
else:
sign = np.ones_like(d)
return d + sign * m / 60 + sign * s / 3600
# Generate the column names, grouped by component
comps = ([f"{c}{a}" for c in ("lon", "lat")] for a in 'DMS')
# Create a single MultiPoint object from the vectorised conversions of the
# longitude and latitude columns
mpoint = asMultiPoint(
vec_dms_to_dec(*(df[c].values for c in cols))
)
# Create a GeoSeries object from the MultiPoints object. Using `list()`
# produces `Point()` objects efficiently, faster than GeoSeries would
# otherwise.
# Interpret the points as using NAD83 == EPSG:4269
coords_nad83 = gpd.GeoSeries(list(mpoint), crs={'init': 'EPSG:4269'})
# Convert the series to NAD27 == EPSG:4267
coords_nad4267 = coords_nad83.to_crs(epsg=4267)
然后您可以自由地将它们再次转换回 D、M、S 表示法中的值：
from shapely.geometry import MultiPoint
def geoseries_to_dms(s, all_signed=True):
fractions, decimals = np.modf(np.array(MultiPoint(s.to_list())))
if not all_signed:
# only the d values signed. Looses information
# for input values in the open range (-1.0, 0.0)
fractions = np.abs(fractions)
fractions, minutes = np.modf(fractions * 60)
seconds = fractions * 60
return pd.DataFrame(
data=np.stack(
(decimals, minutes, seconds), axis=2
).reshape(-1, 6),
columns=loncols + latcols
)
上面用np.modf()从分数中拆分出小数部分，再将分数的绝对值拆分成弧分和弧秒。
如果您仍想使用 GeoDataFrame，请从转换后的 GeoSeries 创建一个，或者像创建 [=28] 一样从 MultiPoints() 对象创建一个=] 来自 MultiPoints() 对象，使用 GeoDataFrame(..., geometry=list(points), ...).
基准测试，或者为什么要矢量化
关于矢量化：我上面的代码将每个度、分和秒列作为三个独立的 numpy 数组，并使用这 3 个数组创建一个十进制度值数组，一步完成所有行。不需要单独调用纬度或经度值，因为 numpy 处理 d、m 和 s 作为数组，并且不关心它们是否只一维或 15.
这 运行 可以大大加快执行速度。为了对此进行基准测试，让我们创建一个具有任意数量的 dms 坐标的新数据框；我发现只生成十进制值并将其转换为 dms 值更容易：
import numpy as np
import pandas as pd
from shapely.geometry import Point, asMultiPoint
def random_world_coords(n):
coords = np.random.random((2, n))
coords[0] = coords[0] * 180 - 90 # lat between -90, 90
coords[1] = coords[1] * 360 - 180 # lon between -180, 180
# convert to d, m, s
fractions, decimals = np.modf(coords)
fractions, minutes = np.modf(fractions * 60)
seconds = fractions * 60
return pd.DataFrame(
data=np.stack((decimals, minutes, seconds), axis=2).reshape(-1, 6),
columns=["lonD", "lonM", "lonS", "latD", "latM", "latS"]
)
并定义将这些值转换为小数点的方法，以适合 GeoSeries() 消费，作为函数。我删除了符号处理，因为 运行dom 数据包含所有 dms 列上的符号，这也使得对标量和数组操作使用相同的转换函数变得微不足道：
def dms_to_dec(d, m, s):
"""convert d, m, s coordinates to decimals"""
return d + m / 60 + s / 3600
def martinvalgur_apply(df):
return df.apply(
lambda row: Point(
dms_to_dec(*row[['lonD', 'lonM', 'lonS']]),
dms_to_dec(*row[['latD', 'latM', 'latS']])
),
axis=1
)
def martijnpieters_vectorised(df):
comps = ([f"{c}{a}" for c in ("lon", "lat")] for a in 'DMS')
return list(asMultiPoint(
dms_to_dec(*(df[c].values for c in comps))
))
此时您可以使用 IPython 的 %timeit 或其他基准测试库来测试 on 的速度：
df100 = random_world_coords(100)
%timeit martinvalgur_apply(df100)
%timeit martijnpieters_vectorised(df100)
# 433 ms ± 15.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 96.2 ms ± 7.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
这是 100 个项目，矢量化速度大约快 4.5 倍。
如果将计数增加到 1000，差异就会变得更加明显：
df1000 = random_world_coords(1000)
%timeit martinvalgur_apply(df1000)
%timeit martijnpieters_vectorised(df1000)
# 4.31 s ± 111 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 35.7 ms ± 909 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
所以在 1000 行时，矢量化仍然只需要几毫秒，而且花费的时间 更少 因为我们现在正在针对更大的数据集进行优化，但使用 [=137] 所花费的时间=] df.apply() 在这 1000 行上已经膨胀到超过 4 秒。
（注意：我还 运行 对使用 DataFrame.copy() 创建的每个测试的输入进行了深度复制，以确保我没有获得 already-processed 数据的优势，但对于 100 -> 1000 行的情况，时间仍然下降，而不是上升。
non-vectorised 版本花费的时间与行数成正比，因此10k 行是可预测的：
df10k = random_world_coords(10_000)
%timeit martinvalgur_apply(df10k)
%timeit martijnpieters_vectorised(df10k)
# 44.1 s ± 1.1 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 331 ms ± 14.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
df.apply() 版本预计需要 44 秒，但我确实需要等待整整 5 分钟才能得到结果，因为 IPython 仍然是 运行 测试7次。
矢量化方法仅耗时 331 毫秒，因此我们可以测试仅 100 万行的版本：
df1m = random_world_coords(1_000_000)
%timeit martijnpieters_vectorised(df1m)
# 3.18 s ± 114 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
因此矢量化方法也可以线性扩展，但它从一个低得多的值开始。大部分时间用于从 MultiPoint() 对象创建 Point() 对象列表，geopandas 项目可以通过 supporting Shapely GeometrySequence objects natively.
改进这一点
我有一个包含 6 列坐标对的数据框：度|分|秒（纬度和经度）。这称为 NAD83 格式。我想将它们转换成一个只有 2 列的十进制格式的新数据框，称为 NAD27。
我常用的库geopy几乎支持所有格式，所以实际上并没有专门的转换函数。我浏览了此处的文档以确保： https://geopy.readthedocs.io/en/1.10.0/
python是否有任何其他方法可以转换为 NAD27？
感谢阅读
假设您的 DataFrame df 包含列 lonD、lonM、lonS、latD、latM 和 latS.
然后下面应该工作，在内部使用 geopandas、shapely 和 pyproj。
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
def dms_to_dec(d, m, s):
sign = 1 - 2 * np.signbit(d)
return d + sign * m / 60 + sign * s / 3600
points = df.apply(lambda row: Point(dms_to_dec(*row[['lonD', 'lonM', 'lonS']]),
dms_to_dec(*row[['latD', 'latM', 'latS']])),
axis=1)
gdf_nad83 = gpd.GeoDataFrame(df, geometry=points, crs={'init': 'EPSG:4269'})
gdf_nad27 = gdf_nad83.to_crs({'init': 'EPSG:4267'})
因为我 运行 也喜欢这个，并且发现 df.apply() 方法太慢，所以我转而使用 MultiPoint() 对象并使用向量化操作，然后将该单个对象转换为Point()s 与 list().
此外，我们需要考虑 DMS 列可能仅在 D 列上包含 - 符号这一事实。如果是这种情况，而且你很幸运，DataFrame 是使用 numpy 浮点数创建的，那么 "-0.0" 可能已存储为 numpy.NZERO（负零），在这种情况下，我们仍然可以使用 [=57 恢复符号=].否则，标志可能会丢失，赤道以南或 zero-th 子午线以西的点将显示为北或东。
明确一点：D、M、S 坐标表示法只是一种不同的纬度和经度坐标表示法，其中 D、M 和 S 代表 度、（弧）分钟和（弧）秒。小数是另一种，它将度值与弧分和弧秒组合成一个数字；弧分是 1/60 度，弧秒是 1/3600 度，因此您可以做一些数学运算将这些值加在一起（保留度数的符号）。 GeoPy 希望使用十进制值，因此您需要将 arc-seconds 和 arc-minutes 折叠成度值。
另一方面，NAD83 和 NAD27 不是 geodetic datums or geodetic systems，并且此类系统与符号无关。它们只是一种标准化的方式，用于指定要使用的坐标系以及坐标系锚定到的参考点。
也就是说，geopandas 可以用于在不同的大地基准之间形成运行。该项目分别接受 CRS strings to define what coordinate system to use when interpreting points (of which the geodetic datum is a component); using a coordinate system database such as https://spatialreference.org/ to find the EPSG codes for NAD83 and NAD27, gives us EPSG:4269 and EPSG:4267。请注意，您不必在此处创建数据框，如果您只需要 conversion.
GeoSeries 就足够了
因此，鉴于您有度、分和秒，您需要将这些值转换为十进制坐标以提供给 geopandas。而且您希望快速高效。您可以通过使用向量化计算来做到这一点（其中 numpy 直接在数据的机器表示上使用非常快速的算术运算将计算应用于所有行，而不是 Python 表示）。
我在这里坚持相同的约定，输入是一个 Pandas DataFrame df，其中包含列 lonD、lonM、lonS， latD、latM 和 latS。使用 geopandas、numpy 和 shapely:
import geopandas as gpd
import numpy as np
from shapely.geometry import asMultiPoint
def vec_dms_to_dec(d, m, s):
"""convert d, m, s coordinates to decimals
Can be used as a vectorised operation on whole numpy arrays,
each array must have the same shape.
Handles signs only present on the D column, transparently.
Note that for -0d Mm Ss inputs, the sign might be have been lost!
However, if it was preserved as np.NZERO, this function will
recover it with np.signbit().
"""
assert d.shape == m.shape == s.shape
# account for signs only present on d
if (m >= 0).all() and (s >= 0).all():
# all s and m values are without signs
# so only d carries this info. Use the sign *bit* so negative
# and positive zero are distinguished correctly.
sign = np.where(np.signbit(d), np.ones_like(d) * -1.0, np.ones_like(d))
else:
sign = np.ones_like(d)
return d + sign * m / 60 + sign * s / 3600
# Generate the column names, grouped by component
comps = ([f"{c}{a}" for c in ("lon", "lat")] for a in 'DMS')
# Create a single MultiPoint object from the vectorised conversions of the
# longitude and latitude columns
mpoint = asMultiPoint(
vec_dms_to_dec(*(df[c].values for c in cols))
)
# Create a GeoSeries object from the MultiPoints object. Using `list()`
# produces `Point()` objects efficiently, faster than GeoSeries would
# otherwise.
# Interpret the points as using NAD83 == EPSG:4269
coords_nad83 = gpd.GeoSeries(list(mpoint), crs={'init': 'EPSG:4269'})
# Convert the series to NAD27 == EPSG:4267
coords_nad4267 = coords_nad83.to_crs(epsg=4267)
然后您可以自由地将它们再次转换回 D、M、S 表示法中的值：
from shapely.geometry import MultiPoint
def geoseries_to_dms(s, all_signed=True):
fractions, decimals = np.modf(np.array(MultiPoint(s.to_list())))
if not all_signed:
# only the d values signed. Looses information
# for input values in the open range (-1.0, 0.0)
fractions = np.abs(fractions)
fractions, minutes = np.modf(fractions * 60)
seconds = fractions * 60
return pd.DataFrame(
data=np.stack(
(decimals, minutes, seconds), axis=2
).reshape(-1, 6),
columns=loncols + latcols
)
上面用np.modf()从分数中拆分出小数部分，再将分数的绝对值拆分成弧分和弧秒。
如果您仍想使用 GeoDataFrame，请从转换后的 GeoSeries 创建一个，或者像创建 [=28] 一样从 MultiPoints() 对象创建一个=] 来自 MultiPoints() 对象，使用 GeoDataFrame(..., geometry=list(points), ...).
基准测试，或者为什么要矢量化
关于矢量化：我上面的代码将每个度、分和秒列作为三个独立的 numpy 数组，并使用这 3 个数组创建一个十进制度值数组，一步完成所有行。不需要单独调用纬度或经度值，因为 numpy 处理 d、m 和 s 作为数组，并且不关心它们是否只一维或 15.
这 运行 可以大大加快执行速度。为了对此进行基准测试，让我们创建一个具有任意数量的 dms 坐标的新数据框；我发现只生成十进制值并将其转换为 dms 值更容易：
import numpy as np
import pandas as pd
from shapely.geometry import Point, asMultiPoint
def random_world_coords(n):
coords = np.random.random((2, n))
coords[0] = coords[0] * 180 - 90 # lat between -90, 90
coords[1] = coords[1] * 360 - 180 # lon between -180, 180
# convert to d, m, s
fractions, decimals = np.modf(coords)
fractions, minutes = np.modf(fractions * 60)
seconds = fractions * 60
return pd.DataFrame(
data=np.stack((decimals, minutes, seconds), axis=2).reshape(-1, 6),
columns=["lonD", "lonM", "lonS", "latD", "latM", "latS"]
)
并定义将这些值转换为小数点的方法，以适合 GeoSeries() 消费，作为函数。我删除了符号处理，因为 运行dom 数据包含所有 dms 列上的符号，这也使得对标量和数组操作使用相同的转换函数变得微不足道：
def dms_to_dec(d, m, s):
"""convert d, m, s coordinates to decimals"""
return d + m / 60 + s / 3600
def martinvalgur_apply(df):
return df.apply(
lambda row: Point(
dms_to_dec(*row[['lonD', 'lonM', 'lonS']]),
dms_to_dec(*row[['latD', 'latM', 'latS']])
),
axis=1
)
def martijnpieters_vectorised(df):
comps = ([f"{c}{a}" for c in ("lon", "lat")] for a in 'DMS')
return list(asMultiPoint(
dms_to_dec(*(df[c].values for c in comps))
))
此时您可以使用 IPython 的 %timeit 或其他基准测试库来测试 on 的速度：
df100 = random_world_coords(100)
%timeit martinvalgur_apply(df100)
%timeit martijnpieters_vectorised(df100)
# 433 ms ± 15.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 96.2 ms ± 7.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
这是 100 个项目，矢量化速度大约快 4.5 倍。
如果将计数增加到 1000，差异就会变得更加明显：
df1000 = random_world_coords(1000)
%timeit martinvalgur_apply(df1000)
%timeit martijnpieters_vectorised(df1000)
# 4.31 s ± 111 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 35.7 ms ± 909 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
所以在 1000 行时，矢量化仍然只需要几毫秒，而且花费的时间 更少 因为我们现在正在针对更大的数据集进行优化，但使用 [=137] 所花费的时间=] df.apply() 在这 1000 行上已经膨胀到超过 4 秒。
（注意：我还 运行 对使用 DataFrame.copy() 创建的每个测试的输入进行了深度复制，以确保我没有获得 already-processed 数据的优势，但对于 100 -> 1000 行的情况，时间仍然下降，而不是上升。
non-vectorised 版本花费的时间与行数成正比，因此10k 行是可预测的：
df10k = random_world_coords(10_000)
%timeit martinvalgur_apply(df10k)
%timeit martijnpieters_vectorised(df10k)
# 44.1 s ± 1.1 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
# 331 ms ± 14.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
df.apply() 版本预计需要 44 秒，但我确实需要等待整整 5 分钟才能得到结果，因为 IPython 仍然是 运行 测试7次。
矢量化方法仅耗时 331 毫秒，因此我们可以测试仅 100 万行的版本：
df1m = random_world_coords(1_000_000)
%timeit martijnpieters_vectorised(df1m)
# 3.18 s ± 114 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
因此矢量化方法也可以线性扩展，但它从一个低得多的值开始。大部分时间用于从 MultiPoint() 对象创建 Point() 对象列表，geopandas 项目可以通过 supporting Shapely GeometrySequence objects natively.