全球人口与洋流数据可视化

全球人口与洋流数据可视化

一、题目要求

示例图

二、准备工作

1. 数据格式分析

   本次实验的数据分为两批，第一批数据为GIS全球人口数据，类型为tif文件，解压后大小为2GB 。其中，TIFF文件以**.tif**为扩展名。其数据格式是一种3级体系结构，Ti内部结构可以分成三个部分，分别是：文件头信息区、标识信息区和图像数据区。其中所有的标签都是以升序排列，这些标签信息是用来处理文件中的图像信息的。

   第二批数据文件为按年平均的海洋表面洋流数据，类型为mat文件，可以直接使用python的sio.loadmat 方法导入。

2. 实验方法

   编程语言：Python3.9

   编译环境：Anaconda3

   编译器：Spyder

3. 基础知识

   大数据是指无法在一定时间内用常规软件工具对其内容进行抓取、管理和处理的数据集合。本次实验中所采用的tif文件，虽然没有达到目前工业和实验室级别的存储单位，但与之前的实验相比，已经显得十分“庞大”。对于大数据，由于设备的内存限制无法直接读取，因此可以采用分块读入的方式。

   洋流，是指海水沿着一定方向有规律的具有相对稳定速度的水平流动，是从一个海区水平或垂直地向另一个海区大规模的非周期性的运动，是海水的主要运动形式。

   在可视化中，洋流的绘制实际上是一个矢量场的绘制与优化问题。

   高中地理知识讲解：洋流 - 知乎 (zhihu.com)

   洋流地图：可视化我们的海洋运动 - 知乎 (zhihu.com)

   

   洋流对照图

三**、关键步骤**

1、大数据：全球人口

1.1 数据预处理

• 数据读取

考虑到本次tif数据过大，无法直接读取，因此先通过meta解析出元数据的信息头，获取到一些关键信息，以此为依据再做进一步的处理。

file_path = "GHS_POP_GPW42000_GLOBE_R2015A_54009_250_v1_0.tif"
with rasterio.open(file_path) as src:
    meta = src.meta
		print(meta)

获取的文件头信息整理如下：

💡 **TIF文件头Info：**

• driver: 文件的驱动程序，表示文件格式。在这种情况下，文件格式是GeoTIFF（GTiff）。
• dtype: 文件中像素数据的数据类型。在这种情况下，数据类型是32位浮点数（float32）。
• nodata: 指定像素数据中的无数据值。在这种情况下，该文件未指定无数据值（None）。
• width: 文件中像素的宽度，以像素为单位。在这种情况下，宽度为141988像素。
• height: 文件中像素的高度，以像素为单位。在这种情况下，高度为60942像素。
• count: 文件中的波段数。在这种情况下，文件只有一个波段。
• crs: 文件中像素的坐标参考系统。在这种情况下，坐标参考系统是一个世界莫尔韦德（World Mollweide）投影，采用WGS 84椭球体，并使用米作为单位。
• transform: 文件中像素坐标和实际坐标之间的转换，以Affine对象的形式表示。在这种情况下，每个像素的大小为250米，而左下角的像素的实际坐标为(-17619594.54744353, 8750529.46186849)。

• 数据分块读入

width和height指明了tif文件像素的尺寸，分别是141988和60942，crs是文件中像素的坐标参考系统，本文件是一个世界莫尔韦德投影。考虑到这是一个以像素为基本单位构成的图像文件，对于大尺寸，可以规定区块大小（block_size)，用尺寸除以区块大小，来获得区间个数，并分块读入。

我们先指定size（height,width)为尺寸大小，则平面直角坐标系下的横向和纵向分块数为：

$$\begin{gathered} nu{m}_{-}bloc{k}_{-}x=ceil(width/bloc{k}_{-}size) \\ nu{m}_{-}bloc{k}_{-}y=ceil(height/bloc{k}_{-}size) \end{gathered}$$

对于每一个分块，x_min为区间左端点， x_max为区间右端点，同理，y_min为区间上端点，y_max为区间下端点，对于第（x,y)个区间：

