一、突变检测法[8]

通过确定基于突变检测原则的单个阈值，该阈值既有显著降低灯光数据的噪声以及沿海城市的灯光溢出的功能，同时也能保留具有连续灯光亮度值的多边形即较大的城市区域。对于分离城市灯光区与非城市灯光区的最佳阈值，此处经验性地认为其能够保持城市核心区完好无损。

突变检测法利用了一种GIS方法，引入了周长和面积参数，计算某个阈值二值化后的影像（大于灯光阈值部分）覆盖范围内所有多边形的周长与面积。在基于该阈值二值化后，我们明显可以观察到覆盖区域多边形内部出现碎片，同时所有多边形的周长测量总和有一个大幅的增加同时面积减少。通常城市建成区“内部一定程度破碎”前的阈值被认为是该研究区提取城市建成区的最佳阈值。

原文献是通过尝试了多个城市内部破裂前的阈值并采取平均的数学方法确定了推广阈值。

二、arcpy库包简概[1]

arcpy是一个必须依附于Arcgis/GeoScene软件的Python库包，同时Arcgis或Arcgis Pro的下载都会携带一个内置Python。在Pycharm中导入arcpy库包需要以arcgis内置Python版本作为解释器。

Arcpy可提供如下能力：

1、以实用高效的方式通过Python执行地理数据分析、数据转换、数据管理和地图自动化；

2、快速调用Arcgis/GeoScene提供的地理处理工具以及其他函数、类和模块、并且可以创建更加灵活可控的工作流；

3、访问ArcGIS/GeoScene软件生成的工程内容、资源、配置，以及进行操作处理；

4、利用第三方Python工具包，与GIS软件尽心互补；

5、封装自定义的脚本为工具包，发布给其他人使用。

在学习arcpy的过程中，我们需要对arcgis的功能使用有一定的了解，如下提供一种较为有效的查阅arcpy使用学习的两种方式：

一、我们可以首先通过构思整体处理流程框架，然后依次去寻找tools box里的工具点击工具帮助，然后查阅arcpy代码示例进行学习；

二、查阅[12][13]。

三、OS库包简概[5][6]

os模块是python标准库中的一个用于访问操作系统功能的模块，os模块提供了其他操作系统接口，可以实现跨平台访问。

四、实现代码

代码采用的是Mann-Kendall突变点检测法（曼-肯德尔法），具体原理详见[9][10]，Python代码参考[2]。（需要注意的是，我们在为文件命名时应保证良好的命名习惯，如以字母开头后接字符）

此段代码可以修改相对路径，将各个函数整合为一个函数，后续改进空间较大。

