MATLAB与Arcgis技术层面实现城市灯光中心转移可视化

一、MATLAB实现夜间灯光数据的灯光中心获取

普遍来说，城市中心被认为是与周边地区相比，拥有较高人口密度、就业规模和经济聚集度的城市区域（Giuliano and Small，1991），是一座城市经济文明的聚集区，是各种社会经济要素在空间上集聚的结果。总的来说，城市中心就是人类活动比周围都要活跃的城市区域（Tammaru，2001），表现在夜间灯光影像上就是高值像元（Chen et al.，2017），即在一定的城市空间范围内，该像元的夜间灯光亮度值比周围的所有像元都要大[6]。

目前大尺度空间的城市中心区域空间边界的确定主要依靠统一的阈值方法，下述方法提取的是中心点[6]。（值得注意的是，根据生活经验来说，独立工厂通常也具有较大的光亮值并且其通常可能被误提取为中心）

（一）MATLAB代码

由于遥感影像本身具备图像值矩阵与坐标元素结构体两部分的复杂性以及研究区四周几乎不存在中心出现的可能性，暂不考虑边界像素的遍历问题，多次采用不同边长方形核实验后，采用25*25核进行遍历。

%获取文件绝对路径
File = 'D:\matlab\code\Plight\20202.tif';
%本处提供核大小更改交互（仅支持奇数边长方形核）
len = input('请输入奇数边方形核边长：');
ck = (len-1)/2;
%读取tiff文件的影像矩阵与坐标系文件
[A,R] = geotiffread(File);
%获取投影坐标系
info = geotiffinfo('D:\matlab\code\Plight\20202.tif');
[r,c] = size(A);
%创建一个与影像矩阵同大小的outPut矩阵
outPut = ones(size(A));
%将output矩阵值全赋为-1，避免与原数据相减时误提取灯光中心
outputR = outPut .*(-1);
%处理过程
for i = (ck+1):(r-ck)
    for j = (ck+1):(c-ck)
        MatC = A((i-ck):(i+ck),(j-ck):(j+ck));
        MatMax = max(MatC(:));
        for k = (i-ck):(i+ck)
            for l = (j-ck):(j+ck)
                outputR(k,l) = MatMax;
            end
        end
    end
end
%输出图像
fprintf('数据处理完毕，正在导出！\n');
geotiffwrite('centralLight2.tif',outputR,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
fprintf('数据导出完毕！');

（二）灯光中心输出影像与原始影像做差值(可基于ENVI中的Band Math来做)

(abs(灯光中心影像-原始影像) lt 0.000009)*原始影像

差值影像中像元值为0（除背景值）的像元即为灯光中心（此处需要考虑到浮点数运算必须有个误差限比较）。