$$\begin{gathered} x_{-}min=x\ast bloc{k}_{-}size \\ {y}_{-}min=y\ast bloc{k}_{-}size \\ {x}_{-}max=min((x+1)\ast bloc{k}_{-}size,width)) \\ {y}_{-}max=min((y+1)\ast bloc{k}_{-}size,height)) \end{gathered}$$

由于文件的尺寸大小并不一定能完全分成整数个区间，因此在x_max、y_max的计算中，我们还需要比较计算值与尺寸的大小关系，取更小的那个值，防止数组访问溢出。

在src源文件中以指定窗口大小的方式读入数据到data。

data = src.read(1, window=rasterio.windows.Window(x_min, y_min, x_max - x_min, y_max - y_min))

• 数据整合

分块读入数据后，相当于将原本的数据压缩了block_size 倍，以此达到缩小数据量的效果。本文采用对读入区域取平均值的方法来进行数据压缩，使(height,width)尺寸的数据缩放到（num_blocks_y,num_blocks_x) 。

伪代码：

输入:待处理的tif数据src；tif数据的元数据meta；分块大小block_size
输出:tif数据分块统计平均值的数组result_array

1.计算x和y方向上分块数
2.遍历每个分块,计算其对应的x和y方向的坐标范围
3.读取对应范围内的栅格数据
4.计算栅格数据的平均值并存入结果数组

1.2 数据可视化

在获得压缩后的数组result_array 后，便可以直接进行可视化了。通过python提供的imshow（）方法可以直接指定cmap颜色映射来绘制热度图。

plt.imshow(np.log(result_array), cmap='Blues')
plt.colorbar()
plt.show()

由于数据本身数值比较接近，可视化结果对比并不明显，因此取对数操作，加强对比度。

1.3 区域高精度展示

• 定义RectangleSelector() 来选择区域

上文已经实现了按照分块来读取tif图绘制热度图，对于更高精度的区域展示，可以利用矩形区域选择工具来划定范围，并到源数据中读取划定部分，以更小的block_size读取，这样就能够实现。

划定区域用的是Python提供的RectangleSelector() 方法，该方法指定了绘制的形状，调用一个响应函数Onselect() ;

rect_selector = RectangleSelector(ax, onselect, drawtype='box', useblit=True,
button=[1], minspanx=5, minspany=5, spancoords='pixels',interactive=True)

• 根据选择区域再次分块读取数据

假设我们划定了一个矩形区域，左上端点在画布中的坐标为（x1，y1），右下端点在画布中的坐标为(x2,y2)。由于分块读入，实际上画布的坐标相比原本在数据中的位置已经缩小了block_size 倍。因此(x1,y1)和(x2,y2)在源数据中的位置为：

$$\begin{gathered} x_1,y_1=x_1\ast bloc{k}_{-}size,y_1\ast bloc{k}_{-}size \\ {x}_2,y_2=x_2\ast bloc{k}_{-}size,y_2\ast bloc{k}_{-}size \end{gathered}$$

然后，定义一个比原本区块大小要更小的block_size1 ，比如说16；再运用类似上文的逻辑来分块读入所划定的区域，在新的画布中绘制出来。

• 成果图

1.4 实时查看和局部缩放

• 工具选择

在上面的实验步骤中，本文已经实现了分块读入、选择区域进行高精度展示。但每次选择后，都要等待大量运算，才能显示出划定区域。如果想要做到实时查看的效果， 只运用Python来对数据进行读取就有些低效。因此，本文考虑了多种工具，如Arcgis、Keplergl、leaflet等，最终选择Leaflet作为绘制工具。

• 思路

对于大数据，唯一的制约即是运算速度，如果再运用分块读入的方式进行展示，那必然不能达到“实时”的效果。因此，对数据进行预处理，输出成其他文件格式再进行读取操作，会更为实际一些。

