譯者 | 朱先忠
審校 | 重樓
我經常看到網上流傳著美麗的人口地圖；然而，我也常常會遇到一些技術問題，比如可視化本文中顯示的其他的地圖片段，或者將大規模光柵數據轉換為更便于計算的向量格式。在本文中，我將通過兩個主要全球人口數據來源的實踐來嘗試克服這其中的一些問題。
另一方面，同樣要注意，除了它們的美學價值外，顯示它們的人口數據和地圖是人們可以為任何城市發展或位置智能任務收集和整合的最基本和有價值的信息之一。它們在規劃新設施、選址和集水區分析、估計城市產品規模或分析不同社區等實踐應用中特別有用。
1.數據來源
在本文試驗中，我將依賴以下兩個細粒度的人口估計數據源，您可以通過所附鏈接來下載這些文件：
- 歐盟委員會的GHSL（全球人類住區層：https://ghsl.jrc.ec.europa.eu/ghs_pop2019.php）——用于測量每個網格單元的人口水平。從該數據源中可以找到整體描述以及我在他們2023年的報告中使用的空間分辨率為100米的特定集合。
- WorldPop中心。我將以德國為例，使用分辨率為100米的受限條件下的獨立國家數據集。你可以在鏈接https://hub.worldpop.org/geodata/listing?id=78處查找到完整的國家列表數據，還可以在鏈接https://hub.worldpop.org/geodata/summary?id=49789處查找到德國數據。
2.可視化全球人類住區層
2.1.導入數據
我第一次看到這個數據集是在“體系結構性能”部分的Datashader教程處。在復制了他們的可視化結果之后，在將其擴展到全球地圖時我遇到了一些麻煩，這些問題促使我開展了本文有關的研究工作。所以，接下來我將向您展示我是如何找到破解上述難題的解決方案的！
首先，我使用xarray包來解析光柵文件，代碼如下：
import rioxarray
file_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 ds
import xarray as xr
from colorcet import palette
from 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 = 1200
asp = 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雖然這段代碼在技術上看起來還可以，但我的2021款帶有16GB RAM的M1 Macbook Pro出現了一個糟糕的內存溢出錯誤。因此，讓我們裁剪一下圖像以便查看數據！為此，我關注上述數據源的“體系結構性能”部分，并專注于歐洲數據，這是我暫時關注的內容，因為這樣的選擇確實有效。
然而，我稍后要回答的主要問題是，盡管內存有限，但我們如何在使用本地機器的情況下可視化整個地球的數據呢？請先等一等！
import datashader as ds
import xarray as xr
from colorcet import palette
from datashader import transfer_functions as tf
import 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 = 1200
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")
img此代碼塊將輸出以下的視覺效果：
歐洲的人口分布圖（作者本人提供的圖片）
此繪制中，使用“火”色圖似乎是一個行業標準，這是有充分理由的；然而，如果你想把事情搞混，你可以在鏈接https://colorcet.holoviz.org/處找到其他配色方案，并使用類似于下面的編程方式：
#創建數據著色器畫布
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 = 20
downsampled_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 = 10
xstep = (maxx-minx) / N
ystep = (maxy-miny) / N
xsteps = list(np.arange(minx, maxx, xstep))
ysteps = list(np.arange(miny, maxy, ystep))現在，在每個x和y步驟上迭代，并創建每個圖像片段，其中每個圖像文件都以其在原始網格中的位置命名。此循環可能需要一段時間。
import os
foldout = '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('SAVED: ', x_coord, y_coord, y_coord+xstep,y_coord+ystep)最后，如果我們擁有所有的圖像片段，我們可以使用以下函數快速地把它們組合到一起。對于這段代碼，我還要求ChatGPT提供一些提示來加快速度，但和往常一樣，這個過程也需要一些手動調整。
from PIL import Image
def 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 = foldout
segment_width = size
segment_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年策劃的100米的分辨率，并向您展示如何從整個國家中劃出一個城市，并使用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 = 1200
asp = 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_union
import geopandas as gpd
admin = 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 pd
df_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)
# 從落入該邊界框的點創建點幾何圖形
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]})
# 構建一個GeoDataFrame
gdf_berlin = gpd.GeoDataFrame(geodata)
gdf_berlin = gpd.overlay(gdf_berlin, admin)然后，將人口可視化為向量數據：
import matplotlib.pyplot as plt
f, 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米分辨率的人口級別，分配給光柵文件中每個像素對應的每個點的幾何體。
總結
在這篇文章中，我探索了兩個全球人口數據集的可視化展示，它們通過結合各種近似、測量和建模方法，使用光柵網格以100米的顯著空間分辨率實現估計人口水平。這類信息對城市發展和位置智能的廣泛應用非常有價值，如基礎設施規劃、選址、社區概況等。從技術層面來看，我展示了三個空間層面的例子，涵蓋了整個地球，然后放大到國家，最后是城市。雖然該方法可以處理更小的分辨率，但這一切都發生在一臺筆記本電腦上已經令人非常滿意。另外注意到，編程實現過程中我使用了幾個強大的Python開源庫，如Xarray、DataShader和GeoPandas。
譯者介紹
朱先忠，51CTO社區編輯，51CTO專家博客、講師，濰坊一所高校計算機教師，自由編程界老兵一枚。
原文標題：Exploring Large-scale Raster Population Data，作者：Milan Janosov
