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

最新推荐文章于 2024-05-07 12:38:08 发布

doll ～CJ

最新推荐文章于 2024-05-07 12:38:08 发布

阅读量1.8k

点赞数 8

分类专栏：数字图像处理文章标签： python 科技

本文链接：https://blog.csdn.net/Jcb1906824038/article/details/126236457

版权

数字图像处理专栏收录该内容

11 篇文章 4 订阅

订阅专栏

该文介绍了通过突变检测法提取城市建成区的方法，结合GIS和Python的arcpy库，处理NPP-VIIRS和DMSP-OLS夜间灯光数据。首先，使用Mann-Kendall突变点检测法确定最佳阈值，然后通过arcpy进行数据处理，包括栅格复制、重分类、属性表导出和Excel转换。最后，应用提取的阈值识别建成区。整个过程涉及文件操作、数据统计分析及地图自动化。

摘要由CSDN通过智能技术生成

一、突变检测法[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表转xls
    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::] == '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]
    # 正序列计算---------------------------------
    # 定义累计量序列Sk
    Sk = []
    # 定义统计量UFk
    UFk = []
    # 定义Sk序列元素s
    s = 0
    Exp_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 + 1
            else:
                s = s + 0
        Sk.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 = 0
    Exp_value2 = []
    Var_value2 = []
    # 按时间序列逆转样本y
    inputdataT = list(reversed(inputdata))
    # 此时UBk无意义，因此公式中，令UBk(1)=0
    for i in range(0, n):
        for j in range(i):
            if inputdataT[i][1] > inputdataT[j][1]:
                s2 = s2 + 1
            else:
                s2 = s2 + 0
        Sk2.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 K

def 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 savepath

tiffcopy(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博客