栅格STI计算-基于R的STI-package

" 1. 程序目的 (1) 计算1960-2014年中国矩形区域的STI 2. 2024年11月29日 Version 1 3. 数据 根路径：'E:/05_Study/STI_Cal/01_Data/' 3.1 输入数据 01_Temperature/china_matrix_mulmodelmean_month_3d_History_tas_1960-2014.tif 3.2 输出数据 02_STI_Result路径下的文件 " # ----------------------相关包的导入------------------------------- library(sp) library(raster) library(STI) library(tiff) # ---------------------函数预定义---------------------------------- # 循环读取栅格数据计算STI sti_raster <- function(tasdata,scale_aim){ tasdim = dim(tasdata) sti_data <- array(-99,dim=c(tasdim[1],tasdim[2],tasdim[3])) for(roww in c(1:(tasdim[1]))){ for(coll in c(1:(tasdim[2]))){ tas_temp <- tasdata[roww,coll,] sti_data[roww,coll,] = sti(tas_temp,scale_aim) #print(coll) } print(roww) } return(sti_data) } # 栅格数据写出 raster_write <- function(aim_data,tif_infor,file_path,out_namepart01,styr,edyr){ for (year in c(styr:edyr)){ for (month in c(1:12)){ # 特定层数据提取 idx = month + (year-styr)*12 aim_data_temp = aim_data[,,idx] # 特定层数据写出 # 文件名 outnametemp = paste(out_namepart01,as.character(year*100+month),sep='_') outname = paste(outnametemp,'tif',sep='.') # 文件完整路径 outfileaim = paste(file_path,outname,sep='/') # 创建raster raster_temp <- raster(extent(tif_infor@extent@xmin,tif_infor@extent@xmax,tif_infor@extent@ymin,tif_infor@extent@ymax), res = 0.5, crs = 4326 ) raster_temp <- setValues(raster_temp,aim_data_temp) # raster的写出 raster_bb <- writeRaster(raster_temp,outfileaim,overwrite=TRUE) } print(year) } } # ---------------------主程序-------------------------------------- # 路径处理和基本变量定义 rootdir = 'E:/05_Study/STI_Cal/01_Data' tiffile <- paste(rootdir,'china_matrix_mulmodelmean_History_Nat_tas_196001.tif',sep='/') # 辅助获取tif数据的基本信息 aim_scale =3 # 计算3个月尺度的STI styear = 1960 # 数据开始年 edyear = 2014 # 数据结束年 # 输入数据 aim_tas = paste(rootdir,'01_Temperature/china_matrix_mulmodelmean_month_3d_History_tas_1960-2014.tif',sep='/') # 输出路径 outpath = paste(rootdir,'02_STI_Result',sep='/') # 输出文件名第一部分 outnamepart01 = paste('china_matrix_mulmodelmean_month_3d_History_sti03',sep='_') # raster读取tif数据,获取tif数据基本信息 tif_temp <- raster(tiffile) # 读取.tif数据为三维数组，此处aim_tas只有一个文件，可以直接读取 tasdata_3d <- readTIFF(aim_tas) # 计算STI sti_3d = sti_raster(tasdata_3d,aim_scale) # STI计算结果写出 raster_write(sti_3d,tif_temp,outpath,outnamepart01,styear,edyear) # ------------------------------------------------------------------------- print('Finished.')