- 积分
- 702
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-9-29
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2018-10-23 08:04:18
|
显示全部楼层
脚本如下:
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/csm/contributed.ncl"
begin
cmap = (/(/255,255,255/),\
(/0,0,0/),\
(/255,255,255/),\
(/186,222,130/),\
(/145,222,7/),\
(/102,209,10/),\
(/30,186,12/),\
(/255,255,0/),\
(/255,166,0/),\
(/255,0,0/)/) /255.0
filename = "3B-HHR.MS.MRG.3IMERG.20170714-S233000-E235959.1410.V05B.HDF5"
h5_file = addfile(filename,"r")
lat1 = 3.86
lat2 = 53.55
lon1 = 73.66
lon2 = 135.03
data_raw = h5_file->/NS/SLV/precipRateESurface
data_raw@_FillValue = data_raw@CodeMissingvalue
printVarSummary(data_raw)
printMinMax(data_raw,0)
longitude=h5_file->/NS/Longitude
latitude=h5_file->/NS/Latitude
longitude@units = "degrees_east"
latitude@units = "degrees_north"
nsizescan=dimsizes(data_raw(:,0))
nsizeray=dimsizes(data_raw(0,:))
print(nsizescan)
print(nsizeray)
sumdata=0.0
do ns=0,nsizescan-1
do nr=0,nsizeray-1
if (latitude(ns,nr).ge.lat1.and.latitude(ns,nr).le.lat2) then
if (longitude(ns,nr).ge.lon1.and.longitude(ns,nr).le.lon2) then
sumdata=sumdata+data_raw(ns,nr)
end if
end if
end do
end do
print(sumdata)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; yearb=h5_file->/NS/ScanTime/Year(0)
; monb=h5_file->/NS/ScanTime/Month(0)
; dayb=h5_file->/NS/ScanTime/DayOfMonth(0)
; hourb=h5_file->/NS/ScanTime/Hour(0)
; minuteb=h5_file->/NS/ScanTime/Minute(0)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; yeare=h5_file->/NS/ScanTime/Year(nsizescan-1)
; mone=h5_file->/NS/ScanTime/Month(nsizescan-1)
; daye=h5_file->/NS/ScanTime/DayOfMonth(nsizescan-1)
; houre=h5_file->/NS/ScanTime/Hour(nsizescan-1)
; minutee=h5_file->/NS/ScanTime/Minute(nsizescan-1)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; date_b_str = sprinti("%0.4i", yearb)+sprinti("%0.2i", monb)+sprinti("%0.2i", dayb)+\
; sprinti("%0.2i", hourb)+sprinti("%0.2i", minuteb)
; print(date_b_str)
;
; date_e_str = sprinti("%0.4i", yeare)+sprinti("%0.2i", mone)+sprinti("%0.2i", daye)+\
; sprinti("%0.2i", houre)+sprinti("%0.2i", minutee)
; print(date_e_str)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
data = data_raw
data@long_name = "precipRate"
data@_FillValue = -9999.9
data@units = data_raw@Units
data@lat2d = latitude
data@lon2d = longitude
printVarSummary(data)
wks = gsn_open_wks("png", "GPM"+"PRE") ; open workstation
gsn_define_colormap(wks,cmap)
res=True
res@gsnMaximize=True ;make plot large
; res@gsnPaperOrientation = "portrait" ;force portrait orientation
;;set map
res@mpSpecifiedFillColors = (/"white","transparent","white"/) ;鑳屾櫙
res@mpLandFillColor = -1 ; transparent
res@mpOutlineOn = True
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "Earth..4"
res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@mpOutlineSpecifiers = (/"China:states","Taiwan"/)
res@mpGeophysicalLineThicknessF= 2.
res@mpNationalLineThicknessF= 2.
res@pmTickMarkDisplayMode = "Always"
; now change the size of the tickmark labels
res@tmXBLabelFontHeightF = 0.02 ; resize tick labels
res@tmYLLabelFontHeightF = 0.02
; now change the thickness of the X-Y axis
res@tmBorderThicknessF=4
res@tmXBMajorThicknessF=4
res@tmXBMajorThicknessF=4
res@tmYLMajorThicknessF=4
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
res@gsnMaximize = True ; maximize plot in frame
res@cnFillOn = True ; turn on contour fill
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False ; turn off line labels
res@gsnAddCyclic = False ; set to False if plotting ;regional ;data
;res@tiMainString = times(nt)+"_1h precipitation"
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnFillMode="RasterFill" ;faster
;rescources@cnMinLevelValF = 0.1
;rescources@cnMaxLevelValF = 250
;rescources@cnLevelSpacingF=0.1-250
res@cnLevels = (/1,2.5,5,10,20,30,40/)
;rescources@cnFillColors =(/62,92,97,80,83,87,18,21,23,48,52,54,36,41/)
res@lbLabelAutoStride= True
res@lbOrientation="vertical" ;vertical labels
res@trGridType = "TriangularMesh"
; res@tiMainString = "GPM"+"PRE"+date_b_str+"-"+date_e_str
res@mpMinLatF = lat1
res@mpMaxLatF = lat2
res@mpMinLonF = lon1
res@mpMaxLonF = lon2
res@lbLabelBarOn = True
res@gsnLeftStringFontHeightF = 12 ; make font smaller
res@gsnRightStringFontHeightF = 12 ; make font smaller
res@gsnLeftString = "surfacePrecipitationRate"
plot=gsn_csm_contour_map_ce(wks,data,res)
mkres = True
mkres@gsMarkerIndex = 2
mkres@gsMarkerColor = "black"
draw(plot)
end |
|