R语言 空间网格 r语言空间插值

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

1. nc文件批量读取
2. 构造网格点，空间插画
3. 从字符串中提取日期，并转换
4. 日均转月均
5. 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]]<-d1
     newdate[[t]]<-mydate[[i]]
    idx[t]<-i
    t<-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)
}

完成。