金字塔 TIFF（Pyramid TIFF）是一种特殊的 TIFF 图像格式，用于高效地处理和显示大型图像。金字塔 TIFF 文件包含多个图像级别（或分辨率），每个级别是前一个级别的缩小版本。这些级别按照从原始图像（最底层）到最小尺寸（最顶层）的顺序，形成一个金字塔状结构。

图像金字塔

图像金字塔是多尺度表达图像的一种措施，它巧妙地避免了实时读取数据，而是指定不同精度等级的数据，只展现出这几个等级，而忽略中间的变化过程，通过级别变换来实现缩放。

瓦片数据（Tile）是一种将地图数据切分为一系列小方块的技术，使得地图可以快速加载和渲染。每个小方块被称为“瓦片”（Tile），通常是256像素x256像素大小的图片，以PNG或JPEG格式存储。这些瓦片按照网格排列，组成了一张完整的地图。当地图进行缩放操作时，地图会从瓦片服务器上请求不同层级、不同位置的瓦片，以保证地图的连续性和流畅性。瓦片数据可以极大地提高地图的加载速度和用户体验，也是很多在线地图服务的核心技术之一。

• 尝试

一开始我试着用Arcgis处理软件对tif文件进行读取，转换成金字塔，并打包成tpk文件，以此直接生成地图。但是似乎在服务器配置的过程中出了一些问题，因此转而使用Leaflet。

• 数据预处理

Leaflet是一个流行的开源JavaScript库，用于在交互式地图上展示地理空间数据。由于在探索的过程中，瓦片数据给予了我很大的启发。虽然我不能生成完美的瓦片数据，但是我可以继续利用我之前根据不同分块数绘制不同精度人口数据图的方法，生成几张分辨率逐级升高的图片，然后利用Leaflet库进行操作。

首先利用Python读取数据，生成block_size = [256,128,64,32,16]的五张分级图片。之所以不再升高分辨率，是因为block_size=8时，内存已经无法承受了。

file_path = "GHS_POP_GPW42000_GLOBE_R2015A_54009_250_v1_0.tif"
output_folder = "tiles1"

if not os.path.exists(output_folder):
    os.makedirs(output_folder)

with rasterio.open(file_path) as src:
    meta = src.meta
    print(meta)

    # 定义不同等级的 block size
    block_sizes = [256,128,64,32,16]

    for level, block_size in enumerate(block_sizes):
        num_blocks_x = int(np.ceil(meta['width'] / block_size))
        num_blocks_y = int(np.ceil(meta['height'] / block_size))
        results = []

        for y in range(num_blocks_y):
            for x in range(num_blocks_x):
                x_min = x * block_size
                y_min = y * block_size
                x_max = min((x + 1) * block_size, meta['width'])
                y_max = min((y + 1) * block_size, meta['height'])
                data = src.read(1, window=rasterio.windows.Window(x_min, y_min, x_max - x_min, y_max - y_min))
                result = np.mean(data)
                results.append(result)
        result_array = np.array(results).reshape(num_blocks_y, num_blocks_x)
        log_result_array = np.log1p(result_array)
        
        # 保存图像
        output_path = os.path.join(output_folder, f"level_5.png")
        plt.imsave(output_path, log_result_array, cmap='Blues')

        print(f"Saved level {level} image: {output_path}")

在输出文件夹下方，生成了5张图片：

可以看到，分辨率逐级上升。然后，打开Idea，新建一个HTML文件，在HTML中输入下列代码：

<html>
<head>
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.1/dist/leaflet.css">
    <style>
        #image-map {
            width: 100%;
            height: 100%;
            border: 1px solid #ccc;
            margin-bottom: 10px;
        }
    </style>
</head>
<body>
<div id="image-map"></div>

