探索大规模栅格人口数据
大规模栅格人口数据探索
使用Python可视化多尺度地理人口数据：全球、国家和城市级数据
我经常在网上看到漂亮的人口地图；然而，我通常在一些技术部分卡住，比如在教程中展示的地图片段之外进行可视化，或者将大规模光栅数据转换为更易于计算的矢量格式。在本文中，我通过一个实践指南克服了其中一些缺点，介绍了两个主要的全球人口数据源。
值得注意的是，除了美学价值之外，人口数据和显示它们的地图是可以收集和整合到任何城市发展或位置智能任务中的最基本和有价值的信息之一。它们在规划新设施、选址和流域分析、估计城市产品规模或分析不同社区等应用场景中特别有用。
1. 数据来源
我依赖以下两个精细化人口估计数据源，您可以通过附带的链接下载文件（发布日期为本文发布日期）：
- 欧洲委员会的GHSL（全球人类住区层）衡量了每个网格单元的人口水平。在此处找到总体描述，并在此处找到我使用的来自2023年报告的空间分辨率为100米的特定数据集。
- WorldPop中心，我将使用受限的个别国家数据集以100米的分辨率以德国为例。在此处找到列表，并在此处找到德国的数据。
2. 可视化全球人类住区层
2.1. 导入数据！
我最初在“体系结构性能”的Datashader教程中遇到了这个数据集。在复现他们的可视化之后，我在将其扩展到全球地图时遇到了一些障碍，这促使我进行了这项工作，现在我将向您展示我发现的如何解决这些障碍的方法！
首先，使用xarray包解析栅格文件。
import rioxarrayfile_path = "GHS_POP_E2030_GLOBE_R2023A_54009_100_V1_0/GHS_POP_E2030_GLOBE_R2023A_54009_100_V1_0.tif" data_array = rioxarray.open_rasterio(file_path, chunks=True, lock=False)data_array该单元格的输出是数据集的详细描述：
2.2. 可视化数据的片段
我们已经可以看到对于大多数标准笔记本电脑来说，这是一个相当具有挑战性的数据量。不管怎样，让我们试着使用Datashader进行可视化，这是一个非常方便的处理这种规模的地理空间数据集的工具：
# 警告：此代码块很可能会导致内存溢出错误import datashader as dsimport xarray as xrfrom colorcet import palettefrom datashader import transfer_functions as tf# 准备绘图数据data_array_p = xr.DataArray(data_array)[0]data_array_p = data_array_p.where(data_array_p > 0)data_array_p = data_array_p.compute()# 获取图像大小size = 1200asp = data_array_p.shape[0] / data_array_p.shape[1]# 创建数据着色器画布cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_p)# 绘制图像cmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img尽管这段代码在技术上看起来没问题，但我搭载有16GB内存的2021款M1 MacBook Pro仍然会出现内存溢出错误。因此，让我们将图像裁剪为仅查看数据！为此，我遵循“体系结构性能”的做法，重点关注欧洲，暂时只关注欧洲，因为这肯定可行。
然而，我将在稍后回答的主要问题是，尽管存在内存限制，我们如何能够在本地机器上可视化整个地球的数据？请耐心等待！
import datashader as dsimport xarray as xrfrom colorcet import palettefrom datashader import transfer_functions as tfimport numpy as np# 裁剪数据数组data_array_c = data_array.rio.clip_box(minx=-1_000_000.0, miny=4_250_000.0, maxx=2_500_000.0, maxy=7_750_000.0)data_array_c = xr.DataArray(data_array_c)# 准备绘图data_array_c = xr.DataArray(data_array_c)[0]data_array_c = data_array_c.where(data_array_c > 0)data_array_c = data_array_c.compute()data_array_c = np.flip(data_array_c, 0)# 获取图像尺寸size = 1200asp = data_array_c.shape[0] / data_array_c.shape[1]# 创建data shader画布cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_c)# 绘制图像cmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img该代码块输出以下可视化结果：
使用“火焰”颜色映射看起来像是一个行业标准，有一个很好的理由；然而，如果您希望改变一下，您可以在此处找到其他颜色方案并应用如下：
# 创建data shader画布cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_c)# 绘制图像cmap = palette["bmw"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img该代码块输出以下可视化结果：
2.3. 可视化整个地球
因此，数据已经存在，但如果您手头只有一台普通的计算机，而且仍然希望以100米的分辨率可视化整个世界，那么我将在这里向您展示一个相对简单的解决方法——将整个栅格图像分割成大约一百个较小的瓦片，以便我的计算机逐一处理它们，然后使用一些图像处理技巧将它们合并成一个图像文件。
然而，在继续之前，简短说明一下。还有一种以以下方式降采样XArray数组的选项，但是我找不到一个适当的降采样方法来处理整个数据集。此外，我不想失去准确性，希望以完整的荣耀看到整个数据集。
# 快速降采样数据的方法downsampling_factor = 20downsampled_data_array = data_array.coarsen(x=downsampling_factor, y=downsampling_factor).mean()downsampled_data_array输出结果如下，值得与先前绘制的data_array进行对比：
要将整个栅格图像分割成网格段，首先获取其边界并定义N为步长。然后，创建图像段边界的列表。
minx = float(data_array.x.min().values)maxx = float(data_array.x.max().values)miny = float(data_array.y.min().values)maxy = float(data_array.y.max().values)N = 10xstep = (maxx-minx) / Nystep = (maxy-miny) / Nxsteps = list(np.arange(minx, maxx, xstep)) ysteps = list(np.arange(miny, maxy, ystep))现在，迭代每个x和y步骤，并创建每个图像段，其中每个图像文件的名称是根据其在原始网格中的位置命名。这个循环可能需要一些时间。
import osfoldout = 'world_map_image_segments'if not os.path.exists(foldout): os.makedirs(foldout)for idx_x, x_coord in enumerate(xsteps): for idx_y, y_coord in enumerate(ysteps): if not os.path.exists(foldout+'/'+str(idx_x)+'_'+str(idx_y)+'.png'): data_array_c = data_array.rio.clip_box( minx=x_coord, miny=y_coord, maxx=x_coord+xstep, maxy=y_coord+ystep) data_array_c = xr.DataArray(data_array_c)[0] data_array_c = data_array_c.fillna(0) data_array_c = data_array_c.where(data_array_c > 0) data_array_c = data_array_c.compute() data_array_c = np.flip(data_array_c, 0) size = 2000 asp = data_array_c.shape[0] / data_array_c.shape[1] cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size)) raster = cvs.raster(data_array_c) cmap = palette["fire"] img = tf.shade(raster, how="eq_hist", cmap=cmap) img = tf.set_background(img, "black") pil_image = img.to_pil() pil_image.save(foldout+'/'+str(idx_x)+'_'+str(idx_y)+ '.png') print('保存成功: ', x_coord, y_coord, y_coord+xstep,y_coord+ystep)最后，如果我们有了所有的图像段，我们可以使用以下函数快速组装它们。对于这段代码，我也向ChatGPT寻求了一些提示，以加快速度，但像往常一样，还需要一些手动调整。
from PIL import Imagedef find_dimensions(image_dir): max_x = 0 max_y = 0 for filename in os.listdir(image_dir): if filename.endswith(".png"): x, y = map(int, os.path.splitext(filename)[0].split("_")) max_x = max(max_x, x) max_y = max(max_y, y) return max_x + 1, max_y + 1 image_dir = foldoutsegment_width = sizesegment_height = int(asp*size)# 确定大图的尺寸large_image_width, large_image_height = find_dimensions(image_dir)# 创建一个空的大图（白色背景）large_image = Image.new("RGB", (large_image_width * segment_width, large_image_height * segment_height), "black")# 循环遍历各个图像段，并将它们粘贴到大图中for filename in sorted(os.listdir(image_dir)): if filename.endswith(".png"): x, y = map(int, os.path.splitext(filename)[0].split("_")) segment_image = Image.open(os.path.join(image_dir, filename)) x_offset = x * segment_width y_offset = large_image_height * segment_height-1*y * segment_height large_image.paste(segment_image, (x_offset, y_offset))# 保存合并后的大图large_image.save("global_population_map.png") 这里是最终结果，整个地球的映射：
3. 可视化和转换WorldPop数据
我想展示给你的第二个数据源是WorldPop人口数据库，它以大洲和国家在不同分辨率下分开。在这个例子中，为了补充前一节的大洲和全球级别，我选择了德国和2020年的100m分辨率，并向你展示如何从整个国家中切出一个城市，并将其转换为易于使用的矢量格式使用GeoPandas。
3.1. 可视化WorldPop数据
使用之前的方法，我们可以快速可视化这个栅格文件：
# 解析数据data_file = 'deu_ppp_2020_constrained.tif'data_array = rioxarray.open_rasterio(data_file, chunks=True, lock=False)# 准备数据data_array = xr.DataArray(data_array)[0]data_array = data_array.where(data_array > 0)data_array = data_array.compute()data_array = np.flip(data_array, 0)# 获取图片尺寸size = 1200asp = data_array.shape[0] / data_array.shape[1]# 创建数据着色画布cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array)# 绘制图片cmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img此代码块输出以下可视化结果：
3.2. 转换WorldPop数据
在可视化整个地球、欧洲大陆和德国国家之后，我想更深入地研究柏林这个城市，并向您展示如何将这种栅格数据转换为矢量格式，并使用GeoPandas轻松操作它。为此，我在此处以geojson格式访问了柏林的行政边界。
此行政文件包含柏林的区。首先，我将它们合并为整个城市。
from shapely.ops import cascaded_unionimport geopandas as gpdadmin = gpd.read_file('tufts-berlin-bezirke-boroughs01-geojson.json')admin = gpd.GeoDataFrame(cascaded_union(admin.geometry.to_list()), columns = ['geometry']).head(1)admin.plot()此代码块输出以下可视化结果：
现在将xarray转换为Pandas DataFrame，提取几何信息，并构建一个GeoPandas GeoDataFrame。一种方法是：
import pandas as pddf_berlin = pd.DataFrame(data_array.to_series(), columns = ['population']).dropna()现在，从中构建一个GeoDataFrame，重点关注柏林：
from shapely.geometry import Point# 为了更容易地选择坐标，找到限定的边界框minx, miny, maxx, maxy = admin.bounds.T[0].to_list()points = []population = df_berlin.population.to_list()indicies = list(df_berlin.index)# 从落入此边界框中的点创建Point几何geodata = []for ijk, (lon, lat) in enumerate(indicies): if minx <= lat <= maxx and miny <= lon <= maxy: geodata.append({'geometry' : Point(lat, lon), 'population' : population[ijk]}) # 构建一个GeoDataFramegdf_berlin = gpd.GeoDataFrame(geodata)gdf_berlin = gpd.overlay(gdf_berlin, admin)然后，将人口可视化为矢量数据：
import matplotlib.pyplot as pltf, ax = plt.subplots(1,1,figsize=(15,15))admin.plot(ax=ax, color = 'k', edgecolor = 'orange', linewidth = 3)gdf_berlin.plot(column = 'population', cmap = 'inferno', ax=ax, alpha = 0.9, markersize = 0.25)ax.axis('off')f.patch.set_facecolor('black')此代码块输出以下可视化结果：
最后，我们有一个标准的GeoDataFrame，每个像素对应的点几何体上分配了100米分辨率的人口水平。
总结
在本文中，我探索了两个基于各种估计方法的全球人口数据集，这些方法结合了各种逼近、测量和建模方法，以出色的空间分辨率（100m）使用栅格网格来估计人口水平。这种类型的信息在城市发展和位置智能的广泛应用中非常有价值，例如基础设施规划、选址、社区剖析等。从技术角度来看，我展示了三个空间层次的示例，包括整个地球，然后缩小到国家，最后是城市。虽然这种方法甚至可以处理更小的分辨率，但这一切都在一台笔记本电脑上使用强大的Python库（如Xarray、DataShader和GeoPandas）完成。
