用GeoPandas比较国家地区大小
原文：
towardsdatascience.com/comparing-country-sizes-with-geopandas-2ce027282ba0
快速成功数据科学
我住在德克萨斯州（耶-哈！），我们喜欢将德克萨斯州与其他州、国家、岛屿、小行星等进行比较。我们甚至会将德克萨斯州与德克萨斯州进行比较：
没错，伙计们，德克萨斯州比德克萨斯州还要大！（作者）
我喜欢比较一个地理区域与另一个地理区域大小的梗图。尽管我对地理有很好的了解，但我经常感到惊讶。例如，澳大利亚是否比巴西大？新西兰是否比意大利小？德克萨斯州与德国、印度和智利相比如何？以下是答案：
巴西略大于澳大利亚（作者）
意大利略大于新西兰（作者）
德克萨斯州比德国大（作者）
他们不是没有理由称印度为次大陆的！（作者）
智利比德克萨斯州更大（更长）（作者）
除了娱乐性，这些比较也是学生们的优秀学习工具。
在这个快速成功数据科学项目中，我们将编写 Python 代码，以便轻松叠加和比较各国或美国各州的地图。这是一种了解地理空间数据集及其如何使用*GeoPandas*库进行操作的好方法。
流程
这里是我们将使用的一般流程：
-
加载每个国家或州的shapefiles（边界多边形）。
-
使用坐标参考系统投影边界。
-
将一个边界的质心平移以匹配另一个。
-
将一个边界平移和旋转，直到它很好地适应另一个。
-
绘制边界。
GeoPandas 库
*GeoPandas*是一个开源的第三方库，旨在支持 Python 中的地理空间映射。它建立在 pandas 数据分析库之上，使得处理地理空间矢量数据与处理表格数据类似。它还使 Python 能够执行通常需要专用地理空间数据库的操作，例如 PostGIS。
GeoPandas 的GeoDataFrame是一个具有特殊“geometry”列的 pandas DataFrame，用于存储位置数据。这个列将几何对象的类型（如点、线字符串、多边形等）和绘制它所需的坐标（经度和纬度）捆绑在一起。
GeoDataFrame 中的几何列（由作者Python Tools for Scientists提供）
您可以在这里找到 GeoPandas 的安装说明。像往常一样，conda 用户操作简单（conda install geopandas）。如果您使用 pip，您还需要安装 shapely、pyproj、fiona、pyogrio 和 rtree 的主要依赖项，具体请参阅说明。
我们还需要 Matplotlib 进行绘图。作为 GeoPandas 的可选依赖项，当采用 conda 路径时，它应该会自动安装。
Shapefiles
对于国家和州边界，我们将依赖shapefiles。shapefile 是一种地理空间矢量数据格式，用于地理信息系统（GIS）软件。它用于在空间上描述矢量特征，如线、点和多边形。文件中的项目包括描述性属性，如名称。
向量数据类型的示例（来自Python Tools for Scientists）
“shapefile”不是一个单独的文件，而是一个集合，位于文件夹中的文件。下面的图示显示了用于主要湖泊轮廓的示例 shapefile。注意这个文件夹中的文件名如何使用文件夹名作为前缀。
示例 shapefile 目录（图片由作者提供）
三个文件是强制性的，文件扩展名为.shp、.shx和.dbf。虽然实际的 shapefile 专门关联到.shp文件，但其他支持文件是捕获形状几何所必需的。有关这些文件的更多详细信息，请查看维基百科页面上的 shapefiles 和这篇文章：
我们需要向 GeoPandas 提供 shapefiles文件夹的目录路径（压缩或解压缩）。GeoPandas 将使用此文件夹中的文件自动在 GeoDataFrame 中生成geometry列。还有什么比这更简单吗？
国家比较代码
以下代码通过在大约相同的比例上叠加地图来比较两个国家的地图。在文章的后面部分，我们将探讨比较美国州。
导入库
GeoPandas，与 pandas 类似，是基于 Matplotlib 构建的，因此我们将导入 Matplotlib 进行绘图。
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import geopandas as gpd
分配常量
为了方便起见，我们将在程序的顶部分配常数，这样它们将很容易找到和更改。这些包括国家名称、旋转角度（以便更好地“适应”一个国家到另一个国家），以及从 3D 地球到 2D 地图的投影坐标系（CRS）。
# Assign constants:
COUNTRY1 = 'China' # Plots in the background (back).
COUNTRY2 = 'Australia' # Plots in the foreground (front).
ROTATE = -17 # Direction and amount to rotate COUNTRY2
CRS_1 = 'EPSG:3415' # Coordinate Reference System for COUNTRY1
CRS_2 = 'EPSG:3112' # Coordinate Reference System for COUNTRY2
坐标系（CRS）对于准确比较一个国家的面积与另一个国家非常重要。你可能熟悉标准 _墨卡托*地图投影以及它们如何扭曲极地附近国家的面积：
世界墨卡托投影（作者提供）
格陵兰与南美洲在大小上不相上下，而南极洲则纯粹是夸张。这是由于经线在极地附近汇聚的结果。因为这些线是非平行的，将它们展平到 2D 地图上会在赤道以外造成越来越大的扭曲：
经线在极地附近汇聚（来自作者提供的Real-world Python）
这种扭曲可能会严重破坏大小比较。当你将澳大利亚与加拿大的墨卡托投影和更好地保留面积的兰伯特投影进行比较时，会发生以下情况：
澳大利亚与加拿大的墨卡托投影（左）与兰伯特投影（右）（作者提供）
幸运的是，有很多方法可以将球体投影到平面上，其中一些方法可以保留面积。不幸的是，这需要付出代价。
地图投影就像挤压气球。如果你掌握了一个参数，比如面积，另一个参数，比如角度或形状，就会在你的手指之间凸起。
虽然有“等面积”投影，如 _兰伯特方位等面积投影或阿尔伯斯等面积圆锥投影，它们可以保留地图上区域的相对大小*，但它们会扭曲国家轮廓。如果人们注意到国家的形状不熟悉，他们会怀疑显示的真实性，我们可不能有这种情况！
我找到的最好解决方案是使用每个国家的推荐本地投影系统进行投影，然后将一个移动并绘制在另一个之上。本地投影最适合在本地保留面积和形状。缺点是必须处理许多本地投影。尽管如此，这比使用单一投影对两国来说更令人满意，尤其是如果它们位于不同的纬度范围内。
我们将使用 EPSG 代码选择投影。根据 维基百科，“EPSG 地理参数数据集（也称为 EPSG 注册表）是一个公共注册表，包含 地理基准、空间参考系统、地球椭球体、坐标转换和相关的 测量单位，由欧洲石油调查组（EPSG）的成员于 1985 年创建。”
要找到国家的 EPSG 代码，请使用诸如 “澳大利亚推荐的 EPGS 代码是什么？” 等查询在线搜索。以下是我在这个项目中使用的代码列表：
本文使用的 ESPG 代码（作者提供）
如果您想使用长列表的代码，考虑将它们作为字典输入，国家名称作为键，代码作为值。
选择和投影形状文件
Geopandas 以前附带 Natural Earth 公有领域 数据集，其中包含国家和州边界形状文件。这已在 2024 年初被弃用，因此我们必须手动下载文件。
幸运的是，这很简单。只需导航到这个 网站 并点击下面高亮显示的黄色下载链接：
点击高亮链接下载 “countries” 形状文件（作者提供）
将压缩文件移动到包含您的 Python 脚本或笔记本的文件夹中，然后运行以下代码：
# Download the shapefile used below from:
# https://www.naturalearthdata.com/downloads/110m-cultural-vectors/
world = gpd.read_file('ne_110m_admin_0_countries.zip')
# Select countries from GeoDataFrames:
front = world[world['SOVEREIGNT'] == COUNTRY2]
back = world[world['SOVEREIGNT'] == COUNTRY1]
# Project to a CRS before calculating the centroid:
front = front.to_crs(CRS_2)
back = back.to_crs(CRS_1)
形状文件被分配给 world 变量，并使用 GeoPandas 的 read_file() 方法加载。这里显示不出来，但如果你想查看内容，输入 world.head()，就像查看 pandas DataFrame 一样。
front 国家将在显示的前景中以红色显示。back 国家将在背景中以黑白显示。
最后两行将形状文件转换为由 EPSG 代码指定的正确坐标参考系统。
旋转和移动边界多边形
以下代码找到每个国家多边形的 质心，如果需要，旋转 front（前景）多边形，并将其移动，使两个多边形具有相同的中心点。
# Calculate the geographical center of each country:
front_center = front.geometry.centroid.iloc[0]
back_center = back.geometry.centroid.iloc[0]
# Rotate COUNTRY2 around center by ROTATE value (0 = no rotation):
front_rotated = front.set_geometry(front.rotate(ROTATE, origin='centroid'))
# Shift foreground country to overlap background country:
front_shifted = front_rotated.set_geometry(front_rotated.translate(
xoff=back_center.x-front_center.x,
yoff=back_center.y-front_center.y))
ROTATE 常量允许您旋转前景国家的多边形，以改善其在背景多边形内的放置。值应在 0 到 180 之间。负值将使多边形 顺时针旋转。
还可以将前景多边形在垂直和水平方向上移动。我们将在关于美国各州的章节中稍后介绍这个选项。
绘制地图并检查面积
下一个单元格绘制地图，然后打印每个国家多边形的面积，作为质量控制步骤。后者有助于确保我们使用的是有效的 EPSG 代码。
# Plot the background country (COUNTRY1):
base = back.plot(color='white', edgecolor='black')
# Plot the foreground country (COUNTRY2):
front_shifted.plot(ax=base, color='red', alpha=0.5)
# Add a legend:
red_patch = mpatches.Patch(color='red', alpha=0.5, label=COUNTRY2)
plt.legend(handles=[red_patch], loc='best')
# Turn off axis labels and ticks:
plt.xticks([])
plt.yticks([])
# Add a title:
plt.title(f"{COUNTRY2} vs. {COUNTRY1} Size Comparison")
# Save plot (optional):
plt.savefig(f'{COUNTRY2} vs {COUNTRY1}.png', dpi=600)
plt.show();
# Print the country area beneath the map:
front['area_km2'] = front.area / 10**6
back['area_km2'] = back.area / 10**6
print(f"{COUNTRY2} area = {front.area_km2.to_string(index=False)} sq km")
print(f"{COUNTRY1} area = {back.area_km2.to_string(index=False)} sq km n")
前两个步骤绘制背景和前景边界。然后我们添加图例。这里的一个关键步骤是使用 Matplotlib 的mpatches.Patch()类制作一个红色矩形来指定前景国家。
然后我们关闭 x 轴和 y 轴的刻度，添加标题，并保存绘图。
最后，我们计算每个多边形的面积并打印结果。然后我们可以在网上搜索每个国家的面积，看看结果是否匹配。由于我们使用的边界多边形是简化的，所以它们不会完全精确，但它们足够接近，以确保我们的 EPSG 代码没有出错。
下面是澳大利亚与中国比较的结果：
澳大利亚与中国（作者）
Australia area = 7.593688e+06 sq km
China area = 1.028958e+07 sq km
关于绘制联合王国的说明
我们的 shapefile 包括与联合王国的福克兰群岛，这可能会弄乱您的质心计算。为了解决这个问题，请在“ADMIN”列下删除“福克兰群岛”，如下所示：
# Remove the Falkland Islands:
world = gpd.read_file('ne_110m_admin_0_countries.zip')
world = world[(world.ADMIN != 'Falkland Islands')]
接下来，我们将比较美国各州与其他州和国家。
美国各州比较代码
我们必须从 Natural Earth 网站下载一个新的 shapefile 来比较美国各州。点击下面突出显示的黄色链接，并将文件移动到包含您的 Python 脚本或笔记本的文件夹中：
点击突出显示的链接下载“州和省”shapefile（作者）
记住，您不需要解压缩文件即可使用它。
各州之间比较
下面是用于比较一个州与另一个州（在这种情况下，加利福尼亚州与德克萨斯州）的完整代码。它与之前基本相同，但使用了新的 shapefile 和一些变量名称的调整：
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import geopandas as gpd
# Assign constants:
STATE1 = 'Texas' # Plots in the background (back).
STATE2 = 'California' # Plots in the foreground (front).
ROTATE = 13 # Direction and amount to rotate STATE2
CRS1 = 'EPSG:5070' # Coordinate Reference System for STATE1
CRS2 = 'EPSG:5070' # Coordinate Reference System for STATE2
# Load the shapefile:
usa = gpd.read_file('ne_110m_admin_1_states_provinces.zip')
# Select states from GeoDataFrame:
front = usa[(usa.name == STATE2)]
back = usa[(usa.name == STATE1)]
# Project to a CRS before calculating the centroid:
front = front.to_crs(CRS2)
back = back.to_crs(CRS1)
# Calculate the geographical center of each state:
front_center = front.geometry.centroid.iloc[0]
back_center = back.geometry.centroid.iloc[0]
# Rotate STATE2 around center by ROTATE value (0 = no rotation):
front_rotated = front.set_geometry(front.rotate(ROTATE, origin='centroid'))
# Shift foreground state to overlap background state:
front_shifted = front_rotated.set_geometry(front_rotated.translate(
xoff=back_center.x-front_center.x,
yoff=back_center.y-front_center.y))
# Plot the background state (STATE1):
base = back.plot(color='white', edgecolor='k', lw=3)
# Plot the foreground state (STATE2):
front_shifted.plot(ax=base, color='red', alpha=0.5)
# Add a legend:
red_patch = mpatches.Patch(color='red', alpha=0.5, label=STATE2)
plt.legend(handles=[red_patch], loc='best')
# Turn off axis labels and ticks:
plt.xticks([])
plt.yticks([])
# Add a title:
plt.title(f"{STATE2} vs. {STATE1} Size Comparison")
plt.savefig(f'{STATE2}_vs_{STATE1}.png', dpi=600)
plt.show();
# Print the area of each beneath the map:
front['area_km2'] = front.area / 10**6
back['area_km2'] = back.area / 10**6
print(f"{STATE2} area = {front.area_km2.to_string(index=False)} sq km")
print(f"{STATE1} area = {back.area_km2.to_string(index=False)} sq km n")
加利福尼亚州与德克萨斯州（作者）
各州与国家比较
要将一个州（或多个州）与另一个国家进行比较，我们必须加载国家和州的 shapefile。以下是一个比较澳大利亚与下 48 个州的示例：
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import geopandas as gpd
# Assign constants:
COUNTRY1 = 'USA' # Plots in the background (back).
COUNTRY2 = 'Australia' # Plots in the foreground (front).
ROTATE = 0 # Direction and amount to rotate COUNTRY2
CRS_1 = 'EPSG:5070' # Coordinate Reference System for COUNTRY1
CRS_2 = 'EPSG:3112' # Coordinate Reference System for COUNTRY2
# Load the shapefiles:
world = gpd.read_file('ne_110m_admin_0_countries.zip')
usa = gpd.read_file('ne_110m_admin_1_states_provinces.zip')
# Exclude Alaska and Hawaii from USA:
usa = usa[(usa.name != 'Alaska') & (usa.name != 'Hawaii')]
# Select regions from GeoDataFrames:
front = world[world['SOVEREIGNT'] == COUNTRY2]
back = usa
# Project to a CRS before calculating the centroid:
front = front.to_crs(CRS_2)
back = back.to_crs(CRS_1)
# Calculate the geographical center of each region:
front_center = front.geometry.centroid.iloc[0]
back_center = back.geometry.centroid.iloc[0]
# Rotate front around center by ROTATE value (0 = no rotation):
front_rotated = front.set_geometry(front.rotate(ROTATE, origin='centroid'))
# Shift front to overlap back:
front_shifted = front_rotated.set_geometry(front_rotated.translate(
xoff=back_center.x-front_center.x,
yoff=back_center.y-front_center.y))
# Plot the background region:
base = back.plot(color='white', edgecolor='black')
# Plot the foreground region:
front_shifted.plot(ax=base, color='red', alpha=0.5)
# Add a legend:
red_patch = mpatches.Patch(color='red', alpha=0.5, label=COUNTRY2)
plt.legend(handles=[red_patch], loc='best')
# Turn off axis labels and ticks:
plt.xticks([])
plt.yticks([])
# Add a title:
plt.title(f"{COUNTRY2} vs. {COUNTRY1} Size Comparison")
plt.show();
澳大利亚与下 48 个州（作者）
这里有一个问题。澳大利亚太北了，不适合进行比较。为了解决这个问题，我们需要将澳大利亚向下移动。
垂直和水平移动边界多边形
我们已经通过将前景多边形的质心设置为背景多边形的质心来移动前景多边形。为了进一步在南北和东西方向上调整它，我们可以分配额外的常数：
SHIFT_SOUTH = 550_000 # in meters
SHIFT_WEST = 350_000 # in meters
与ROTATE常数一样，显示的值是通过试错得到的。
我们接着编辑将前多边形移动以包含新常数的代码行：
front_shifted = front_rotated.set_geometry(front_rotated.translate(
xoff=back_center.x-front_center.x-SHIFT_WEST,
yoff=back_center.y-front_center.y-SHIFT_SOUTH))
这将澳大利亚向下移动 550 公里，向左移动 350 公里：
澳大利亚与下 48 个州对比，澳大利亚向西和向南移动（作者绘制）
这样可以做出更好的比较。
有趣的是，我在 20 世纪 90 年代初住在墨尔本时，并没有意识到澳大利亚的面积有美国大陆那么大。然后我从墨尔本东南部飞往西海岸的珀斯。这次飞行略超过 4 小时，与从亚特兰大到圣地亚哥的旅行相似。一些简单的数学计算（4 小时 x 500 英里/小时）显示我飞行了 2000 多英里。难怪澳大利亚被归类为大陆！
摘要
Python 的GeoPandas库使得操作地理空间数据以创建有趣的地图变得容易。GeoPandas 的表格GeoDataFrame在流行的 pandas DataFrame的基础上构建，通过包含一个“Geometry”列来存储如多边形和线这样的矢量数据的坐标。这使得 GeoDataFrames“空间感知”。
“Geometry”列通常使用一种名为shapefile的地理空间矢量数据格式来填充。使用 GeoDataFrame，你可以使用坐标参考系统将地理空间数据投影到平面上，旋转和移动数据，计算面积，绘制数据等等。
谢谢！
感谢阅读，并请关注我，未来将有更多快速成功数据科学项目。
扫一扫
1510
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
