使用Python构建尼泊尔的地形图

Python构建尼泊尔地形图

介绍

曾经想过你的国家的地形如何影响经济和政治发展吗？地形图-使用等高线进行可视化的地球表面地图-可以帮助回答这些问题！我们将使用Python创建尼泊尔的地形图，尼泊尔是一个具有有趣地形环境的国家。您将学习如何读取描述国家地形的地理空间数据，如何解释这些数据，以及如何进行可视化。生成的地图可以与其他感兴趣的数据结合在一起，以了解国家地形如何影响其经济和/或政治发展的极低的子国家层次。本博客文章将教您如何生成一个非常有趣的工具，可以为政策和私营部门发展提供信息！

学习目标

• 掌握数字高程数据的数据分析技术。
• 学习如何在Python中使用地理空间数据和相关分析工具。
• 掌握制图技术的知识。
• 发展有效的数据可视化沟通技巧。
• 了解高程对不平等和贫困的重要性。

本文作为数据科学博文马拉松的一部分发表。

什么是地形图？

地形图是地球表面的地图，使用等高线进行可视化。地形图是一种宝贵的工具，用于在陌生地形中导航，并作为城市规划和灾害管理的输入。它们经常用于了解政策或私营部门项目在基础设施发展周围的空间背景，以确定易受自然灾害影响或面临有限教育、医疗保健和基础设施等基本服务的地区，或用于自然资源管理。最终，这些地图可以作为基于证据的决策的输入。在本博客文章中，我们将使用Python为尼泊尔创建一个具有非常有趣地形环境的地形图。

数据描述

为了生成我们的地图，我们将依赖美国地质调查局（USGS）发布的数据。USGS是美国联邦政府的科学机构，生成有关自然资源、地质学、地理学、水资源和自然灾害的数据和研究。要访问他们的数据页面，请在Google中键入“USGS数据”或点击指向他们的Earth Explorer的链接。Earth Explorer是一个在线工具和数据门户，允许您搜索、访问和下载各种地球科学数据。您必须设置一个帐户并登录以充分使用数据。

• 技术专家开始怀疑ChatGPT，人工智能的“幻觉”是否会消失
• 对生成式人工智能负责任使用的紧迫需求
• “AI模型可以帮助确定患者的癌症起源”

数据下载

由于尼泊尔具有独特的地形特征，本博客文章将以尼泊尔为例。尼泊尔拥有世界上最具挑战性和有趣的地形之一。尼泊尔拥有14座海拔超过8000米的山峰中的8座（Trekking Trail Nepal），该国分为三个非常不同的地形区域：山脉、山丘和特拉伊（或平原） (DHS)。虽然这些特征使该国独特和有趣，但一些研究表明尼泊尔的地形使其难以连接该国、向其人口提供基本服务，并对可持续发展道路施加风险和障碍。

为此，我们将在搜索条件中过滤尼泊尔，如下图所示。选择了尼泊尔后，我们选择了我们感兴趣的数据集。要这样做，请单击数据集选项卡，选择数字高程。数字高程数据有几个选项，虽然您可以使用其中几个数据集，但我们将使用全球多分辨率地形高程数据2010（GMTED2010）数据。该数据提供了地球表面的全球覆盖，具有多种分辨率（从7.5弧秒（约250米）到30弧秒（约1千米））。它是由航天器和飞机遥感数据生成的，包括卫星测高、立体影像和地形图。

选择数据后，点击结果选项卡。您现在可以通过单击带有下载选项的符号来下载数据。您还可以通过足迹图标显示数据。我们以最高分辨率（7.5弧秒）下载数据。重要的是，为了覆盖尼泊尔的所有区域，我们需要下载两个不同的马赛克（部分）的底层数据，并稍后将它们合并。您将看到生成的数据集以tif格式呈现，这表示栅格数据。

使用Python的地理空间分析工具准备数据

Python提供了几个用于地理空间分析的工具。在本博客文章中，我们将依赖于Rasterio库，该库可以读取和写入地理空间栅格数据（网格数据）。让我们开始并读取我们之前下载到Jupyter Notebook中的第一个数据镶嵌（部分）：

#导入相关库（在安装后）
import rasterio
import matplotlib.pyplot as plt
import numpy as np

#读取数据并显示数据集的形状
file = rasterio.open(r'path\10n060e_20101117_gmted_mea075.tif')
dataset = file.read()
print(dataset.shape)

我们还要上传第二个数据镶嵌并将它们合并。为此，我们按照Python中的标准栅格数据读取和操作技术进行操作，如下所示：

#上传第二个数据集并显示数据集的形状
file2 = rasterio.open(r'path\30n060e_20101117_gmted_mea075.tif')
dataset2 = file2.read()
print(dataset2.shape)


#合并两个数据集
from rasterio.merge import merge
from rasterio.plot import show

#创建空列表
src_files_to_mosaic = []

#将列表添加到两个文件
src_files_to_mosaic.append(file)
src_files_to_mosaic.append(file2)
src_files_to_mosaic

#合并两个文件
mosaic, out_trans = merge(src_files_to_mosaic)

#复制元数据
output_meta = file.meta.copy()

#更新元数据
output_meta.update(
    {"driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": out_trans,
    }
)

#写入目标
#将镶嵌栅格写入磁盘
out_fp = r"path\Nepal_Mosaic.tif"

with rasterio.open(out_fp, "w", **output_meta) as dest:
        dest.write(mosaic)

