譯者 | 朱先忠
審校 | 重樓
我經(jīng)?？吹骄W(wǎng)上流傳著美麗的人口地圖；然而，我也常常會(huì)遇到一些技術(shù)問(wèn)題，比如可視化本文中顯示的其他的地圖片段，或者將大規(guī)模光柵數(shù)據(jù)轉(zhuǎn)換為更便于計(jì)算的向量格式。在本文中，我將通過(guò)兩個(gè)主要全球人口數(shù)據(jù)來(lái)源的實(shí)踐來(lái)嘗試克服這其中的一些問(wèn)題。
另一方面，同樣要注意，除了它們的美學(xué)價(jià)值外，顯示它們的人口數(shù)據(jù)和地圖是人們可以為任何城市發(fā)展或位置智能任務(wù)收集和整合的最基本和有價(jià)值的信息之一。它們?cè)谝?guī)劃新設(shè)施、選址和集水區(qū)分析、估計(jì)城市產(chǎn)品規(guī)?；蚍治霾煌鐓^(qū)等實(shí)踐應(yīng)用中特別有用。
1.數(shù)據(jù)來(lái)源
在本文試驗(yàn)中，我將依賴以下兩個(gè)細(xì)粒度的人口估計(jì)數(shù)據(jù)源，您可以通過(guò)所附鏈接來(lái)下載這些文件：
- 歐盟委員會(huì)的GHSL（全球人類住區(qū)層：https://ghsl.jrc.ec.europa.eu/ghs_pop2019.php）——用于測(cè)量每個(gè)網(wǎng)格單元的人口水平。從該數(shù)據(jù)源中可以找到整體描述以及我在他們2023年的報(bào)告中使用的空間分辨率為100米的特定集合。
- WorldPop中心。我將以德國(guó)為例，使用分辨率為100米的受限條件下的獨(dú)立國(guó)家數(shù)據(jù)集。你可以在鏈接https://hub.worldpop.org/geodata/listing?id=78處查找到完整的國(guó)家列表數(shù)據(jù)，還可以在鏈接https://hub.worldpop.org/geodata/summary?id=49789處查找到德國(guó)數(shù)據(jù)。
2.可視化全球人類住區(qū)層
2.1.導(dǎo)入數(shù)據(jù)
我第一次看到這個(gè)數(shù)據(jù)集是在“體系結(jié)構(gòu)性能”部分的Datashader教程處。在復(fù)制了他們的可視化結(jié)果之后，在將其擴(kuò)展到全球地圖時(shí)我遇到了一些麻煩，這些問(wèn)題促使我開(kāi)展了本文有關(guān)的研究工作。所以，接下來(lái)我將向您展示我是如何找到破解上述難題的解決方案的！
首先，我使用xarray包來(lái)解析光柵文件，代碼如下：
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此代碼片斷的輸出結(jié)果是對(duì)數(shù)據(jù)集的一段詳細(xì)描述：
2.2.可視化數(shù)據(jù)段
我們已經(jīng)看到，對(duì)于大多數(shù)標(biāo)準(zhǔn)筆記本電腦來(lái)說(shuō)，這是一個(gè)非常具有挑戰(zhàn)性的數(shù)據(jù)量。無(wú)論如何，讓我們?cè)囍褂肈atashader來(lái)可視化數(shù)據(jù)，這是一個(gè)非常方便的工具，適用于這種規(guī)模的地理空間數(shù)據(jù)集的展示：
#警告：此處代碼塊很可能會(huì)導(dǎo)致您的計(jì)算機(jī)內(nèi)存溢出錯(cuò)誤
import datashader as ds
import xarray as xr
from colorcet import palette
from datashader import transfer_functions as tf
# 準(zhǔn)備繪圖
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]
#創(chuàng)建數(shù)據(jù)著色器畫(huà)布
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雖然這段代碼在技術(shù)上看起來(lái)還可以，但我的2021款帶有16GB RAM的M1 Macbook Pro出現(xiàn)了一個(gè)糟糕的內(nèi)存溢出錯(cuò)誤。因此，讓我們裁剪一下圖像以便查看數(shù)據(jù)！為此，我關(guān)注上述數(shù)據(jù)源的“體系結(jié)構(gòu)性能”部分，并專注于歐洲數(shù)據(jù)，這是我暫時(shí)關(guān)注的內(nèi)容，因?yàn)檫@樣的選擇確實(shí)有效。
然而，我稍后要回答的主要問(wèn)題是，盡管內(nèi)存有限，但我們?nèi)绾卧谑褂帽镜貦C(jī)器的情況下可視化整個(gè)地球的數(shù)據(jù)呢？請(qǐng)先等一等！
import datashader as ds
import xarray as xr
from colorcet import palette
from datashader import transfer_functions as tf
import numpy as np
# 裁剪數(shù)據(jù)陣列
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)
# 準(zhǔn)備繪圖
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]
# 創(chuàng)建數(shù)據(jù)著色器畫(huà)布
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此代碼塊將輸出以下的視覺(jué)效果：
歐洲的人口分布圖（作者本人提供的圖片）
此繪制中，使用“火”色圖似乎是一個(gè)行業(yè)標(biāo)準(zhǔn)，這是有充分理由的；然而，如果你想把事情搞混，你可以在鏈接https://colorcet.holoviz.org/處找到其他配色方案，并使用類似于下面的編程方式：
#創(chuàng)建數(shù)據(jù)著色器畫(huà)布
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此代碼塊輸出以下形式的視覺(jué)效果：
歐洲的人口分布圖另一種色圖（作者本人提供的圖片）
2.3.可視化全球人口數(shù)據(jù)
到此，我們已經(jīng)取得了全球人口數(shù)據(jù)，但如果你手邊僅有一臺(tái)普通的電腦，仍然想以100米的分辨率可視化整個(gè)世界，那該怎么辦呢？我將在這里向您展示的解決方法非常簡(jiǎn)單——我將整個(gè)光柵圖像分割成大約一百個(gè)較小的片斷。這樣一來(lái)，我的計(jì)算機(jī)就可以一個(gè)接一個(gè)地很好地處理它們，然后使用一些圖像處理技巧將它們合并到一個(gè)圖像文件中。
然而，在繼續(xù)介紹之前，還有一個(gè)細(xì)節(jié)需要說(shuō)明。我們可以通過(guò)以下方式對(duì)XArray數(shù)組降低采樣率——然而，我找不到一個(gè)合適的降低尺度的方法來(lái)處理整個(gè)數(shù)據(jù)集。此外，我不想失去準(zhǔn)確性，希望看到整個(gè)數(shù)據(jù)集的真實(shí)面貌。
# 對(duì)數(shù)據(jù)進(jìn)行降低采樣率的快速方法
downsampling_factor = 20
downsampled_data_array = data_array.coarsen(x=downsampling_factor, y=downsampling_factor).mean()
downsampled_data_array最后的輸出結(jié)果不錯(cuò)，不亞于之前繪制的data_array：
為了實(shí)現(xiàn)將整個(gè)光柵圖像分割為網(wǎng)格段，首先，獲取其邊界并將N定義為步長(zhǎng)。然后，創(chuàng)建圖像片段邊界列表。
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))現(xiàn)在，在每個(gè)x和y步驟上迭代，并創(chuàng)建每個(gè)圖像片段，其中每個(gè)圖像文件都以其在原始網(wǎng)格中的位置命名。此循環(huán)可能需要一段時(shí)間。
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)最后，如果我們擁有所有的圖像片段，我們可以使用以下函數(shù)快速地把它們組合到一起。對(duì)于這段代碼，我還要求ChatGPT提供一些提示來(lái)加快速度，但和往常一樣，這個(gè)過(guò)程也需要一些手動(dòng)調(diào)整。
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)
# 創(chuàng)建一個(gè)空的大圖像（白色背景）
large_image = Image.new("RGB", (large_image_width * segment_width, large_image_height * segment_height), "black")
# 循環(huán)瀏覽各個(gè)圖像片段并將它們粘貼到大圖像中
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") 最后的結(jié)果是，整個(gè)全球數(shù)據(jù)都被成功繪制出來(lái)了：
全球人口分布（作者本人的圖片）
3.可視化和轉(zhuǎn)換WorldPop數(shù)據(jù)
我想向大家展示的第二個(gè)數(shù)據(jù)來(lái)源是WorldPop人口數(shù)據(jù)庫(kù)，它以不同的分辨率分別列出了各大洲和相應(yīng)的國(guó)家。在這個(gè)例子中，為了補(bǔ)充前面繪制的所有大陸和全球級(jí)別的數(shù)據(jù)情況，我在這一部分中主要集中在各國(guó)家和相應(yīng)城市級(jí)別數(shù)據(jù)上。例如，我選擇了德國(guó)和2020年策劃的100米的分辨率，并向您展示如何從整個(gè)國(guó)家中劃出一個(gè)城市，并使用GeoPandas庫(kù)將其轉(zhuǎn)化為易于使用的向量格式。
3.1.可視化WorldPop數(shù)據(jù)
使用與前面相同的方法，我們可以再次快速可視化這一部分的光柵文件：
#分析數(shù)據(jù)
data_file = 'deu_ppp_2020_constrained.tif'
data_array = rioxarray.open_rasterio(data_file, chunks=True, lock=False)
# 準(zhǔn)備數(shù)據(jù)
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]
#創(chuàng)建數(shù)據(jù)著色器畫(huà)布
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此代碼片斷將輸出以下視覺(jué)效果：
德國(guó)的人口分布圖（作者本人的圖片）
3.2.轉(zhuǎn)換WorldPop數(shù)據(jù)
在可視化了整個(gè)地球、歐洲大陸和德國(guó)之后，我想更多地了解一下柏林市信息，并向您展示如何將這些光柵數(shù)據(jù)轉(zhuǎn)換為向量格式，并使用GeoPandas輕松操作。為此，我在這里以geojson格式訪問(wèn)柏林的行政邊界。
這個(gè)管理文件包含柏林的行政區(qū)，所以首先，我將它們合并為一個(gè)整體。
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()此代碼塊輸出以下所示的視覺(jué)效果：
柏林的行政邊界圖（作者本人的圖片）
現(xiàn)在，將xarray轉(zhuǎn)換為Pandas DataFrame，提取幾何體信息，并構(gòu)建GeoPandas GeoDataFrame。一種方法是：
import pandas as pd
df_berlin = pd.DataFrame(data_array.to_series(), columns = ['population']).dropna()現(xiàn)在，在此基礎(chǔ)上構(gòu)建一個(gè)GeoDataFrame，重點(diǎn)關(guān)注柏林信息：
from shapely.geometry import Point
#找到限制邊界框以便于坐標(biāo)選擇
minx, miny, maxx, maxy = admin.bounds.T[0].to_list()
points = []
population = df_berlin.population.to_list()
indicies = list(df_berlin.index)
# 從落入該邊界框的點(diǎn)創(chuàng)建點(diǎn)幾何圖形
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]})
# 構(gòu)建一個(gè)GeoDataFrame
gdf_berlin = gpd.GeoDataFrame(geodata)
gdf_berlin = gpd.overlay(gdf_berlin, admin)然后，將人口可視化為向量數(shù)據(jù)：
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')此部分代碼塊輸出以下視覺(jué)效果：
柏林的人口分布圖（作者本人的圖片）
最后，我們得到了一個(gè)標(biāo)準(zhǔn)的GeoDataFrame，它具有100米分辨率的人口級(jí)別，分配給光柵文件中每個(gè)像素對(duì)應(yīng)的每個(gè)點(diǎn)的幾何體。
總結(jié)
在這篇文章中，我探索了兩個(gè)全球人口數(shù)據(jù)集的可視化展示，它們通過(guò)結(jié)合各種近似、測(cè)量和建模方法，使用光柵網(wǎng)格以100米的顯著空間分辨率實(shí)現(xiàn)估計(jì)人口水平。這類信息對(duì)城市發(fā)展和位置智能的廣泛應(yīng)用非常有價(jià)值，如基礎(chǔ)設(shè)施規(guī)劃、選址、社區(qū)概況等。從技術(shù)層面來(lái)看，我展示了三個(gè)空間層面的例子，涵蓋了整個(gè)地球，然后放大到國(guó)家，最后是城市。雖然該方法可以處理更小的分辨率，但這一切都發(fā)生在一臺(tái)筆記本電腦上已經(jīng)令人非常滿意。另外注意到，編程實(shí)現(xiàn)過(guò)程中我使用了幾個(gè)強(qiáng)大的Python開(kāi)源庫(kù)，如Xarray、DataShader和GeoPandas。
譯者介紹
朱先忠，51CTO社區(qū)編輯，51CTO專家博客、講師，濰坊一所高校計(jì)算機(jī)教師，自由編程界老兵一枚。
原文標(biāo)題：Exploring Large-scale Raster Population Data，作者：Milan Janosov