<script src="https://unpkg.com/leaflet@1.3.1/dist/leaflet.js"></script>
<script>
    // Using leaflet.js to pan and zoom a big image.
    // create the slippy map
    var map = L.map('image-map', {
        minZoom: 1,
        maxZoom: 5,
        center: [0, 0],
        zoom: 1,
        crs: L.CRS.Simple
    });

    // dimensions of the image
    var w = 6000,
        h = 3500;

    // calculate the edges of the image, in coordinate space
    var southWest = map.unproject([0, h], map.getMaxZoom()-1);
    var northEast = map.unproject([w, 0], map.getMaxZoom()-1);
    var bounds = new L.LatLngBounds(southWest, northEast);

    // Define an array of image URLs for different zoom levels
    var imageURLs = [
        'tiles/level_0.png',
        'tiles/level_1.png',
        'tiles/level_2.png',
        'tiles/level_3.png',
        'tiles/level_4.png'
    ];

    // Set initial image overlay
    var imageOverlay = L.imageOverlay('tiles/level_0.png', bounds).addTo(map);
    // L.imageOverlay('tiles/level_0.png', bounds).addTo(map);

    // tell leaflet that the map is exactly as big as the image
    map.setMaxBounds(bounds);

    // Function to update image overlay based on current zoom level
    function updateImageOverlay() {
        var zoom = map.getZoom();
        var newImageOverlay = L.imageOverlay(imageURLs[zoom - 1], bounds);
        map.removeLayer(imageOverlay);
        imageOverlay = newImageOverlay.addTo(map);
    }

    // Attach zoom event to the map
    map.on('zoomend', function() {
        updateImageOverlay();
    });
</script>
</body>
</html>

在这段代码中，我为不同的缩放级别定义了一个名为 imageURLs 的数组。数组中的每个元素都是一个 URL，指向对应缩放级别的图像文件。使用 map.on('zoomend', updateImageOverlay) 为地图添加了一个事件侦听器，当缩放操作结束时，将调用 updateImageOverlay 函数。updateImageOverlay 函数会根据当前的缩放级别更新图像覆盖层。这样，就可以实现伪地图的效果。

1.5 效果图展示

等级1

等级2

等级3

等级4

等级5

2、向量场可视化

1.1 数据预处理

• 数据读取

首先，使用SciPy库中的**sio.loadmat函数读取drifter_annualmeans.mat**文件中的数据，提取经度、纬度、经向流速和纬向流速信息。

data = sio.loadmat('drifter_annualmeans.mat')
lon = data['Lon']
lat = data['Lat']
u = data['U']
v = data['V']

洋流.mat

• 数据预处理

仔细观察读取到的U和V数据，发现有许多空值。这来源于浮标数据采集时未采集到的空缺数据，可以不做处理，在绘制图像时再加以判断即可。

• 添加地图特征

接下来，使用Matplotlib和Cartopy库创建一个地图投影，并添加地理特征，如陆地、海岸线和国界。

fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='none', alpha=0.3)
ax.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=0.5)
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=0.5)

1.2 streamplot () 方法直接绘制流线图

python的matplotlib 库自带了许多绘制矢量场的工具，如streamplot、quiver和contour等。其中streamplot绘制流线，quiver绘制方向场而contour则绘制等高线。

ax.streamplot(lon, lat, u, v, transform=ccrs.PlateCarree(), density=4)

其中，transform指定地理坐标的映射方式，density负责控制绘制流线的密度。

Streamplot with density = 8

1.3 复现streamplot()方法

• 参数定义

为了绘制洋流图像，本文定义了一个名为simple_streamplot(ax, lon, lat, u, v, density=20, color='k', linewidth=1, start_points = None)的函数。该函数接收地图轴、lat、lon、u、v等参数作为数据；density控制流线的绘制密度，start_points控制起始点。

• 流线生长的逻辑

选择起始点是绘制洋流的第一步，由于提供的数据点包含了海洋和陆地的经纬度坐标，因此，我们选择数据的边界点作为起始点，随机选择1000个点，添加到border_points数组中。

生成起始点后，对每一个起始点(x0,y0)进行遍历，并在读取的数据中寻找最接近当前起始点的经纬度点，读取该点对应的u、v，记为du、dv。

$$\begin{gathered} i,j=min(dist(lon,x0)).index,min(dist(lat,y0)).index \\ du=u[j,i] \\ dv=v[j,i] \end{gathered}$$

