【GIS】夜光遥感数据时序分析

基于夜光遥感的旅游产业发展分析

---

夜光遥感技术是指遥感平台获取夜间灯光数据的遥感技术。其具有空间清晰的与城市化有关的社会经济活动强度的感知信息,为探索城市和城市化问题提供了新的研究途径。除却反映城市夜间灯光活动，夜光遥感数据在监测森林火灾、反映人类活动、衡量经济发展等诸多领域也崭露头角，展现出其独特的价值。

本文基于格洛纳斯(GLONASS)卫星导航系统获取到的2018-2021年夜光遥感数据，分析丽江市旅游产业在疫情时代下的发展趋势。(注：该数据已经过预处理，如果是珞珈数据或是其他数据格式可能需要根据其提供的帮助文档进行预处理)

在开始之前，我们先介绍一些常见的夜光遥感数据的获取途径：

• 1992年——2013年DMSP-OLS数据：https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html

• 2012年——2019年NPP/VIIRS数据：https://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html

• 2010年DMSP-OLS数据、2013年——2019年NPP/VIIRS数据：https://www.lightpollutionmap.info

• 珞珈一号：

  http://59.175.109.173:8888/app/login.html

---

一、项目准备工作

1.1 研究数据

在本项目中，我们使用的数据有：

• 2018年9月-2018年12月 云南丽水逐月平均夜光遥感影像图
• 2019年1月-2019年12月 云南丽水逐月平均夜光遥感影像图
• 2020年1月-2020年12月 云南丽水逐月平均夜光遥感影像图
• 2021年1月-2021年12月 云南丽水逐月平均夜光遥感影像图
• 丽江旅游产业相关POI的Excel表数据

夜光影像图的空间分辨率为500m。

1.2 环境搭建与数据预处理

为了方便工作的进行，我们将数据统一放入一个文件夹中，并在此新建一个新的文件地理数据库。

我们以2018年十月的数据为例，进行数据探查。在ArcGIS中，其显示为

该数据的值域范围为[0.06,181.83]，以float栅格形式存储。为了更加清晰地表征数据，我们可以对数据的色彩直方图做一定的线性拉伸，该操作能够更好地分配区间内色彩值的映射，也就是更好的可视化表达。直方图均衡化也常见于CV中灰度图像的处理，有利于更好地挖掘图像的潜在特征表达。

在直方图均衡化后，结果相对清晰、分明。

遥感图像由于受到噪声点影响，会出现负值的情况，本文将负值统一修正为0，以便满足其在语义上的不变性。其操作可以在栅格计算器中完成。

处理完栅格数据之后，我们需要加载兴趣点数据。本文获取的POI是以Excel表的形式存储，其数据结构如下：

在Excel转表工具中，可以将其转为关系表：

接着将其作为点shp添加进文件数据库中：

至此，我们已经完成了初步的环境搭建与数据处理部分。

---

二、评价体系

为了评估不同时期夜光数据的变化，本文采用了环比增长率和同比增长率进行处理。其公式表达如下：
$$S_i=\frac{S_i-S_{i-1}}{S_{i-1}}$$
其中，计算同比增长率时$i$表示年度数据，计算环比增长率时$j$表示月度数据。

对于增长率为正值的区域，可认为该区域处于发展状态，而对于增长率为负值的区域，则认为该区域处于萧条状态。

---

三、实验操作

3.1 增长率计算

环比与同比的计算方法是一样的，这里我们拿同比增长率来举例。计算同比增长率需要相隔一年的同月数据，这里我们选择2019年10月和2020年10月的数据进行计算。将数据加载进ArcMap后，选择栅格计算器，输入公式1。

于是我们得到了一张同比增长率的图。关于这幅图我们等等细说，为了方便研究增长率，我们可以创建二值图像，并对信息做一些统计。

创建好的二值图像如下图所示，其中绿色表示增长率为负数的区域，粉色则是增长率为正数的区域，值表示对应栅格个数。

虽然该方法不能评估整体的经济发展状态，但可以得出大部分地区相较于去年，都呈积极发展态势。再来看看下面这张图，结合同比增长率变化图，我们可以发现，左下高度发展区域在一年后反而呈现衰退趋势，而右上和右下两个高度发展区域发展态势依旧良好。

当然，这些图都是在表征整个研究区域的情况，我们的目的是研究旅游变化趋势，因而需要进一步细化研究区域。

3.2 旅游焦点区域

前面提到本文初步计算得到的结果是全局的结果，为了追踪旅游行业的变化趋势，我们选取了带有相关语义的POI数据，这些数据表现为与旅游产业呈现较高的相关性。在前面的数据预处理中，我们将这些数据进行了格式转换，并导入了ArcMap中。

仅有POI还是很难体现空间上的相关性，因此，我们为这些POI增设一个服务范围。该服务范围具有一个服务搜索半径，该半径可以是固定值，也可以是自适应增长的。关于自适应增长的搜索半径接下来我也会发一篇文章重点讲讲，本文选择了较为简单的固定方法。

此时我们选择三个固定阈值作为POI的搜索半径，分别是$1KM,3KM,5KM$。借助这些半径值构建的环形缓冲区，我们将其称为兴趣点的服务范围。实际含义为服务场所所能提供相关服务的影响范围。

我们以景区的POI数据为例，在ArcMap中构建$3KM$的缓冲区。

由于距离原因，难免会出现缓冲区上的叠置。但本篇文章的论述重点并不在某个景点或是某个酒店，故本文将所有的相关语义点缓冲区视为相同的共享对象。换而言之，就是把带有旅游语义的兴趣点视为空间上的一个点。(可以简单理解为图论中的超级源点)

