Python ArcPy实现批量计算多时相遥感影像的各项元平均值
作者:疯狂学习GIS 发布时间:2022-06-04 14:45:51
在遥感应用中,我们经常需要对某一景遥感影像中的全部像元的像素值进行平均值求取——这一操作很好实现,基于ArcMap软件或者简单的Python代码就可以实现;但有时候,我们会需要结合同一地区、不同时相的多景遥感影像,求取每一个像元在全部时相中像素值的平均值——这一需求的实现较之前者就有些麻烦,本文对此加以介绍。
首先,我们来明确一下本文的具体需求。现有一个存储有大量.tif
格式遥感影像的文件夹,其中每一个遥感影像的文件名中都包含有该图像的成像时间,如下图所示。且其中除了.tif
格式的遥感影像文件外,还具有其它格式的文件。
我们希望,对于同一年成像的遥感影像进行逐像元平均值的求取。例如,上图中具有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) 中运行代码。代码运行完毕后,我们看一下结果文件夹。可以看到,其中的图像已经是按照成像时间,分别完成平均值求取后的结果了。
在最后,还需要说明一点——用以上代码来求取长时间序列遥感影像的像元平均值,对于任意一个像元,只要该像元在任意一个时相的图像中是无效值(即为NoData),那么该像元在最终求出的平均值结果图中,像素值也将会是无效值NoData。针对这一问题的解决,我们将在下一篇博客中介绍。
来源:https://www.cnblogs.com/fkxxgis/p/17330411.html


猜你喜欢
- 例如:JSON字符串:var str1 = '{ "name": "cxh", "
- 本文实例为大家分享了celery实现订单超时取消的具体代码,供大家参考,具体内容如下Celery官方文档中关于定时任务使用的说明项目目录结构
- 关于骨架屏介绍骨架屏的作用主要是在网络请求较慢时,提供基础占位,当数据加载完成,恢复数据展示。这样给用户一种很自然的过渡,不会造成页面长时间
- TensorFlow训练时,遇到内存不断增长,最终导致内存不足,进程被杀死。在这里我不准备对造成这一现象的所有原因进行探讨,只是记录一下我在
- 今天在继续学习Python时,打开Pycharm后,发现有一个项目下的项目文件名是红色的,如下图:刚开始我以为是我升级 Pycharm导致的
- CPU活动展示导入模块,创建画板,创建画笔进行绘画出cpu的数据,一定要用线程,负责会卡住哦实现代码import tkinterfrom t
- python-docx的简单使用'''设置表格所有单元格的四个边为0.5磅,黑色,实线可以使用返回值,也可以不使用&
- 个人开发的 flask 论坛进入尾声,还剩最后一个上传图片更换头像功能,搞了一整天,最后终于解决了所有问题,现在记录下解决方案。1. 上传文
- 1、安装依赖包yum -y install gcc-c++ ncurses-devel cmake make perl gcc autoco
- 正在看的ORACLE教程是:Oracle9iPL/SQL编程的经验小结。平时在PL/SQL中的编程中遇到一些问题,这里以问答的形式来进行把它
- 如下所示:<span style="font-size:18px;"># -*- coding:utf-8
- 动态变量名赋值在使用 tkinter 时需要动态生成变量,如动态生成 var1...var10 变量。使用 exec 动态赋值exec 在
- 本文实例讲述了Python设计模式之MVC模式。分享给大家供大家参考,具体如下:一.简单介绍mvc模式 the mo
- Pandas 按周、月、年、统计数据介绍将日期转为时间格式 并设置为索引import pandas as pddata=pd.read_ex
- 一、事件流 IE中是冒泡型事件,即从最特定的事件目标到最不特定的事件目标。 Netscape Navigator使用的是捕获型事件,这个跟I
- 具体如下: 1>如我们知道开始时间,要加减一个时间,得出一个结果时间,可以用以下代码 $time1='2008-10-1 12
- select CONVERT(varchar, getdate(), 120 ) 200
- 函数很简单, 主要是针对字符串和数字两种类型的传入数据分别进行了处理,具体用法:字符类型的strUsername = C
- 如果有一个字符串 eg: "sun,star,moon,clouds",想要在MS SQL中根据给定的分隔符',
- 一、环境介绍操作系统: win10 64位python版本: 3.8IDE: 采用vscode用到的相关安装包CSDN打包下载地址: htt