# encoding:UTF-8
# 编码声明import arcpy
from arcpy.sa import *
import os
import pandas as pd
import numpy as np# 本次灯光数据源自美国科罗拉多矿业大学制作的NPP-VIIRS月合成数据以及DMSP-OLS年合成数据，并且数据已经过
# 年内校正、年际校正、传感器相互校正、去除异常值操作# tiff文件复制导出
# 一般来说，在ENVI获取到的TIFF文件在arcgis中数据图层加入的时候会缺失属性源的统计值
# 所以我们直接通过GetRasterProperties_management（）函数获取统计值时会产生error“没有可用的统计数据而失败”
# 解决办法即为重新导出文件或者再次复制数据，此时其数据属性源就能够携带统计数据
# rootfile为tiff文件存储的根目录，savepath为自己创建的文件夹目录
def tiffcopy(rootfile, savepath):# 为复制导出的栅格tiff文件创建专属存储文件夹if not os.path.exists(savepath):os.mkdir(savepath)# 批量文件复制操作for root, dirs, files in os.walk(rootfile):for f in files:path = os.path.join(root, f)# 以绝对路径字符串后缀名为tif判断tif文件if f[len(f)-3::] == 'tif':copypath = os.path.join(savepath, f)print(copypath)arcpy.CopyRaster_management(in_raster=path, out_rasterdataset=copypath, format="TIFF")print("tiff文件导出结束！")return savepath# 浮点数数据按规定间隔重分类为整数并导出像元值统计
def datastatistic(rootfile, savepath):# 为重分类栅格tiff文件创建专属存储文件夹if not os.path.exists(savepath):os.mkdir(savepath)# 根据tiff栅格文件统计像元最小值、最大值，定义重分类间隔列表# 批量文件操作for root, dirs, files in os.walk(rootfile):for f in files:if f[len(f)-3::] == 'tif' and root == rootfile:# 记录tiff文件的绝对路径absfpath = os.path.join(root, f)# 开始获取tiff栅格文件像元最大值和最小值(实际上预处理过后的灯光值一般最小均为0)Rmaxinfo = arcpy.GetRasterProperties_management(absfpath, "MAXIMUM")dmax = int(eval(Rmaxinfo.getOutput(0))+1)Rmininfo = arcpy.GetRasterProperties_management(absfpath, "MINIMUM")dmin = int(eval(Rmininfo.getOutput(0)))print("dmax:{0}     dmin:{1}".format(dmax, dmin))# 对数据构造重分类间隔列表originalSTR = "["for i in range(dmin, dmax):if i == dmax - 1:originalSTR = originalSTR + "[{0},{1},{2}]]".format(dmax - 1, dmax, dmax)else:originalSTR = originalSTR + "[{0},{1},{2}],".format(i, i+1, i+1)# 生成分割列表myRemapRange = RemapRange(eval(originalSTR))outputTIFF = Reclassify(absfpath, "VALUE", myRemapRange)# 保存重分类文件path = os.path.join(savepath, f)outputTIFF.save(path)print(path)print("批量TIFF数据重分类已结束!")return savepath# 对重分类数据将属性表dbf导出
def recCount(rootfile, savepath):# 为dbf数据创建专属存储文件夹if not os.path.exists(savepath):os.mkdir(savepath)# 首先批量导出栅格属性表（利用表至表的方法由栅格数据提取栅格属性表为.dbf）for root, dirs, files in os.walk(rootfile):for f in files:path = os.path.join(root, f)# 以绝对路径字符串后缀名为tif判断tif文件if f[len(f) - 3::] == 'tif':fname = f[:-4] + ".dbf"arcpy.TableToTable_conversion(path, savepath, fname)print(savepath + fname)print("批量dbf属性表数据转换完成！")return savepath# dbf数据表转换xls
def dbf2xls(rootfile, savepath):# 为excel数据创建专属存储文件夹if not os.path.exists(savepath):os.mkdir(savepath)# 首先批量dbf表转xlsfor root, dirs, files in os.walk(rootfile):for f in files:path = os.path.join(root, f)# 以绝对路径字符串后缀名为tif判断tif文件if f[len(f) - 3::] == 'dbf':fname = f[:-4] + ".xls"outpath = os.path.join(savepath, fname)arcpy.TableToExcel_conversion(path, outpath)print(outpath)print("批量dbf表格转excel表已完成！")return savepath# Mann-Kendall突变检测法
def Kendall_change_point_detection(inputdata):# numpy.ndarray数据结构.shape函数作用是# 读取矩阵的长度，或矩阵在某一维上的长度。shape[0]返回矩阵的行数,shape[1]返回矩阵的列数n = inputdata.shape[0]# 正序列计算---------------------------------# 定义累计量序列SkSk = []# 定义统计量UFkUFk = []# 定义Sk序列元素ss = 0Exp_value = []Var_value = []# 因为根据统计量UFk公式，i=0时，Sk(0)、E(0)、Var(0)均为0,此处注意我们在重分类时认为0-1区间归为1，因此我们仍考虑从i=0开始# 此时UFk无意义，因此公式中，令UFk(0)=0# 正序列计算for i in range(0, n):for j in range(i):if inputdata[i][1] > inputdata[j][1]:s = s + 1else:s = s + 0Sk.append(s)# Sk[i]的均值Exp_value.append(((i + 1) * (i + 2) / 4))# Sk[i]的方差Var_value.append((i + 2) * (i + 1) * (2 * (i + 2) + 5) / 72)UFk.append((Sk[i] - Exp_value[i]) / np.sqrt(Var_value[i]))# 逆序列计算# 定义逆序累计量序列Sk2，长度与inputdata一致Sk2 = []# 定义逆序统计量UBk，长度与inputdata一致UBk = []UBk2 = []s2 = 0Exp_value2 = []Var_value2 = []# 按时间序列逆转样本yinputdataT = list(reversed(inputdata))# 此时UBk无意义，因此公式中，令UBk(1)=0for i in range(0, n):for j in range(i):if inputdataT[i][1] > inputdataT[j][1]:s2 = s2 + 1else:s2 = s2 + 0Sk2.append(s2)# Sk[i]的均值Exp_value2.append((i + 1) * (i + 2) / 4)# Sk[i]的方差Var_value2.append((i + 1) * (i + 2) * (2 * (i + 2) + 5) / 72)UBk.append((Sk2[i] - Exp_value2[i]) / np.sqrt(Var_value2[i]))UBk2.append(-UBk[i])# 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，# 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反# 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))# 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势#  UBk(i)=0-(Sk2(i)-E)/sqrt(Var)# ------------------------------逆序列计算# 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量# 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此# 再按时间序列逆转结果统计量UBk，得到时间正序的UBkT，UBkT = list(reversed(UBk2))diff = np.array(UFk) - np.array(UBkT)K = list()# 找出交叉点for k in range(0, n-1):if diff[k] * diff[k + 1] < 0:K.append(k)# 返回序列号，需要在mk函数中转换对应的DN值return Kdef mk(rootfile, savepath):# 创建正反交点图专属文件夹if not os.path.exists(savepath):os.mkdir(savepath)# 计算分割阈值函数主体list_MK = []for name in os.listdir(rootfile):if name.endswith('.xls'):name_path = os.path.join(rootfile, name)# 只读取excel中的DN值像元计数列以及对应阈值excel = pd.read_excel(name_path, usecols=[1, 2])df = pd.DataFrame(excel)# pandas中 head()函数只能读取前五行数据# df.drop()作用为从行或列中删除指定的标签。# df.drop()把Value字段删除了，实际上其与直接按字段读取结果一致df1 = df.drop('Value', axis=1)# numpy和pandas数据结构转换,<class “pandas.core.frame.DataFrame”>转换为<numpy.ndarray>inputdata = np.array(df)K = Kendall_change_point_detection(inputdata)list_MK.append(K[1])# 阈值列表获取代码函数MK = []for i in list_MK:MK.append(inputdata[i][0])print("时序分割阈值如下:")print (MK)return MK# 分割阈值导出建成区
# 此处的rootfile为存储重分类栅格的文件夹,dnList为阈值提取列表
def extract(rootfile, dnList, savepath):# 创建文件夹存储建成区数据if not os.path.exists(savepath):os.mkdir(savepath)# 将栅格数据raster与提取阈值对应成字典rasters = []for root, dirs, files in os.walk(rootfile):for f in files:path = os.path.join(root, f)# 以绝对路径字符串后缀名为tif判断tif文件if f[len(f) - 3::] == 'tif' and root == rootfile:rasters.append(path)# zip函数可详见参考资料[15]# zip()函数用于将可迭代的对象作为参数，将对象中对应的元素打包成一个个元组，然后返回由这些元组组成的列表DNdicts = dict(zip(rasters, dnList))# 开始批量按属性提取for raster in rasters:sql = "Value >= {}".format(DNdicts[raster])urban = ExtractByAttributes(raster, sql)path = os.path.join(savepath, raster[-9::])urban.save(path)print(path)print("批量提取建成区已完成！")return savepathtiffcopy(r"G:\esri\arcpyEXP\data", r"G:\esri\arcpyEXP\re\newtiff")
datastatistic(r"G:\esri\arcpyEXP\re\newtiff", r"G:\esri\arcpyEXP\re\newtiff\reclass")
recCount(r"G:\esri\arcpyEXP\re\newtiff\reclass", r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf")
dbf2xls(r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf", r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\xls")
mk(r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\xls", r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\xls\graphic")
extract(r"G:\esri\arcpyEXP\re\newtiff\reclass", mk(r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\xls", r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\xls\graphic"), r"G:\esri\arcpyEXP\re\newtiff\reclass\dbf\urban")

参考资料：

[1][ArcPy百科]第一节：什么是ArcPy

[2]python处理批量根据夜间灯光数据采用突变检验法提取建成区_哔哩哔哩_bilibili

[3]在PyCharm中导入和使用arcpy_小云要努力的博客-CSDN博客_pycharm导入arcpy

[4]ArcGIS Pro Python 参考—ArcGIS Pro | 文档

[5]Python入门之os模块详解_ly2020_的博客-CSDN博客_python中os是什么模块

[6]python标准库 —— os模块_wakeyo_J的博客-CSDN博客_python的os库

[7]余柏蒗，王丛笑，宫文康，陈佐旗，施开放，吴宾，洪宇辰，李乔玄，吴健平.2021.夜间灯光遥感与城市问题研究：数据、方法、应用和展望.遥感学报，25（1）：342-364.

[8]Imhoff M L, Lawrence W T, Stutzer D C and Elvidge C D. 1997. A technique for using composite DMSP/OLS“city lights”satellite data to map urban area. Remote Sensing of Environment, 61(3): 361-370 [DOI: 10.1016/S0034-4257(97)00046-1.

[9]【4.4 Mann-Kendall检验】百度文库

[10]【曼－肯德尔法_百度百科】百度百科

[11]【MATLAB Mann-Kendall突变检验 (mk突变检验)】CSDN编程社区

[12]arcgis 10.2 arcpy帮助文档

[13]ArcGIS 帮助 10.1（网页版帮助文档）

[14]（二十三）arcpy开发&利用GetRasterProperties_management获取栅格数据相关信息_yGIS的博客-CSDN博客

[15]Python os.work()函数_万wu皆可爱的博客-CSDN博客_os.work python（os库分享优秀博主）

[16]Python之glob.glob() 函数解析_lindergarten的博客-CSDN博客_glob.glob() 函数

[17]Python中zip()、zip(*zipped)、*zip()函数总结_何极光的博客-CSDN博客_python zipped

[18]pandas写入excel_Pandas Excel教程：如何读取和写入Excel文件_cumei1658的博客-CSDN博客

Pycharm中利用arcpy实现灯光遥感数据的建成区批量提取相关推荐

hopper_如何利用卫星收集的遥感数据轻松对蚱hopper中的站点进行建模
hopper 建筑学与数据科学 (Architectonics and Data Science) Understanding the site and topography are crucial ...
Suomi NPP VIIRS夜间灯光遥感数据简介与下载（一）——数据介绍，FTP下载与hdf5读取
失踪人口回归+爷青回,本篇介绍下Suomi NPP VIIRS夜间灯光遥感数据下载. 文章目录 1 夜间灯光遥感数据简介 2 夜间灯光遥感数据下载 1 年和月尺度产品下载 2 日尺度产品下载 3 NP ...
Suomi NPP VIIRS夜间灯光遥感数据简介与下载
一.夜间灯光遥感数据简介当前常用的夜间灯光遥感数据主要是两个卫星,一个是DMSP(Defense Meteorological Sate-llite Program)是美国国防部的极轨卫星计划,传感 ...
Sql Server 中利用游标对table 的数据进行分组统计式输出…
Sql Server 中利用游标对table 的数据进行分组统计式输出- Table Name: Tb_Color Create table Tb_Color(id int identity(1,1) ...
在pycharm中使用arcpy
相关软件配置介绍: windows10,64位 pycharm2020.1 64位 arcgis10.6 32位 anconda 4.10.1 64位实现目标:在pycharm中实现arcpy的调用 ...
npp夜光数据介绍 viirs_1.夜间灯光遥感数据概述
1.夜间灯光遥感数据概述: 夜间灯光遥感数据主要来自于由美国国防气象卫星搭载的可见光成像线性扫描业务系统(DMSP/OLS)和国家极轨卫星搭载的可见光近红外成像辐射仪(NPP/VIIRS)获取的夜间灯 ...
三、DMSP/OLS、NPP/VIIRS夜间灯光数据之建成区提取——阈值确定（2）
https://blog.csdn.net/weixin_44725365/article/details/112169151https://blog.csdn.net/weixin_44725365 ...
二、DMSP/OLS、NPP/VIIRS夜间灯光数据之建成区提取——阈值确定
一.DMSP/OLS夜间灯光数据之建成区提取--理论https://mp.csdn.net/mp_blog/creation/editor/112169085 一.前言前文提到基于DMSP/OLS夜 ...
在pycharm中利用labelme标注生成语义分割文件
源程序: https://github.com/wkentaro/labelme 可以在自己创建的环境下安装labelme或者重新创建一个虚拟环境安装labelme 1.重新创建环境:在py ...

Pycharm中利用arcpy实现灯光遥感数据的建成区批量提取