我们首先通过按掩膜提取的方式获取焦点区域的数据值：

接下来可以通过符号系统里面的分类方式简单获取数据的统计信息。

3.3 结果可视化

由于在ArcMap中，栅格数据矢量表的构建仅限于单波段整形数据，本文为了保证数据的两位精度，先将栅格扩大$1e2$倍，再转换为整形栅格：

此时栅格数据已经构建了保留两位精度的属性表，我们将其导出为csv格式的数据，再进行处理。

我们采用高级箱线图展现数据，该箱线图能够在各个采样区间评估数据，越宽所代表的样本量就越多。

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
path=r"C:\Export_Output.csv"
data=pd.read_csv(path)

# 此时需要构建一个统计信息
# 会占用大量空间
collect=[]

for i in range(data.shape[0]):
    collect.extend([data.iloc[i,1]/100.0]*int(data.iloc[i,2]))

# 再通过seaborn绘制盒须图
new_data=pd.DataFrame(collect,columns=["A"])
sns.set_theme(style="whitegrid")
ax=sns.boxenplot(y=new_data["A"])
ax.set_title("2020_10 Night Value")
ax.set_ylabel("Value")
plt.show()

或者是小提琴图:

ax=sns.violinplot(y="A",data=new_data)
sns.despine(left=True)

当然，本文最主要的还是做比对，除却小提琴图外，我们依旧需要绘制总量上的一个变化和变化率曲线图。(以下是环比数据的比较)

通过比对数据，我们发现在2020年10月到12月三个月中，带有旅游语义的夜光遥感总值呈现下降、上升趋势，而最发达的地区进一步发展，大量欠发展地区进入更高的发展阶段。11月可能受到疫情和政策的影响，旅游产业总体发展呈现下滑趋势，但低发展区域更多地进入了新一发展阶段，而最高值仍有提升。中间段反而呈现下降趋势，中小旅游产业受到疫情的冲击较大，而微型和大型旅游产业依旧稳固。

---

四、最终结果展示

各季度夜光数据分布图

各季度亮度总值

各季度环比增长率

各季度同比增长率

可视化代码

import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import argparse

def parse_args():
    parser=argparse.ArgumentParser("Read CSV")
    parser.add_argument("-root","--r",help="The root of your csv-files",default=r"..\csv",type=str)
    return parser.parse_args()

def main(args):
    path=args.root
    
    list_path=[ i for i in os.listdir(path) if i .endswith(".csv")]
    
    collect=[]
    sum_list=[]
    label=[]
    
    for lp in list_path:
        data=pd.read_csv(path+"\\"+lp)
        tem=[]
        for i in range(data.shape[0]):
            tem.extend([data.iloc[i, 1] / 100.0] * int(data.iloc[i, 2]))
            label.extend([lp.split(".")[0]] * int(data.iloc[i, 2]))
        sum_list.append(sum(tem))
        collect.extend(tem)
    
    data=np.array([collect,label]).T
    new_data=pd.DataFrame(data,columns=["A","B"])
    new_data['A'] = new_data['A'].astype('float')
    sns.set_theme(style="whitegrid")
    # ax=sns.boxenplot(x="B",y="A",data=new_data)
    new_data= new_data.explode('A')
    
    ax=sns.violinplot(x="B",y="A",data=new_data)
    sns.despine(left=True)
    ax.set_title("Night Values")
    ax.set_ylabel("Value")
    ax.set_xlabel("Times")
    plt.show()
    
    plt.plot([i.split(".")[0] for i in list_path],sum_list)
    plt.xlabel("Times")
    plt.ylabel("Value")
    plt.title("Sum")
    plt.show()
    
    record=sum_list[2:]
    
    # 环比
    for i in range(1,len(sum_list)):
        sum_list[i-1]=(sum_list[i]-sum_list[i-1])/sum_list[i-1]
    sum_list=sum_list[:-1]
    plt.plot(range(len(sum_list)),sum_list)
    plt.xlabel("Times")
    plt.ylabel("rate")
    plt.title("Month-on-month Rate")
    plt.show()
    
    # 同比
    X=["First","Second","Third","Fourth"]
    data=[]
    for j in range(2):
        tem=[]
        for i in range(j*4,j*4+4):
            tem.append((record[i+4]-record[i])/record[i])
        data.append(tem)
    x=range(4)
    plt.plot(x,data[0],"r",label="2019-2020")
    plt.plot(x,data[1],"b",label="2020-2021")
    plt.legend()
    plt.xlabel("Quarter")
    plt.ylabel("Rate")
    plt.xticks(x,X)
    plt.title("Year-on-Year")
    plt.show()

if __name__ == '__main__':
    args=parse_args()
    main(args)

---

五、时空格局分析

为了分析夜光强度的空间格局变化，我们先给夜光数据统一量纲标准。在重分类工具中，我们将夜光强度以每30个标准单位进行划分。

由于重心分析需要矢量数据，所以我们利用按属性提取工具将重分类值为2的数据提取出来(每个都要提取，这里举例)

利用栅格转面工具将其转为shp图。

当然，由于栅格大小的原因，在矢量化的时候会丢失较多的信息，可以通过重采样的方式增加栅格的分辨率，或者是通过GAN生成高分辨率图像。

在空间统计工具中，找到标准差椭圆工具。

为其绘制标准差椭圆。

此时，我们就得到了2018年第四季度中等强度的空间分布格局，我们再对2021年第四季度做同样的事情：

该图便显现了三年中，中等强度地区的发展倾向。原先的重心开始往东移动，且更加聚集。