在模仿中精进数据分析与可视化01——颗粒物浓度时空变化趋势（Mann–Kendall Test）

本文是在模仿中精进数据分析与可视化系列的第一期——颗粒物浓度时空变化趋势（Mann–Kendall Test），主要目的是参考其他作品模仿学习进而提高数据分析与可视化的能力，如果有问题和建议，欢迎在评论区指出。若有其他想要看的作品，也欢迎在评论区留言或通过公众号发消息给出相关信息。

扫描下面的二维码，关注公众号GeodataAnalysis，搜索本篇文章获取文中示例数据与代码：

一、简介

本次要模仿的作品来自论文Investigating the Impacts of Urbanization on PM2.5 Pollution in the Yangtze River Delta of China: A Spatial Panel Data Approach，研究区域为上海、安徽、浙江和江苏，所用数据为 2002–2017该区域PM2.5浓度栅格数据，数据来源于 Dalhousie University Atmospheric Composition Analysis Group开发的年均PM2.5数据集V4.CH.03，空间分辨率为0.01°×0.01°（原论文采用数据的空间分辨率为1km×1km，但我在该网站上找不到，可能是不提供下载了）。

二、数据下载和处理

数据下载格式为.asc，使用gdal将其转为.tif格式，所用代码如下。

# -*- coding: utf-8 -*-

import os
from osgeo import gdal

inpath = "./ASCII" #待转换的栅格的存储路径，会转换该路径下的所有栅格
outpath = "./TIF" #输出栅格的路径，最好是空路径

def ASCII2TIF(filepath, outfilepath):
    ds = gdal.Open(filepath)
    tif = gdal.Translate(outfilepath, ds)
    tif = None

print("Starting Convert!")
for filename in os.listdir(inpath):
    if filename.endswith(".asc"):
        filepath = os.path.join(inpath, filename)
        outfilepath = os.path.join(outpath, filename.replace(".asc", ".tif"))
        ASCII2TIF(filepath, outfilepath)

print("Convert Over!")

三、Mann–Kendall趋势分析

Mann–Kendall趋势分析的具体计算方法这里不再赘述，原文作者采用R语言的trend package计算的，本文采用python的pymannkendall计算，github项目地址为https://github.com/mmhs013/pyMannKendall。

原文的趋势分析包括两部分，一部分是计算slope值，slope值为正，则表明具有上升的趋势，反之亦然；另一部分是计算p值，p值越小趋势越显著，0.01<p<0.05说明趋势显著，p<0.01说明趋势非常显著。这两部分的结果都可以使用pymannkendall的original_test函数计算，pymannkendall的简单用法介绍如下。

A quick example of pyMannKendall usage is given below. Several more examples are provided here.

import numpy as np
import pymannkendall as mk

# Data generation for analysis
data = np.random.rand(360, 1)

result = mk.original_test(data)
print(result)

Output are like this:

Mann_Kendall_Test(trend='no trend', h=False, p=0.9507221701045581, z=0.06179991635055463, Tau=0.0021974620860414733, s=142.0, var_s=5205500.0, slope=1.0353584906597959e-05, intercept=0.5232692553379981)

Whereas, the output is a named tuple, so you can call by name for specific result:

print(result.slope)

or, you can directly unpack your results like this:

trend, h, p, z, Tau, s, var_s, slope, intercept = mk.original_test(data)

四、计算并保存结果

这里依然使用python作为分析计算的工具，所用代码如下所示。

由于 pymannkendall较为臃肿，计算速度很慢，并且暂不支持numba加速，有需要大量计算的可根据其源码重新编写函数，实现numba加速，如本文的get_slope函数，在使用numba加速后计算pvalues仅需4秒，使用pymannkendall的original_test则需要几分钟的时间。

# -*- coding: utf-8 -*-
import os
import numba
import numpy as np
from glob import glob
import rasterio as rio
import geopandas as gpd
import pymannkendall as mk
from rasterio.mask import mask

