问题引入
最近需要获得黄土高原2018逐月气温栅格数据,网上直接获取难度较大。所以就到气象数据网下载原始数据,通过ArcGIS插值来生成所需上数据。
在气象网上直接下载的数据是2018年全国166个站点逐日气温txt格式数据,如下图:
处理流程
首先,用R读入txt格式数据,逐站点按月求得逐月气温,接着转化成指定excel格式数据导入ArcGIS,然后通过插值生成逐月栅格图。
1.将txt数据导入R
由于txt原始数据字符以空格为间隔,但空格数不一致,使用read.table()方法不能奏效。但scan()函数可以逐行获取txt文件中的字符生成向量,通过循环将12个月气温数据整合。
#导入需要的包
library(dplyr)
library(reshape2)
#设置工作环境为气象数据所在文件夹
setwd('E:/1data/站点数据/气象网数据/气温')
#读入所以TXT格式数据,引号里的内容是匹配文件名的正则表达式
tem.name<-list.files(pattern = '*.TXT')
#scan()函数读入数据并存入向量x,由于数据都是13列,因此将x转化为13列的矩阵,将12个月的数据通过循环拼接到一起
tem.data<-data.frame()
for (f in tem.name) {
x=scan(f)
y=matrix(x,nrow = length(x)/13, byrow = TRUE)
tem.data=rbind(tem.data,y)
}
#可以选择输出成csv
write.csv(tem.data,'E:/1data/站点数据/气象网数据/气温/气温.csv')
初步得到如下数据框
> tem.data
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
1 50353 5144 12638 1739 2018 1 1 -141 -110 -184 9 9 9
2 50353 5144 12638 1739 2018 1 2 -120 -102 -160 9 9 9
3 50353 5144 12638 1739 2018 1 3 -177 -102 -295 9 9 9
4 50353 5144 12638 1739 2018 1 4 -215 -124 -261 9 9 9
5 50353 5144 12638 1739 2018 1 5 -211 -94 -283 9 9 9
2.求得各站点逐月降水
#由于我们只需要平均温度数据,因此删除无用列
tem.data<-tem.data[,-4]
tem.data<-tem.data[,-8:-12]
#设置列名
names(tem.data<-c('station','weidu','jingdu','year','month','day','tem')
#剔除缺失值(温度值大于30000的为缺失值)
tem1<-tem.data[tem.data$tem<30000,]
#group_by与summarise是dplyr包中的函数,结合使用可以轻松对数据进行分组求和。dcast函数将月份因子变成变量(比较抽象,但是尝试一下就能理解了)
by_day<-group_by(tem1,station,weidu,jingdu,year,month)
by_month<-summarise(by_day,tem=mean(tem))
x1<-dcast(by_month,station+weidu+jingdu+year~month)
x2<-paste('x',201801:201812,sep = '')
names(x1)<-c('station','weidu','jingdu','year',x2)
head(x1)
write.csv(x1,'E:/1data/站点数据/气象网数据/气温/2018.csv')
最终形式如下,导出为csv格式
> head(x1)
station weidu jingdu year x201801 x201802 x201803 x201804 x201805 x201806 x201807 x201808
1 50353 5144 12638 2018 -25.41935 -24.39286 -9.016129 4.033333 12.78065 17.27333 20.71935 17.94516
2 50434 5029 12141 2018 -30.40000 -27.43571 -12.351613 1.626667 10.02258 14.62667 18.00323 15.17097
3 50527 4915 11942 2018 -27.19355 -24.60714 -9.187097 4.520000 13.48387 18.63333 21.12258 19.09355
4 50557 4909 12514 2018 -23.95806 -22.81429 -7.577419 6.990000 14.45806 18.92667 22.50000 19.52903
5 50564 4925 12720 2018 -22.16129 -20.79286 -6.432258 6.073333 13.82581 17.77333 22.44194 19.61290
6 50658 4805 12547 2018 -22.89032 -19.21786 -6.351613 6.953333 14.89032 19.23000 22.76774 19.89032
ps:由于每个月有些站点经纬度居然会不一样(可能因为我气象网权限很低下载的数据质量很差),导致一个站点的数据被导出两遍,由于数量不多,所以就手动删除了。
3.将csv格式数据转化成excel并导入ArcGIS
将导出的csv格式数据在excel中另存为Excel 97-2003 工作表,命名为tem。
在ArcMAP中选择添加xy数据,这里有三个注意事项:
1.x字段填经度,y是纬度
2.输入坐标系必须是地理坐标系,不能是投影坐标系
3.检查经纬度是否是标准形式
4.插值生成栅格数据
得到全国气温站点图之后,通过spline插值,投影栅格,掩膜提取,最后得到黄土高原气温栅格图。
小结
由于完全没有专业基础,在自学得过程中遇到了很多问题。希望可以给和我一样从零开始的小伙伴一些参考,不足之处也请多指正。