R语言空间插值/R语言离散数据网格化/R语言空间点插值/R语言nc日均转月均、日期转换

在我们日常使用数据时，我们获取的数据并不是连续的网格，而是网格状的离散数据。
这时，我们就需要利用拥有的离散网格数据，自定义将其插值。
通常插值我们需要shp文件给予边界,不过，只要知道了所需区域的经纬度范围，我们可以自行构建网格将点数据插值在自己定义的经纬度范围中。
本文包含：

nc文件批量读取
构造网格点，空间插画
从字符串中提取日期，并转换
日均转月均
nc格式导出

数据说明与读取

本次数据为OCO系列卫星数据的CO2数据，空间分辨率：约为1.75km网格状的离散数据，存储方式为nc。

library('ncdf4')
library('lubridate')
library('sp')
library('raster')
library('magrittr')
library('dplyr')
library('gstat')
library('stringr')
ncfiles<-list.files(pattern = ".nc")
data<-nc_open(ncfiles[1])
#sink('ncinfo.txt',split = TRUE)
print(data)#查看nc文件基本信息,并输出保存
#根据数据信息，开始提取对应变量
c_data<-list()
for (i in 1:length(ncfiles)) {data<-nc_open(ncfiles[i])lat<-ncvar_get(data,'latitude')lon<-ncvar_get(data,'longitude')cpf<-ncvar_get(data,'xco2')c_data[[i]]<-data.frame(lon,lat,cpf)
}

此读取的数据长这样：

由于数据是离散网格状，每个nc文件的数据长度不同，lat、lon也不同，这样，我们若要进行批量的计算话，很有必要将其处理成等维度的大小。

数据筛选

该数据为全球范围，我们需要的部分为中国范围，可以根据中国的经纬度范围提取对应的数据，提高运算速率。

#中国经纬度范围：3.86至53.55°N，73.66至135.05°E
#由于数据本身为全球范围离散数据，插值时可直接构造等经纬度格点
#数据分辨率约1.75KM,约为0.0175°，
#由于分辨率过高插值耗时较长，此处插值为0.1°格点，分辨率为10Km
#根据自身需要，修改c_lat与c_lon中seq中的by即可
c_lat=seq(3.86,53.55,0.1)
c_lon=seq(73.66,135.05,0.1)```#从原始数据中筛选出中国地区数据
newdata<-list()
newdate<-list()
idx<-vector(mode="numeric",length=0)
t<-1
for (i in 1:length(c_data)) {d<-c_data[[i]]d1<-subset(d,lon>=73.66&lon<=135.05&lat>=3.86&lat<=53.55,select=c(lon,lat,cpf))if(dim(d1)[1]>0){newdata[[t]]<-d1newdate[[t]]<-mydate[[i]]idx[t]<-it<-t+1}
}

网格构建与空间插值

在插值前，我们需要构建自己需要的网格。
由于我们并不涉及栅格投影转换问题，此处CRS等信息无需考虑。
由于插值后的数据太大，我们也需要把数据重新改成数据框形式。

#空间插值
#构造插值网格点
c_grid <- expand.grid(x=c_lon,y=c_lat)
coordinates(c_grid) <- ~x+y
gridded(c_grid) <- TRUE
#插值，IDW插值inverse distance weighted interpolation
idw_data<-list()
for (i in 1:length(newdata)) {d<-newdata[[i]]c_sp<-SpatialPointsDataFrame(coords = cbind(x =d$lon, y =d$lat),data=d)#将数据转为空间栅格类型idw<-idw(formula =cpf~1,locations = c_sp,newdata=c_grid)idw_output <- as.data.frame(idw)names(idw_output)[1:3]<-c("long","lat","cpf")idw_data[[i]]<-idw_output[,c("long","lat","cpf")]print(i)
}

日期转换

虽然nc文件里有date，但从文件里读取速度较慢，我们根据文件名提取日期，再进行转换。

mydate<-list()
for (i in 1:length(ncfiles)) {chdate<-str_extract(ncfiles[i],"[0-9]{6}")mydate[[i]]<-as.Date(chdate,'%y%m%d')
}

日均转月均

主要需要注意数据转换。
使用lappy和do.call也可以实现，但速度慢，我将其转为matrix，再用apply函数计算。

newdate<-as.Date(unlist(newdate),origin = '1970/1/1')
mon<-month(newdate)
ave_mon<-list()
for (ii in 1:12) {idx<-which(mon==ii)cpf_mon<-array(0,dim = c(length(idx),305158,3))for (j in 1:length(idx)) {cpf_mon[j,,]<-as.matrix(idw_data[[idx[j]]])}ave_mon[[ii]]<-apply(cpf_mon, c(2,3), mean)
}

nc文件导出

#导出为nc文件
coord<-c_grid@coords
y<-ncdim_def('y',units = 'degrees',vals =coord[,2] )
x<-ncdim_def('x',units = 'degrees',vals = coord[,1])
for (i in 1:12) {t<-ncdim_def('time',units = 'months',vals =i)ave_cof<-ncvar_def( name = 'xco2',units = 'ppm', dim = list(x),prec = 'double' )lat<-ncvar_def( name = 'lat',units = 'degrees', dim = list(x),prec = 'double' )lon<-ncvar_def( name = 'lon',units = 'degrees', dim = list(x),prec = 'double' )ncnew <- nc_create(filename =paste('ave_cof/avemon_',as.character(i),'.nc',seq='',collapse = ''),list(ave_cof,lat,lon))ncvar_put(ncnew, varid = ave_cof, vals =ave_mon[[i]][,3] )ncvar_put(ncnew, varid = lat, vals =ave_mon[[i]][,2] )ncvar_put(ncnew, varid = lon, vals =ave_mon[[i]][,1] )nc_close(ncnew)
}

完成。