#打开合并的栅格数据
file_mosaic = rasterio.open(out_fp)

#读取数据
dataset_mosaic = file_mosaic.read()
print(file_mosaic.shape)

#显示数据
plt.imshow(dataset_mosaic[0], cmap='Spectral')
plt.show()

全球多分辨率地形高程数据

现在我们有了一个合并的全球多分辨率地形高程数据2010 GMTED2010，包括尼泊尔的所有数据，但文件还覆盖了尼泊尔之外的大部分周边地区。让我们通过使用尼泊尔的一个形状文件将区域限制在尼泊尔范围内。我们将使用一个包含世界各国边界的形状文件。您可以在此处下载此数据集。然后，我们使用mask函数剪裁栅格数据和形状文件。我们只使用形状文件的第一行和几何列。此操作的结果存储在clipped_array中，它是剪裁后的栅格数据，以及clipped_transform，它表示剪裁后的栅格的变换信息。

import geopandas as gpd
from shapely.geometry import mapping
from rasterio import mask as msk

#上传包含世界国家边界的形状文件
df = gpd.read_file(r'path/world-administrative-boundaries.shp')

#限制为尼泊尔
nepal = df.loc[df.name=="Nepal"]
nepal.head()

#剪裁数据
clipped_array, clipped_transform = msk.mask(file_mosaic, [mapping(nepal.iloc[0].geometry)], crop=True)

#

还有一个问题。栅格数据中的无数据值会严重破坏数据。因此，会扭曲我们的地图可视化，因为这些值是值范围的一部分。

理解问题

让我们按照以下方式解决这个问题，如本博客文章所述：

• 我们构建一个函数来处理无数据值。我们构造一个no-data参数来指定被认为是剪切数组中的“无数据”值。在这种情况下，它设置为（np.amax(clipped_array[0]) + 1），这意味着它等于剪切数组中的最大值加一。该值将被视为“无数据”值。
• 通过将剪切数组中的最小值的绝对值添加到剪切数组的第一个波段（索引0）来调整剪切数组。这一步确保剪切数组中的所有值都变成非负数。
• 我们还计算剪切数组的值范围。它将最大值和剪切数组中最小值的绝对值相加。value_range变量将保存计算得到的值范围。
• 使用基于现有颜色值字典（seismic字典）手动构建颜色值字典，并为“无数据”值定义背景颜色。
• 在最后一步中，我们使用新的颜色范围new_seismic绘制地图。

#让我们调查无数据值

nodata_value = file_mosaic.nodata 
print("无数据值:", nodata_value)
#无数据值: -32768.0

#将无数据值更改为最大高程的值加一
def clip_raster(gdf, img):
    clipped_array, clipped_transform = msk.mask(img, [mapping(gdf.iloc[0].geometry)], crop=True)
    clipped_array, clipped_transform = msk.mask(img, [mapping(gdf.iloc[0].geometry)],
                                                           crop=True, nodata=(np.amax(clipped_array[0]) + 1))
    clipped_array[0] = clipped_array[0] + abs(np.amin(clipped_array))
    value_range = np.amax(clipped_array) + abs(np.amin(clipped_array))
    return clipped_array, value_range

nepal_topography, value_range = clip_raster(nepal, file_mosaic)


#检查是否成功
print(value_range)


#让无数据值有新的背景颜色
from matplotlib import cm
from matplotlib.colors import ListedColormap,LinearSegmentedColormap

#Seismic
new_seismic = cm.get_cmap('seismic', 8828)

#定义背景颜色
background_color = np.array([0.9882352941176471, 0.9647058823529412, 0.9607843137254902, 1.0])

#使用颜色映射
newcolors = new_seismic(np.linspace(0, 1, 8828))

#将背景颜色添加到newcolors数组的最后一行。
newcolors = np.vstack((newcolors, background_color))

#使用新的Italy颜色映射
new_seismic = ListedColormap(newcolors)

#创建最终地图并保存
plt.figure(figsize=(10,10))
c = plt.imshow(nepal_topography[0], cmap = new_seismic)
clb = plt.colorbar(c, shrink=0.4)
clb.ax.set_title('高程（米）',fontsize=10)

plt.savefig(r'path\Topographic_Map_Nepal.png', bbox_inches='tight')
plt.show()

Voilá！我们有了一张能清楚显示尼泊尔不同海拔和三个地形区的地形图。

结论

您学会了使用来自美国地质调查局（USGS）的地理空间数据，在Python中生成地形图。您还了解到在最终数据集中处理缺失值对可视化的重要性。

决策者和从业人员现在可以将此地图用于进一步的分析，例如将其与贫困地图或自然灾害地图等其他地图结合起来，以分析是否存在某种联系。我们生成了一个有价值的工具，可以为政治中的基于证据的决策提供信息支持！

主要要点

• 地形图对基于证据的决策具有重要作用。
• 地形和海拔在城市规划、服务交付和不平等中起着至关重要的作用。
• Python具有分析地理空间数据的有用工具。
• 在此类数据中处理无数据值对于可视化至关重要。
• 可视化地理空间数据可以在细分级别上生成有价值的信息。

希望您觉得本文有帮助。欢迎在LinkedIn上与我联系。让我们共同努力，利用数据实现积极的变革。

常见问题

本文中显示的媒体不归Analytics Vidhya所有，仅由作者自行决定使用。