Python ArcPy实现批量计算多时相遥感影像的各项元平均值

作者:疯狂学习GIS 时间:2022-06-04 14:45:51 

在遥感应用中,我们经常需要对某一景遥感影像中的全部像元的像素值进行平均值求取——这一操作很好实现,基于ArcMap软件或者简单的Python代码就可以实现;但有时候,我们会需要结合同一地区、不同时相多景遥感影像,求取每一个像元全部时相中像素值的平均值——这一需求的实现较之前者就有些麻烦,本文对此加以介绍。

首先,我们来明确一下本文的具体需求。现有一个存储有大量.tif格式遥感影像的文件夹,其中每一个遥感影像的文件名中都包含有该图像的成像时间,如下图所示。且其中除了.tif格式的遥感影像文件外,还具有其它格式的文件。

Python ArcPy实现批量计算多时相遥感影像的各项元平均值

我们希望,对于同一年成像的遥感影像进行逐像元平均值的求取。例如,上图中具有2001年第185天成像、第193天成像、第201天成像……等等遥感影像8幅,每一幅都是这一年不同时间在同一空间位置的成像;同时,还有2005年不同时间成像的遥感影像9幅。我们希望,首先将2001年成像的8幅遥感影像加以逐像元平均值的求取,即求取每一个像元在这8景图像中像素值的平均;随后再对2005年成像的9幅遥感影像加以逐像元平均值的求取,以此类推。

明确了需求后,我们就可以开始具体的操作。首先,本文所需用到的代码如下。

# -*- coding: utf-8 -*-
"""
Created on Sat Apr 16 10:48:37 2022

@author: fkxxgis
"""

import arcpy
from arcpy.sa import *

tif_file_path="E:/LST/Data/MODIS/05_Resample/"
average_file_path="E:/LST/Data/MODIS/06_Average/"
arcpy.env.workspace=tif_file_path

tif_file_name=arcpy.ListRasters("*","tif")
tif_file_year=tif_file_name[0][0:4]
one_year_tif_list=[]
sum_pic=0

for tif_file in tif_file_name:
   if tif_file[0:4]==tif_file_year:
       one_year_tif_list.append(tif_file)
       tif_file_temp=tif_file
       if tif_file==tif_file_name[len(tif_file_name)-1]:
           pic_num=len(one_year_tif_list)
           for tif_file_new in one_year_tif_list:
               sum_pic=sum_pic+Raster(tif_file_new)
           (sum_pic/pic_num).save(average_file_path+tif_file_year+"_Ave.tif")
   else:
       pic_num=len(one_year_tif_list)
       for tif_file_new in one_year_tif_list:
           sum_pic=sum_pic+Raster(tif_file_new)
       (sum_pic/pic_num).save(average_file_path+tif_file_year+"_Ave.tif")
       one_year_tif_list=[]
       sum_pic=0
       one_year_tif_list.append(tif_file)
       tif_file_year=tif_file[0:4]

其中,tif_file_path是原有计算平均值前遥感图像的保存路径,average_file_path是我们新生成的求取平均值后遥感影像的保存路径,也就是结果保存路径。

在这里,和我们前期的博客Python ArcPy批量拼接长时间序列栅格图像类似,需要首先在资源管理器中,将tif_file_path路径下的各文件以“名称”排序的方式进行排序;随后,利用arcpy.ListRasters()函数,获取路径下原有的全部.tif格式的图像文件,并截取第一个文件的部分文件名,从而获取其成像时间的具体年份。

接下来,遍历tif_file_path路径下全部.tif格式图像文件。其中,我们通过一个简单的判断语句if tif_file[0:4]==tif_file_year:,来确定某一年的遥感影像是否已经读取完毕——如果已经读取完毕,例如假如2001年成像的8幅遥感影像都已经遍历过了,那么就对这8景遥感影像加以逐像元的平均值求取,并开始对下一个年份(即2005年)成像的遥感影像继续加以计算;如果还没有读取完毕,例如假如2001年成像的8幅遥感影像目前仅遍历到了第5幅,那么就不求平均值,继续往下遍历,直到遍历完2001年成像的8幅遥感影像。

这里相信大家也看到了为什么我们要在前期先将文件夹中的文件按照“名称”排序——是为了保证同一年成像的所有遥感影像都排列在一起,遍历时只要遇到一个新的年份,程序就知道上一个年份的所有图像都已经遍历完毕了,就可以将上一个年份的所有栅格图像加以平均值求取。

在这里,逐像元的平均值求取其实也非常简单——我们对每一个像元分别执行以下操作:首先将该像元在当前年份里所有遥感影像的像素值相加,随后除以这一年份的遥感影像的数量,得到的就是该像元在这一年中像素值的平均值

最后,通过if tif_file==tif_file_name[len(tif_file_name)-1]:这个判断,来确认是否目前已经遍历到文件夹中的最后一个图像文件。如果是的话,就需要将当前成像年份的所有图像进行平均值的求取,并宣告代码完成运行。

在 IDLE (Python GUI) 中运行代码。代码运行完毕后,我们看一下结果文件夹。可以看到,其中的图像已经是按照成像时间,分别完成平均值求取后的结果了。

Python ArcPy实现批量计算多时相遥感影像的各项元平均值

在最后,还需要说明一点——用以上代码来求取长时间序列遥感影像的像元平均值,对于任意一个像元,只要该像元在任意一个时相的图像中是无效值(即为NoData),那么该像元在最终求出的平均值结果图中,像素值也将会是无效值NoData。针对这一问题的解决,我们将在下一篇博客中介绍。

来源:https://www.cnblogs.com/fkxxgis/p/17330411.html

标签:Python,ArcPy,遥感数据
0
投稿

猜你喜欢

  • SQL Server的FileStream和FileTable深入剖析

    2023-07-17 01:17:24
  • 301转向和网址规范化

    2007-09-26 14:00:00
  • python创建学生成绩管理系统

    2023-08-09 04:19:38
  • PHP自定义函数格式化json数据示例

    2023-07-17 07:17:45
  • 微信小程序实现页面跳转传递参数(实体,对象)

    2023-08-09 06:30:20
  • python开发之list操作实例分析

    2023-04-17 16:07:33
  • Django ORM 多表查询示例代码

    2021-07-25 05:22:02
  • python pandas.DataFrame.loc函数使用详解

    2023-10-04 07:01:58
  • 设计的技术含量

    2009-01-12 18:20:00
  • 百度“有啊”首页首次曝光 以绿色调为主

    2008-10-20 12:52:00
  • python基础教程之csv格式文件的写入与读取

    2021-05-24 09:20:12
  • Python 不设计 do-while 循环结构的理由

    2021-08-04 11:55:19
  • 教程javascript的function(函数)

    2007-09-30 13:38:00
  • 在SQL Server中编写通用数据访问方法

    2009-01-20 11:35:00
  • 在Go语言程序中使用gojson来解析JSON格式文件

    2023-07-02 14:19:52
  • python使用pandas按照行数分割表格

    2021-06-25 03:58:51
  • Python 进程操作之进程间通过队列共享数据,队列Queue简单示例

    2023-06-01 12:17:58
  • Python timeit模块原理及使用方法

    2023-09-22 08:19:18
  • 基于python图像处理API的使用示例

    2022-12-19 02:14:53
  • 略谈美国雅虎首页改版

    2009-07-26 18:34:00
  • asp之家 网络编程 m.aspxhome.com