inpath = r"./TIF"  # .tif文件的保存路径
p_path = r"./pvalues.tif"  # p-values的输出路径
slope_path = r"./slopes.tif"  # slopes的输出路径
trend_path = r"./trends.tif"  # 原图左图中不同的趋势
border_path = r"./Shapefiles/border.shp"  # 研究区域

# 获取2002-2017年的栅格数据的路径
def get_raster_paths(inpath):
    paths = []
    for year in range(2002, 2018):
        year_path = glob(os.path.join(inpath, "*"+str(year)+"*.tif"))
        if year_path:
            paths.append(year_path[0])
        else:
            print("can't find raster of {} year!".format(year))
    return paths

# 裁剪栅格，并将结果转为numpy数组
def clip_raster_to_array(paths, border, dtype='float64'):
    src = rio.open(paths[0])
    gdf = gpd.read_file(border)
    array, gt = mask(src, [gdf.dissolve().geometry.__geo_interface__],
                     crop=True, nodata=src.nodata)
    array[array == src.nodata] = np.NAN
    arrays = np.full(shape=(len(paths), array.shape[0], array.shape[1]),
                     fill_value=np.NAN, dtype=dtype)
    arrays[0] = array

    for i in range(1, len(paths)):
        src = rio.open(paths[i])
        array, gt = mask(src, [gdf.dissolve().geometry.__geo_interface__],
                        crop=True, nodata=src.nodata)
        array[array == src.nodata] = np.NAN
        arrays[i] = array

    return arrays, src.crs, src.transform

# 使用rasterio将numpy数组转为tif
def array_to_raster(path, array, crs, transform, nodata=None):
    with rio.open(path, 'w', driver='GTiff', nodata=nodata,
                  height=array.shape[0], width=array.shape[1],
                  count=1, dtype=array.dtype, crs=crs,
                  transform=transform, compress='lzw') as dst:
        dst.write(array, 1)

# 计算slope值，使用numba加速
@numba.njit
def get_slope(x):
    if np.isnan(x).any():
        return np.NAN
    idx = 0
    n = len(x)
    d = np.ones(int(n*(n-1)/2))

    for i in range(n-1):
        j = np.arange(i+1, n)
        d[idx: idx + len(j)] = (x[j] - x[i]) / (j - i)
        idx = idx + len(j)

    return np.median(d)

# 计算slope和p值
def mk_(x):
    if np.isnan(x).any():
        return np.NAN
    result = mk.original_test(x.data)
    return result.slope, result.p


paths = get_raster_paths(inpath)
arrays, crs, gt = clip_raster_to_array(paths, border_path)
print("clip raster to array over!")

print("calculate slope and p-value over!")
slopes, pvalues = np.apply_along_axis(mk_, 0, arrays)
# 使用numba加速的示例
# slopes = np.apply_along_axis(get_slope, 0, arrays)

# 计算有显著和非常显著趋势的区域
trends = np.full(shape=slopes.shape, fill_value=0, dtype=np.int8)
trends[~np.isnan(slopes)] = 0  # nodata值及不显著区域设为0
# 比较显著增加的区域设为1
trends[(slopes > 0) & ((0.01 < pvalues) & (pvalues < 0.05))] = 1
# 显著增加的区域设为2
trends[(slopes > 0) & (pvalues < 0.01)] = 2
# 比较显著减少的区域设为3
trends[(slopes < 0) & ((0.01 < pvalues) & (pvalues < 0.05))] = 3
# 显著减少的区域设为4
trends[(slopes < 0) & (pvalues < 0.01)] = 4

# 保存栅格
array_to_raster(p_path, pvalues, crs, gt, np.NAN)
array_to_raster(slope_path, slopes, crs, gt, np.NAN)
array_to_raster(trend_path, trends, crs, gt)
print("save rasters over!")

五、结果绘图

由于QGIS软件打开和一些相关操作的速度都要比ArcGIS快的多，而且QGIS内置的取色器的功能也方便绘图时设置颜色，因此本文使用QGIS绘制结果图，如下图所示。