令当前的(x0,y0) =(x_current,y_current) ,新建line列表，该列表是个不断向其中填充点的点集。

设置迭代次数step，在step次循环中，都令：

$$\begin{gathered} x_{-}current+=du\ast 0.1 \\ {y}_{-}current+=dv\ast 0.1 \end{gathered}$$

以此来不断延伸流线的范围，直到控制当前点的坐标越界为止，退出；若未越界，则迭代完成退出。这样，一条以(x0,y0) 为起始点的流线就生长完成。

• 绘制流线

对于每一个起始点的循环，最终都会生成一个line列表，其中添加了所有绘制本条流线的控制点。用zip(*line)将line转换成绘制折线所需要的数据形式，最后用plot绘制。由于每一步更新生长点的步长很小（0.1），因此绘制的折线可以近似看成曲线。

ax.plot(*zip(*line), color=color, linewidth=linewidth)

• 伪代码

1.如果start_points为None，则根据lon和lat的边界和密度生成起点数组border_points和random_points，并拼接成起点数组start_points。
2.将lon和lat转为一维数组。
3.遍历起点数组中的每个起点：
	1.在经度数组中查找当前起点的索引i，在纬度数组中查找当前起点的索引j。
	2.如果i或j超出数组范围或者在边界上，则跳过当前起点。
	3.根据当前起点在u和v数组中查找相应的du和dv。
	4.根据du和dv计算下一个点的坐标x_current和y_current，并将该点添加到流线line中。
	5.在经度数组中查找当前点的索引i，在纬度数组中查找当前点的索引j。
	6.如果i或j超出数组范围，则跳出循环。
	7.根据当前点在u和v数组中查找相应的du和dv。
4.在绘图对象ax中绘制所有流线。

• 成果图

效果仍然有待改进。

1.4 区域高精度向量场展示

由于都是区域高精度展示，绘制的逻辑与上面的逻辑相仿，都是先划定区域，再提高精度的方法。只不过在流场绘制中需要控制的精度为density 。

def plot_selected_area(x1, y1, x2, y2):
    fig2 = plt.figure(figsize=(12, 6))
    ax2 = fig2.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    gl2 = ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
    gl2.top_labels = False
    gl2.right_labels = False

    ax2.coastlines()
    ax2.set_extent([x1, x2, y1, y2], crs=ccrs.PlateCarree())

    ax2.streamplot(lon, lat, u, v, transform=ccrs.PlateCarree(), density=1)
    plt.show()

本次高精度展示直接采用了Python自带的绘制流线方法，因为该方法可以根据划定区域自动重新计算流线密度，效果更好。

• 成果图：

原图

选择区域1

高精度1

选择区域2

高精度2

选择区域3

高精度3

四、可视化故事

• 人口数据

  我们先查看中国区域的人口数据热度图，可以发现，中国的人口分布大部分集中在东部或中部地区。其中，东北地区，人数的分布较为平均，且人口密度普遍较大；东南地区沿海地区，可以看到一个明显的趋势，即有几个城市人口的热度颜色突出，但周边城市热度颜色较为平淡，因此呈现出蛛网型分布。

  

  除此之外，可以发现，大部分城市的人口分布都是呈现一点尤为突出，然后向周边散开的形势，这说明当今世界的人口分布还是以大城市周边都市圈为主。而也有一些逆城市化的国家，如英国、德国，比起中国、美国或者日本则人数分布更加均匀一些。

• 洋流

  洋流对照图

  首先我们可以查看一下全局的洋流分布，总体上，与资料图较为匹配。包括两道赤道暖流和赤道逆流的分布、西风漂流等等。

  洋流是由大海中的各种力量所驱动形成的。其中，风力和水的密度差异是影响洋流形成和演变的主要因素。例如，赤道附近的东西向风力会产生压强差，从而推动海水向东流动，形成了赤道逆流。

  西风带位于南极洲和南半球的海洋之间，是全球最强劲的风带之一。在这里，强劲的西风把海水往东吹，形成了一个巨大的涡旋，从而形成了西风漂流。