数据基础

你有一个基准经纬度的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)

}

文章链接

评论可见,请评论后查看内容,谢谢!!!评论后请刷新页面。