数据基础
你有一个基准经纬度的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)
}
文章链接
发表评论