||
## 读取数据 # 数据格式:降水、最高气温、最低气温、风速。txt文件 vic_clim_forcing_read <- function(txtpath,staion_lat,station_lon,time_start,time_end, savepath){ # 变量的初始化 path <- txtpath lat <- staion_lat lon <- station_lon from <- time_start to <- time_end # 定位文件所在的位置 filename <- dir(path) filePath <- paste(path,filename,sep='/') # 批量读取数据 time <- seq(as.Date(from), as.Date(to),by = 'day') data <- list() for (i in 1:length(filePath)) { data[[i]] <- read.table(filePath[i],header = F,sep = '') time <- seq(as.Date(from), as.Date(to),by = 'day') long <- rep(lon[i],length(data[[i]])) latg <- rep(lat[i],length(data[[i]])) data[[i]] <- data.frame(latg,long,time,data[[i]]) } a <- data.frame() for (k in 1:length(time)){ # 时间 for (j in 1:length(data)){ # 数据 a[j,1:7] <- subset(data[[j]], data[[j]]$time == time[k]) filename <- paste(time[k], '.csv',sep = '') filepath <- savepath file <- paste(filepath, filename, sep = '/') write.csv(assign(paste('data_', time[k], sep = ''),a),file = file, sep = ',') message('正在保存', k) # assign(paste('data_', time[k], sep = ''), 分配变量名给数据 } } } ## 反距离插值 vic_clim_forcing_idw <- function(data_path,point_path,save_path){ library(sp) library(gstat) path<-data_path filename <- dir(path) filepath <- paste(path,filename,sep = '/') n <- length(filepath) dataCsv <- list() for (i in 1:n) { dataCsv[[i]]<-read.csv(filepath[i],header = T,sep = ',') } # 插值求值idw unknow<-read.csv(point_path,header = T,sep = ',')# 待插值点的经纬度 coordinates(unknow) <- ~ lat+lon e <- data.frame() # 获取每一天139个点的插值结果 for (j in 1:length(dataCsv)) { # 循环每一天对每一天进行插值 for (k in 4:7) { # 提取要插值的列 手动修改需要插值的数据所在的列 f <- dataCsv[[j]] # 提取第一天 coordinates(f) <- ~ latg + long data1 <- f@data # 准备每一天的数据 idwmodel = idw(data1[,k-1] ~1, f, unknow, maxdist = Inf, idp = 2) # 对每一列进行插值 e[1:139, k-3] <- idwmodel@data$var1.pred # 保存每一列的数据到数据框 filename1 <- paste(substr(filename[j],1,10), '.csv', sep = '') filepath <- save_path file <- paste(filepath, filename1, sep = '/') write.csv(e, file = file, sep = '') message('正在保存', j) } } } ##数据的写出 vic_clim_forcing_outTxt <- function(data_path,point_path,save_path){ path1<- data_path filename1 <- dir(path1) filepath1 <- paste(path1,filename1,sep = '/') n<-length(filepath1) dataCsv1 <- list() for (i in 1:n) { dataCsv1[[i]] <- read.csv(filepath1[i],header = T,sep = ',') } # 插值点的经纬度数据,用于输出文件的name unknow <- read.csv(point_path,header = T,sep = ',')# d <- data.frame() for (m in 1:length(unknow[,1])) { for (ii in 1:length(dataCsv1)) { d[ii,1:5] <- subset(dataCsv1[[ii]],dataCsv1[[ii]]$X == m) name<-paste('data',unknow[m,2],unknow[m,1],sep = '_') fil <- paste(save_path, name,sep = '/') write.table(round(d[,2:5],3),row.names = F,col.names = F,file = fil,sep = ' ') message('正在保存',ii) } } } ##函数的调用方式 # 由于R的内核是单线程单核心计算方式,所哟采用并行提高计算的速度 library(parallel) staion_lat <-c(37.1833, 37.3833 , 38.8, 37.2 , 37.9167, 38.2333 ) station_lon <- c(104.05, 101.6167, 101.0833, 102.8667, 102.6667 , 101.9667 ) system.time({ c1 <- detectCores(logical = F)-1 c2<-makeCluster(c1) clusterEvalQ(c2, library(xts)) clusterMap(c2, vic_clim_forcing_read, txtpath<-'C:/Users/zhufe/Desktop/dd/上游流域', MoreArgs=list(staion_lat,station_lon, time_start<-'2010/1/1', time_end<-'2012/12/31', savepath<-'C:/Users/zhufe/Desktop/dd/data')) stopCluster(c2) }) # system.time({ c1 <- detectCores(logical = F)-1 c2<-makeCluster(c1) clusterEvalQ(c2, library(xts)) clusterMap(c2,vic_clim_forcing_idw , data_path<-'C:/Users/zhufe/Desktop/dd/data', MoreArgs=list(point_path<-'C:/Users/zhufe/Desktop/point.csv', save_path<-'C:/Users/zhufe/Desktop/dd/data1')) stopCluster(c2) }) # system.time({ c1 <- detectCores(logical = F)-1 c2<-makeCluster(c1) clusterEvalQ(c2, library(xts)) clusterMap(c2,vic_clim_forcing_outTxt , data_path<-'C:/Users/zhufe/Desktop/dd/data1', MoreArgs=list(point_path<-'C:/Users/zhufe/Desktop/point.csv', save_path<-'D:/data4')) stopCluster(c2) }) 插值的版本是最低的!!! 可以完善修改!!! 2019-3-2
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 08:14
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社