問題引入
最近需要獲得黃土高原2018逐月氣溫柵格數(shù)據(jù)朵栖,網(wǎng)上直接獲取難度較大。所以就到氣象數(shù)據(jù)網(wǎng)下載原始數(shù)據(jù)柴梆,通過ArcGIS插值來生成所需上數(shù)據(jù)陨溅。
在氣象網(wǎng)上直接下載的數(shù)據(jù)是2018年全國166個站點逐日氣溫txt格式數(shù)據(jù),如下圖:
處理流程
首先绍在,用R讀入txt格式數(shù)據(jù)门扇,逐站點按月求得逐月氣溫雹有,接著轉(zhuǎn)化成指定excel格式數(shù)據(jù)導入ArcGIS,然后通過插值生成逐月柵格圖臼寄。
1.將txt數(shù)據(jù)導入R
由于txt原始數(shù)據(jù)字符以空格為間隔霸奕,但空格數(shù)不一致,使用read.table()方法不能奏效吉拳。但scan()函數(shù)可以逐行獲取txt文件中的字符生成向量质帅,通過循環(huán)將12個月氣溫數(shù)據(jù)整合。
#導入需要的包
library(dplyr)
library(reshape2)
#設置工作環(huán)境為氣象數(shù)據(jù)所在文件夾
setwd('E:/1data/站點數(shù)據(jù)/氣象網(wǎng)數(shù)據(jù)/氣溫')
#讀入所以TXT格式數(shù)據(jù)留攒,引號里的內(nèi)容是匹配文件名的正則表達式
tem.name<-list.files(pattern = '*.TXT')
#scan()函數(shù)讀入數(shù)據(jù)并存入向量x煤惩,由于數(shù)據(jù)都是13列,因此將x轉(zhuǎn)化為13列的矩陣炼邀,將12個月的數(shù)據(jù)通過循環(huán)拼接到一起
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/站點數(shù)據(jù)/氣象網(wǎng)數(shù)據(jù)/氣溫/氣溫.csv')
初步得到如下數(shù)據(jù)框
> 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.求得各站點逐月降水
#由于我們只需要平均溫度數(shù)據(jù)魄揉,因此刪除無用列
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包中的函數(shù),結合使用可以輕松對數(shù)據(jù)進行分組求和拭宁。dcast函數(shù)將月份因子變成變量(比較抽象洛退,但是嘗試一下就能理解了)
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/站點數(shù)據(jù)/氣象網(wǎng)數(shù)據(jù)/氣溫/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:由于每個月有些站點經(jīng)緯度居然會不一樣(可能因為我氣象網(wǎng)權限很低下載的數(shù)據(jù)質(zhì)量很差)杰标,導致一個站點的數(shù)據(jù)被導出兩遍兵怯,由于數(shù)量不多,所以就手動刪除了腔剂。
3.將csv格式數(shù)據(jù)轉(zhuǎn)化成excel并導入ArcGIS
將導出的csv格式數(shù)據(jù)在excel中另存為Excel 97-2003 工作表媒区,命名為tem。
在ArcMAP中選擇添加xy數(shù)據(jù)桶蝎,這里有三個注意事項:
1.x字段填經(jīng)度驻仅,y是緯度
2.輸入坐標系必須是地理坐標系,不能是投影坐標系
3.檢查經(jīng)緯度是否是標準形式
4.插值生成柵格數(shù)據(jù)
得到全國氣溫站點圖之后登渣,通過spline插值噪服,投影柵格,掩膜提取胜茧,最后得到黃土高原氣溫柵格圖粘优。
小結
由于完全沒有專業(yè)基礎,在自學得過程中遇到了很多問題呻顽。希望可以給和我一樣從零開始的小伙伴一些參考雹顺,不足之處也請多指正。