提取WRF的nc数据

今天又破案了,应该用wrfout数据与观测数据进行对比,具体是:
气温:2米处的气温

T2.png

相对湿度:2米的相对湿度
wrf_user_getvar (ucar.edu)
RH2.png

风:10米处的风
U10和V10.png

气压:海平面气压
slp.png

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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 200,667评论 5 472
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,361评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 147,700评论 0 333
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,027评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,988评论 5 361
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,230评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,705评论 3 393
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,366评论 0 255
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,496评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,405评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,453评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,126评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,725评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,803评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,015评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,514评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,111评论 2 341