NCL处理nc就是会更方便,但是避免不了一些用不了NCL的情况
写了一个简单用matlab提取nc文件站点数据的脚本,通过判断经纬度与格点经纬度距离,选取最近格点值。
后面也附上了提取wrfout中站点数据的NCL脚本,如果是再分析资料的话,wrf系列函数不知道是否可以,可能还是需要通过距离判断。
matlab code:
clc; %清屏
clear all; %清空
%%
%nc批量读取的数据
datadir1='C:\Users\sunxy\Desktop\'; %指定批量数据所在的文件夹
filelist1=dir([datadir1,'200901.nc']); %指定批量数据的类型
%%
%读取数据
lat_station=32.5;lon_station=90.8;
filename1=[datadir1,filelist1.name];
ncid1=netcdf.open(filename1,'NC_NOWRITE'); %打开nc文件
lat=ncread(filename1,'latitude'); %读入变量lat
lon=ncread(filename1,'longitude'); %读入变量lon
[lat_2D lon_2D]=meshgrid(lat,lon);
distance=sqrt((lat_2D-lat_station).^2+(lon_2D-lon_station).^2); %算出站点与每个格点距离
[ilon ilat]=find(distance==min(min(distance))); %寻找距离最近格点
lat_back=lat(ilat); %输出经纬度检查是否为希望读取站点
lon_back=lon(ilon,:);
blh= ncread(filename1,'blh');
blh_station=blh(ilon,ilat,:);
netcdf.close(ncid1); %关闭nc文件
%
% % ncdisp('C:\Users\sunxy\Desktop\200901.nc','/','full')
NCL code:
row = numAsciiRow("/datadir1/station.txt")
col = numAsciiCol("/datadir1/station.txt")
f_h = asciiread("/datadir1/station.txt",(/row,col/),"float")
station=f_h(:,0)
lon=f_h(:,1)
lat=f_h(:,2)
station_num=dimsizes(station)
;read wrfout u v data
DATADir=(/"/datadir1/zhuyan/wrfout/2020/expt1/18/"/)
FILE= systemfunc (" ls -1 " + DATADir + "wrfout_d01_2019*")
filenum=dimsizes(FILE)
a = addfiles(FILE + ".nc", "r")
u = wrf_user_getvar(a,"U10",-1)
output_u=new((/station_num,filenum+5/), float)
output_u(:,0)=station
output_u(:,1)=lon
output_u(:,2)=lat
do istation=0,station_num-1
loc = wrf_user_ll_to_ij(a,lon(istation), lat(istation), True)-1
output_u(istation,3)=loc(0)
output_u(istation,4)=loc(1)
output_u(istation,5::)=u(:,loc(1),loc(0))
end do
opt =True
M=filenum+5
fWidth = 7 ; specify the format width
fDec= 2 ; specify the number to the right of decimal point
fmtf= M + "f" + fWidth + "." + fDec ; same as "7f7.2"
opt@fout ="output_u_202001.txt"
write_matrix (output_u, fmtf, opt)