R语言提取站点的nc文件时间序列数据

数据基础

你有一个基准经纬度的nc文件
你有一个多站点的经纬度文件
你有多年包含月天小时的nc文件

# 安装和加载所需的包
#install.packages("ncdf4")
library(ncdf4)
library(readxl)
library(openxlsx)
library(raster)
library(geosphere)
rm(list=ls()) 
setwd("E:/Program1/input/surfout-beijing-202007")

# 读取包含经纬度的NetCDF文件
#nc_open 将数据读入我称为 nc_data 的数据结构中
nc_coords <- nc_open('../wrfout_d03_2020-08-02_00_00_00.nc')
crs(nc_coords)
names(nc_coords $var)
# 读取经纬度变量lat"通常用来表示纬度（Latitude）
lat <- ncvar_get(nc_coords, "XLAT")
lon <- ncvar_get(nc_coords, "XLONG")
dim(lat)
# 关闭文件
nc_close(nc_coords)
#读取站点信息
Stations=read_excel("../Beijing_staions.xlsx",sheet = 1, col_names = c("beijing","stations","city","lat","lon"), col_types = NULL, na = "", skip = 0)
Stations <- as.data.frame(Stations)
indexs=c("ALBEDO","T2","Q2","U10","V10","PSFC")

# 提取WRDF对应站点模拟结果
for (station in 1:18){
  city=Stations[station,2]
  station_lat=Stations[station,]$lat
  station_lon=Stations[station,]$lon
  
  # 计算站点与每个网格点的距离
  #但是经纬度不等同于米的计算所以不可以
  #计算栈距离最近的网格
  distance=data.frame(matrix( nrow = 150, ncol = 150))
  for(i in 1:150){
    for(j in 1:150){
      #distance[i,j] <- distGeo(c(lon[i,j],lat[i,j]), c(station_lon,station_lat))
      distance[i,j] <- distm(c(lon[i,j],lat[i,j]), c(station_lon,station_lat))
    }
  }
  #distance是米的距离
  # 找到最小距离的索引按列平铺
  #min_index <- which.min(distances)
  #返回最小值的行列号
  #arr.ind = TRUE 则告诉函数返回的索引是数组形式的，而不是向量形式的
  value_hl=which(distance==min(distance),arr.ind=T)#71  30
  # 打开每日的nc文件，提取数据
  station_city <- data.frame(matrix(nrow = 31 * 24, ncol = length(indexs) + 2))
  colnames(station_city) <- c("day", "t","ALBEDO","T2","Q2","U10","V10","PSFC")

  #将nc读取合并
    for (day in 1:31) {
    # 假设站点数据存储在另一个文件中
    julyday <- sprintf("%02d", day)
    nc=paste("surfout_d03_2020-07-", julyday, "_00_00_00", sep = "")
    ncfile <- nc_open(nc) 
    for(dex in 1:6){
      index=indexs[dex]
      index_data=ncvar_get(ncfile,index)
      for(t in 1:24){
        index_time=index_data[,,t]
        index_value=index_time[70,30]
        station_city[(day - 1) * 24 + t, "t"] <- t
        station_city[(day - 1) * 24 + t, "day"] <- day
        station_city[(day - 1) * 24 + t, index] <- index_value
       }
    }
  }
  nc_close(ncfile)
  # 将结果写入Excel文件
 output_folder <- "../nc/"
 station_city_file <- paste(output_folder,"station", city, "_data.xlsx", sep = "")
  write.xlsx(station_city,station_city_file , rowNames= FALSE)
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76