402com永利平台|402com永利1站|55.402com永利网址

您的位置:402com永利平台 > 科学研究 > 批量提取逐日tif文件的像元值,边学边写博客

批量提取逐日tif文件的像元值,边学边写博客

2019-08-26 08:48

批量提取逐日tif文件的像元值——基于arcgis and R

代码语言居然没有R##arcgis准备提取点的shp文件#研究区域的shp文件#使用arcgis栅格重采样工具确定DEM的分辨率#extractbymask工具确定研究区域的栅格范围#rastertopoint工具实现栅格转shp#栅格、shp统一的地理坐标系,所有的栅格影像要统一投影坐标系#在shp文件的属性表中,添加经度和纬度字段,获取每一个点的经纬度。##R中实现批量提取数据#定义函数extr_tif<-function(Tifpath,SHPpath,savepath){#Tifpathtif文件所在的文件夹路径#SHPpathsho文件所在的文件夹路径#savepath数据保存到本地的路径librarylibrarylibrary#定义数据框和列表data1<-list()b<-data.frame()#变量的初始化pathshp<-SHPpathtifpath<-Tifpathfilepath<-savepathsetwd##读取shp文件ogrInfop<-readOGR(dsn=pathshp)#shp文件#读取tif文件file<-list.files(tifpath,pattern=".tif$",full.names=FALSE)#modis转tif的数据,所以要确定每一期影像的日期#modis的日期转化成年—月—日a<-extractDate(file,asDate=TRUE)date<-a$inputLayerDatesa<-p@data#提取tif数据for(jin1:length{data1[[file[j]]]<-extract(raster,p)b<-as.data.frame(data1[[j]])e<-data.frame(a$lat,a$lon,a$grid_code,b)#保存数据filename<-paste(date[j],'.csv',sep='')file1<-paste(filepath,filename,sep='/')write.csv(assign(paste('data_',date[j],sep='_'),e),file=file1,row.names=F)}}#并行计算librarysystem.time({cl<-detectCores(logical=F)-1c2<-makeClusterclusterEvalQ(c2,libraryclusterMap(c2,extr_tif,Tifpath<-'',#添加路径,MoreArgs=list(SHPpath<-'',#添加路径savepath<-''))#添加路径stopCluster2019-2-22

##读取数据#数据格式:降水、最高气温、最低气温、风速。txt文件vic_clim_forcing_read<-function(txtpath,staion_lat,station_lon,time_start,time_end,savepath){#变量的初始化path<-txtpathlat<-staion_latlon<-station_lonfrom<-time_startto<-time_end#定位文件所在的位置filename<-dirfilePath<-paste(path,filename,sep='/')#批量读取数据time<-seq(as.Date,as.Date,by='day')data<-list()for(iin1:length){data[[i]]<-read.table(filePath[i],header=F,sep='')time<-seq(as.Date,as.Date,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(kin1:length{#时间for(jin1:length{#402com永利1站,数据a[j,1:7]<-subset(data[[j]],data[[j]]$time==time[k])filename<-paste(time[k],'.csv',sep='')filepath<-savepathfile<-paste(filepath,filename,sep='/')write.csv(assign(paste('data_',time[k],sep=''),a),file=file,sep=',')message#assign(paste('data_',time[k],sep=''),分配变量名给数据}}}##反距离插值vic_clim_forcing_idw<-function(data_path,point_path,save_path){librarylibrarypath<-data_pathfilename<-dirfilepath<-paste(path,filename,sep='/')n<-lengthdataCsv<-list()for{dataCsv[[i]]<-read.csv(filepath[i],header=T,sep=',')}#插值求值idwunknow<-read.csv(point_path,header=T,sep=',')#待插值点的经纬度coordinates<-~lat lone<-data.frame()#获取每一天139个点的插值结果for(jin1:length{#循环每一天对每一天进行插值for{#提取要插值的列手动修改需要插值的数据所在的列f<-dataCsv[[j]]#提取第一天coordinates<-~latg longdata1<-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_pathfile<-paste(filepath,filename1,sep='/')write.csv(e,file=file,sep='')message}}}##数据的写出vic_clim_forcing_outTxt<-function(data_path,point_path,save_path){path1<-data_pathfilename1<-dirfilepath1<-paste(path1,filename1,sep='/')n<-length(filepath1)dataCsv1<-list()for{dataCsv1[[i]]<-read.csv(filepath1[i],header=T,sep=',')}#插值点的经纬度数据,用于输出文件的nameunknow<-read.csv(point_path,header=T,sep=',')#d<-data.frame()for(min1:length(unknow[,1])){for(iiin1:length){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的内核是单线程单核心计算方式,所哟采用并行提高计算的速度librarystaion_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)-1c2<-makeClusterclusterEvalQ(c2,libraryclusterMap(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#system.time({c1<-detectCores(logical=F)-1c2<-makeClusterclusterEvalQ(c2,libraryclusterMap(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#system.time({c1<-detectCores(logical=F)-1c2<-makeClusterclusterEvalQ(c2,libraryclusterMap(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插值的版本是最低的!!!可以完善修改!!!2019-3-2

