中国气象局grd数据，CMORPH逐小时做图

woshiwzm

和官网给的图对比完全一致，借鉴了前人做的0.5*0.5的数据，改一改就成0.1*0.1的了

1. ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
3. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   
4. load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
   
5. begin
   
6.   
   
7.   ;fName = "2019082000.grd"  
   
8.   myfile = systemfunc("ls *.grd")
   
9.   nfile = dimsizes(myfile)
   
10.                      
    
11.   nlat  = 450                               ; YDEF
    
12.   mlon  = 700
    
13.   ntime = nfile
    
14.                             ; XDEF
    
15. 
    
16.                                           ; create an array to contain data
    
17. rain   = new ( (/nlat,mlon/), double)
    
18. hour_rain = new ( (/ntime,nlat,mlon/), double)
    
19.   ;rain@units     = "0.1mm"
    
20. 
    
21.     do i=0, nfile-1                              
    
22.      setfileoption("bin","ReadByteOrder","LittleEndian")
    
23.                                          
    
24.      rain(:,:) = fbindirread(myfile(i),0, (/nlat,mlon/), "float")
    
25.      hour_rain(i,:,:) = rain(:,:)
    
26.      rain = 0
    
27.      ;print(myfile(i))
    
28.     end do
    
29.    
    
30.   printVarSummary(hour_rain)
    
31.   printMinMax(hour_rain, 1)
    
32. 
    
33. lat=15+ispan(0,nlat-1,1)*0.1
    
34. lon=70+ispan(0,mlon-1,1)*0.1
    
35. time = ispan(0, nfile-1, 1)
    
36. 
    
37. lat!0= "lat"
    
38. lat@long_name = "latitude"
    
39. lat@units = "degrees_north"
    
40. lon!0= "lon"
    
41. lon@units = "degrees_east"
    
42. lon@long_name = "longitude"
    
43. rain!0 = "lat"
    
44. rain!1 = "lon"
    
45. rain&lat = lat
    
46. rain&lon = lon
    
47. rain@_FillValue =-999
    
48. time!0= "time"
    
49. time@long_name = "time"
    
50. time@units = "hour"
    
51. hour_rain!0 = "time"
    
52. hour_rain!1 = "lat"
    
53. hour_rain!2 = "lon"
    
54. hour_rain&time = time
    
55. hour_rain&lat = lat
    
56. hour_rain&lon = lon
    
57. hour_rain@_FillValue =-999  
    
58. sumrain = dim_sum_n(hour_rain, 0)
    
59. copy_VarMeta(rain, sumrain)
    
60. printMinMax(sumrain, 1)
    
61. delete(hour_rain)
    
62. ;print(lat)
    
63. ;smthrain = smth9(sumrain, 0.5, 0.25, False)
    
64. ;copy_VarMeta(rain, smthrain)
    
65. delete(rain)
    
66.   wks = gsn_open_wks("pdf","00")     ; Open a workstation
    
67.   gsn_define_colormap(wks,"precip3_16lev")      ; define a different colormap.
    
68.   res = True
    
69.   res@gsnAddCyclic  = False      
    
70.   res@gsnDraw               = False
    
71.   res@gsnFrame = False
    
72.   res@gsnRightString = "mm/day"
    
73.       
    
74.   ; Don't draw yet
    
75.   res@mpDataSetName         = "Earth..4"   ; This new database contains
    
76.                                            ; divisions for other countries.
    
77.   res@mpDataBaseVersion     = "MediumRes"  ; High resolution database
    
78.   res@mpOutlineOn           = True         ; Turn on map outlines
    
79.   res@mpOutlineSpecifiers   = (/"China:states","Taiwan"/)       ;China:states
    
80.   res@mpGeophysicalLineThicknessF= 2.      ; double the thickness of geophysical boundaries
    
81.   res@mpNationalLineThicknessF= 2.         ; double the thickness of national boundaries
    
82.   ;res@mpProjection = "LambertConformal"
    
83.   res@mpLambertMeridianF = 105.0
    
84.   res@mpLimitMode = "LatLon"
    
85.   res@cnFillOn = True
    
86.   res@cnSmoothingOn = True
    
87.   res@cnSmoothingDistanceF = 0.01
    
88.   res@cnSmoothingTensionF = 2.5
    
89.   res@cnLevelSelectionMode = "ExplicitLevels"
    
90.   res@cnLevels = (/0.1,10,20,30,40,50,60,70,80,90,100,110/)
    
91.   ;res@cnFillColors = (/0,128,96,64,48,32,144,160,176,192,208,244,255/)
    
92.   res@mpMinLatF             =  20        ; Asia limits
    
93.   res@mpMaxLatF             =  30
    
94.   res@mpMinLonF             =  100
    
95.   res@mpMaxLonF             =  112
    
96. 
    
97.   res@mpAreaMaskingOn = True   
    
98.   res@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/)   ;China:states
    
99.   res@mpOceanFillColor = 0  
    
100.   res@mpInlandWaterFillColor = 0  
     
101.   res@cnFillOn      = True
     
102.   res@cnFillMode="AreaFill"
     
103.   ;res@cnRasterSmoothingOn = True
     
104.   res@cnConstFEnableFill=False
     
105.   res@cnLinesOn     =False
     
106.   res@cnLineLabelsOn = False
     
107.   ;res@cnFillDrawOrder = "PreDraw"         ; draw contours first   
     
108.   res@lbLabelBarOn = True       ;
     
109.   res@tiMainFont            = "helvetica"
     
110.   res@tiMainOffsetYF        = 0.02  ;set place for main title along   Y,offset
     
111.   res@tiMainFontHeightF     = 0.02   ;set main title font size
     
112.   res@tiMainString          = "rain"
     
113.   plot= gsn_csm_contour_map(wks,sumrain,res)
     
114.   lat1 = (/25.7,27.3,27.3,25.7,25.7/)
     
115.   lon1 = (/103,103,105.5,105.5,103/)
     
116.   plres = True
     
117.   plres@gsLineThicknessF = 2.0
     
118.   plres@gsLineColor = "grey20"
     
119.   plot_1 = gsn_add_polyline(wks, plot, lon1, lat1, plres)
     
120.   draw(plot)
     
121. 
     
122.   ;frame(wks)
     
123. 
     
124.    end
     
125. 
     

复制代码

pokerface

大佬，代码我看了，不怎么懂，想问一下如何读取这个数据每个网格用到的站点数？我在matlab中只能读取降水数据

woshiwzm

pokerface 发表于 2020-5-31 22:13
大佬，代码我看了，不怎么懂，想问一下如何读取这个数据每个网格用到的站点数？我在matlab中只能读取降水数 ...

这不是站点数据，这是网格数据

立夏

做出来只有一张空白图，没有降水数据呢

woshiwzm

立夏 发表于 2020-11-16 10:31
做出来只有一张空白图，没有降水数据呢

你的图片格式是什么。png要frame（wks）

啦啦海绵

{:5_213:}

啦啦海绵

谢谢楼主的分享！

小迪

请问中国自动站和cmorph融合降水资料，最新到哪个时段？谢谢，想要2020年10月份的数据。

wangwangxiao

感谢楼主分享

kare

pokerface 发表于 2020-5-31 22:13
大佬，代码我看了，不怎么懂，想问一下如何读取这个数据每个网格用到的站点数？我在matlab中只能读取降水数 ...

大佬，matlab怎么打开这个grd的数据啊，我读取怎么都是乱码啊