今天又破案了,应该用wrfout数据与观测数据进行对比,具体是:
气温:2米处的气温
相对湿度:2米的相对湿度
wrf_user_getvar (ucar.edu)
风:10米处的风
气压:海平面气压
NCL: Writing CSV (comma-separated values) files (ucar.edu)
write_csv_5.ncl
这个脚本非常实用,我贴在这里方便学习。
技术点包括:
①读多个文件
②根据多个站点经纬度找格点位置
③将结果写成csv文件输出
;----------------------------------------------------------------------
; write_csv_5.ncl
;
; Concepts illustrated:
; - Writing a CSV file with a header using write_table
; - Appending data of mixed types to a CSV file inside a loop
; - Writing select WRF-ARW data to a CSV file
;----------------------------------------------------------------------
; This example calculates temperature at 2m from a WRF-ARW output
; file, and writes a subset of the data based on an array of
; lat / lon values.
;----------------------------------------------------------------------
; These files are loaded by default in NCL V6.4.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
dir = "./"
files = systemfunc (" ls -1 " + dir + "wrfout_d01_2008-09* ")
a = addfiles(files+".nc","r")
times = wrf_user_list_times(a) ; "2008-09-29_18:30:00", etc
tk2 = wrf_user_getvar(a,"T2",-1) ; T2 in Kelvin
times = wrf_user_list_times(a)
ntimes = dimsizes(times)
print("ntimes = " + ntimes)
;---Calculate i,j locations of data closest to set of lat/lon points
lats = (/ 30, 40, 50/)
lons = (/130,135,140/)
nlatlon = dimsizes(lats)
loc = wrf_user_ll_to_xy(a, lons, lats, True) ; 2 x nlatnlon
; loc = wrf_user_ll_to_ij(a, lons, lats, True) ; 2 x nlatnlon
; loc = loc - 1 ; wrf_user_ll_to_ij returns 1-based indexes
;---Set up CSV file and header information for the file
csv_filename = "wrf_2m_temperature.csv"
system("rm -f " + csv_filename) ; Remove file in case it exists.
fields = (/"TIME", "LAT", "LON", "TEMPERATURE (degC)"/)
;---Create a header line for CSV file
dq = str_get_dq()
fields = dq + fields + dq ; Pre/append quotes to field names
header = [/str_join(fields,",")/] ; Header is field names separated
; by commas.
;
; Format to use for writing each variable to CSV file.
; If you don't want spaces in CSV file, use the following
; format string:
; format = "%s,%g,%g,%g"
;
format = "%s,%6.2f,%7.2f,%6.2f"
;
; Loop through each time step and desired list of lat/lon values,
; and write a single line of data to CSV file.
;
write_table(csv_filename, "w", header, "%s") ; Write header to CSV file.
do it = 0,ntimes-1
do nl = 0,nlatlon-1
nln = loc(0,nl)
nlt = loc(1,nl)
lat1 = a[0]->XLAT(0,nlt,nln) ; nearest grid point
lon1 = a[0]->XLONG(0,nlt,nln)
alist = [/times(it),lat1,lon1,tk2(it,nlt,nln)/] ; Store data to be written in a list.
write_table(csv_filename, "a", alist, format) ; Write list to CSV file.
end do
end do
end
我自己的脚本
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ; These two libraries are automatically
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ; loaded from NCL V6.4.0 onward.
;
begin
dir = "./201912/"
files = systemfunc (" ls -1 " + dir + "wrfout_d03* ") + ".nc"
a = addfiles(files,"r")
times = wrf_user_list_times(a)
ntimes = dimsizes(times)
print("ntimes = " + ntimes)
lats = (/39.13/)
lons = (/115.8/)
nlatlon = dimsizes(lats)
loc = wrf_user_ll_to_xy(a, lons, lats, True)
idlon = loc(0)
idlat = loc(1)
loc2 = (/(/idlon-1,idlon,idlon+1,idlon-1,idlon,idlon+1,idlon-1,idlon,idlon+1/),(/idlat-1,idlat,idlat+1,idlat,idlat+1,idlat-1,idlat+1,idlat-1,idlat/)/)
nlatlon = 9
print(loc2)
;idlon = 145
;idlat = 109
time = -1
slp = wrf_user_getvar(a,"slp",time) ; slp
tc2 = wrf_user_getvar(a,"T2",time) ; T2 in Kelvin
u10 = wrf_user_getvar(a,"U10",time) ; u at 10 m
v10 = wrf_user_getvar(a,"V10",time) ; v at 10 m
rh2 = wrf_user_getvar(a,"rh2",time) ; rh at 2 m
lat2d = wrf_user_getvar(a, "lat", time)
lon2d = wrf_user_getvar(a, "lon", time)
;---Set up CSV file and header information for the file
csv_filename = "wrfout_TRHWindPress2.csv"
system("rm -f " + csv_filename) ; Remove file in case it exists.
fields = (/"TIME", "LAT", "LON","slp", "TEMPERATURE (degC)","RH2","wspd","wdir"/)
;---Create a header line for CSV file
dq = str_get_dq()
fields = dq + fields + dq ; Pre/append quotes to field names
header = [/str_join(fields,",")/] ; Header is field names separated
; by commas.
;
; Format to use for writing each variable to CSV file.
; If you don't want spaces in CSV file, use the following
; format string:
; format = "%s,%g,%g,%g"
;
format = "%s,%6.2f,%7.2f,%6.2f,%6.2f,%7.2f,%6.2f,%6.2f"
;
; Loop through each time step and desired list of lat/lon values,
; and write a single line of data to CSV file.
;
write_table(csv_filename, "w", header, "%s") ; Write header to CSV file.
do it = 0,ntimes-1
do nl = 0,nlatlon-1
nln = loc2(0,nl)
nlt = loc2(1,nl)
lat1 = a[0]->XLAT(0,nlt,nln) ; nearest grid point
lon1 = a[0]->XLONG(0,nlt,nln)
wspd = wind_speed(u10(:,nlt,nln), v10(:,nlt,nln))
opt = -999
wdir = wind_direction(u10(:,nlt,nln), v10(:,nlt,nln), opt)
alist = [/times(it),lat1,lon1,slp(it,nlt,nln),tc2(it,nlt,nln),rh2(it,nlt,nln),wspd(it),wdir(it)/] ; Store data to be written in a list.
write_table(csv_filename, "a", alist, format) ; Write list to CSV file.
end do
end do
end