有几个excel表,但是每个表中的数据都不一样,所以读进R里面进行简单的处理。过程很简单,基本是数据筛选,重命名,添加列还有合并。还包含了读取表和存储表的注意事项。

有几个excel表,但是每个表中的数据都不一样,所以读进R里面进行简单的处理。过程很简单,基本是数据筛选,重命名,添加列还有合并。还包含了读取表和存储表的注意事项。

#百度
setwd("e:/baidu")#设置工作目录
a<-read.table("20151125.csv",sep=",",fill=T,skip=7,header = T)#读取文件,跳过空白
a[,3]<-"baidu"#修改第三列的内容
b<-read.table("20151125.csv",sep=",",fill=T,skip=7,header = T)
b[,3]<-"baixingwang"
c<-rbind#合并两表
baidu<-data.frame(c$日期,c$账户,c$推广计划,c$展现,c$点击,c$消费)#挑选有用变量

#百度
setwd("e:/baidu")#设置工作目录
a<-read.table("20151125.csv",sep=",",fill=T,skip=7,header = T)#读取文件,跳过空白
a[,3]<-"baidu"#修改第三列的内容
b<-read.table("20151125(1).csv",sep=",",fill=T,skip=7,header = T)
b[,3]<-"baixingwang"
c<-rbind(a,b)#合并两表
baidu<-data.frame(c$日期,c$账户,c$推广计划,c$展现,c$点击,c$消费)#挑选有用变量

colnames=c("日期","账户","推广计划","展现","点击","消费")#重新命名
#搜狗
setwd("e:/sogou")
d<-read.table("[搜狗推广服务]定制报告-搜狗推广报告-2015-11-25_全设备.csv",sep=",",fill=T,header = T)
d<-d[-1,]
sogou<-data.frame(d$日期,d$账户,d$推广计划,d$展示数,d$点击数,d$消耗)
sogou[,2]<-"sogou"
colnames=c("日期","账户","推广计划","展现","点击","消费")
#360点睛
setwd
e<-read.table("2015-11-24推广计划数据报告.csv",sep=",",fill=T,header = T)
日期<-c(1:length
账户<-c(1:length
e<-data.frame
e$日期<-"2015/11/24"
e$账户<-"360点睛"
dianjing<-data.frame(e$日期,e$账户,e$推广计划,e$展示次数,e$点击次数,e$总费用)
dianjing[,2]<-"360"
colnames=c("日期","账户","推广计划","展现","点击","消费")
#神马=============
setwd("e:/shenma")
library
f<-read.xlsx("计划推广报告-2015-11-25~2015-11-25#2015-11-26 15-46-38.xlsx",1,encoding='UTF-8')
f<-f[,1:6]
shenma<-data.frame
shenma[,2]<-"shenma"
colnames=c("日期","账户","推广计划","展现","点击","消费")
#合并============
all<-rbind(baidu,sogou,dianjing,shenma)#纵向合并
setwd("e:/output")
write.csv(all,"四渠道",row.names = FALSE)
#与四渠道合并
setwd
#quote="/""认为只有双引号才分隔,这样Xi’an可以在一起
g<-read.table("Analytics 全站数据 二手房整站页 20151030-20151129.csv",sep=",",fill=T,header = T,skip=6,quote = """)
m<-merge(g,all,by="推广计划",all=T)
setwd("e:/output")
write.csv(m,"四渠道匹配.csv")
#读取GA数据=============================================
#全站pv导航============================
setwd
b<-read.table("Analytics 全站数据 PV 导航 20151116-20151122.csv",fill=T,sep=",",header = F,skip=7)
b<-b[1:(which[1]-1),]#根据第一列的长度来截取数据。因为这个表笔记特殊,上面是四列,下面两列,两列的数据用不到。
colnames=c("来源","会话","新会话百分比","新用户","跳出率","每次会话浏览页数","平均会话时长","用户数","浏览量")

colnames(baidu)=c("日期","账户","推广计划","展现","点击","消费")#重新命名
#搜狗
setwd("e:/sogou")
d<-read.table("[搜狗推广服务]定制报告-搜狗推广报告-2015-11-25_全设备.csv",sep=",",fill=T,header = T)
d<-d[-1,]
sogou<-data.frame(d$日期,d$账户,d$推广计划,d$展示数,d$点击数,d$消耗)
sogou[,2]<-"sogou"
colnames(sogou)=c("日期","账户","推广计划","展现","点击","消费")
#360点睛
setwd("e:/360")
e<-read.table("2015-11-24推广计划数据报告.csv",sep=",",fill=T,header = T)
日期<-c(1:length(e$推广计划))
账户<-c(1:length(e$推广计划))
e<-data.frame(e,日期,账户)
e$日期<-"2015/11/24"
e$账户<-"360点睛"
dianjing<-data.frame(e$日期,e$账户,e$推广计划,e$展示次数,e$点击次数,e$总费用)
dianjing[,2]<-"360"
colnames(dianjing)=c("日期","账户","推广计划","展现","点击","消费")
#神马=============
setwd("e:/shenma")
library(xlsx)
f<-read.xlsx("计划推广报告-2015-11-25~2015-11-25#2015-11-26 15-46-38.xlsx",1,encoding='UTF-8')
f<-f[,1:6]
shenma<-data.frame(f)
shenma[,2]<-"shenma"
colnames(shenma)=c("日期","账户","推广计划","展现","点击","消费")
#合并============
all<-rbind(baidu,sogou,dianjing,shenma)#纵向合并
setwd("e:/output")
write.csv(all,"四渠道",row.names = FALSE)
#与四渠道合并
setwd("e:/ga")
#quote="/""认为只有双引号才分隔,这样Xi’an可以在一起
g<-read.table("Analytics 全站数据 二手房整站页 (计划) 20151030-20151129.csv",sep=",",fill=T,header = T,skip=6,quote = """)
m<-merge(g,all,by="推广计划",all=T)
setwd("e:/output")
write.csv(m,"四渠道匹配.csv")
#读取GA数据=============================================
#全站pv导航============================
setwd("e:/ga")
b<-read.table("Analytics 全站数据 PV 导航 20151116-20151122.csv",fill=T,sep=",",header = F,skip=7)
b<-b[1:(which(b$V1=="")[1]-1),]#根据第一列的长度来截取数据。因为这个表笔记特殊,上面是四列,下面两列,两列的数据用不到。
colnames(b)=c("来源","会话","新会话百分比","新用户","跳出率","每次会话浏览页数","平均会话时长","用户数","浏览量")

本文由402com永利平台发布于科学研究,转载请注明出处:批量提取逐日tif文件的像元值,边学边写博客

关键词: