; Plot(2) is Only the CONVECTIVE PIXELS load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" begin Plot = new(4,graphic) ; 1. Read file and Print Meta Info res = True res@mpMinLatF = 13.0 res@mpMaxLatF = 14.0 res@mpMinLonF = 75.0 res@mpMaxLonF = 76.0 res@mpFillOn = False res@mpOutlineBoundarySets = "National" res@mpGeophysicalLineColor = "Red" res@mpGeophysicalLineThicknessF= 3.5 res@mpOceanFillColor = 8 res@gsnDraw = True res@gsnFrame = True res@gsnMaximize = True ; ::::::: TICK MARK:::: res@tmXTOn = False res@tmYROn = False res@tmXBMode = "Explicit" ; res@tmXBValues = (/"73","74","75","76","77","78","79","80","81","82","83"/) ; res@tmXBLabels = (/"73~S~o~N~E","74~S~o~N~E","75~S~o~N~E","76~S~o~N~E","77~S~o~N","78~S~o~N~E","79~S~o~N~E","80~S~o~N~E","81~S~o~N~E","82~S~o~N~E","83~S~o~N~E"/) res@tmXBValues = (/"75","75.1","75.2","75.3","75.4","75.5","75.6","75.7","75.8","75.9","76"/) res@tmXBLabels = (/"75~S~o~N~E","75.1~S~o~N~E","75.2~S~o~N~E","75.3~S~o~N~E","75.4~S~o~N","75.5~S~o~N~E","75.6~S~o~N~E","75.7~S~o~N~E","75.8~S~o~N~E","75.9~S~o~N~E","76~S~o~N~E"/) res@tmXBMajorLengthF = 0.01 res@tmXBMinorLengthF = 0.01 res@tmYLMode = "Explicit" ; res@tmYLValues = (/"10","11","12","13","14","15","16","17","18","19"/) ; res@tmYLLabels = (/"10~S~o~N~N","11~S~o~N~N","12~S~o~N~N","13~S~o~N","14~S~o~N~N","15~S~o~N~N","16~S~o~N~N","17~S~o~N~N","18~S~o~N~N","19~S~o~N~N"/) res@tmYLValues = (/"13","13.1","13.2","13.3","13.4","13.5","13.6","13.7","13.8","13.9","14"/) res@tmYLLabels = (/"13~S~o~N~N","13.1~S~o~N~N","13.2~S~o~N~N","13.3~S~o~N","13.4~S~o~N~N","13.5~S~o~N~N","13.6~S~o~N~N","13.7~S~o~N~N","13.8~S~o~N~N","13.9~S~o~N~N","14~S~o~N~N"/) res@tmYLMajorLengthF = 0.01 res@tmYLMinorLengthF = 0.01 res@tmXBMinorOn = False res@tmYLMinorOn = False ; ;:::::::::::::::::::::::::: RES for PLOT(0) ; Rmaxz = res Rmaxz@gsnRightString = "Max Z" Rmaxz@cnFillOn = True Rmaxz@cnFillMode = "RasterFill" Rmaxz@cnLevelSelectionMode = "ExplicitLevels" Rmaxz@cnLevels = (/17,20,36,40,44,52/) ; Colors_res = (/"white","LemonChiffon","LightSkyBlue","medium spring green","maroon","green yellow","Yellow","Red"/) ; <17 17-20 20-28 28-36 36-40 40-44 44-52 , >52 ; Rmaxz@cnFillPalette = Colors_res Rmaxz@cnFillPalette = "CBR_wet" ; :::::::RES for the PLOT(1) ; Rtyperes = res Rtyperes@gsnRightString = "Rain type" Rtyperes@cnFillOn = True Rtyperes@cnFillMode = "RasterFill" Rtyperes@cnLevelSelectionMode = "ExplicitLevels" Rtyperes@cnLevels = (/100,170,200,292,300/) ; Colors = (/"white","green yellow","green yellow","red","Lightgray","white"/) ; <100 , 100-170 , 170-200 , 200-292, 292-300 , >300 ; Rtyperes@cnFillPalette = Colors Rtyperes@cnFillPalette = "CBR_wet" ;::::::::RES for PLOT(2):::: ; Resconv = Rmaxz Resconv@gsnRightString = "Convective pixels" Resconv@cnFillPalette = "CBR_wet" Resconv@cnFillMode = "RasterFill" ; ;::::::::::::::::: Read File n date::::::::::::::::::::::: ; DATE_Of_FILE = "2A25.20110414.76394" file_name = "./DATA_FNL/"+DATE_Of_FILE+".7.HDF" extract_date = str_split(file_name,".") ; DATE print(extract_date) hdf4_file = addfile(file_name, "r") Rmaxz@gsnLeftString = "(a) "+extract_date(1) ; DATE res for PLOT(0) Rtyperes@gsnLeftString = "(b) "+extract_date(1) Resconv@gsnLeftString = "(c) "+extract_date(1) longitude = hdf4_file->Longitude latitude = hdf4_file->Latitude longitude@units = "degrees_east" latitude@units = "degrees_north" X = getfiledimsizes(hdf4_file) ; print (X) Total_Npix = X(0) Total_Nscan = X(2) Total_Level = X(5) WKS_FIG4 = gsn_open_wks("png","test."+extract_date(1)) ; ; Read data. Radar obs "Reflectivity" is short ;---------------------------------------------------- Z_short = hdf4_file->correctZFactor Z_short@_FillValue = toshort(-9999) Z_short = where(Z_short.lt. 0, Z_short@_FillValue,Z_short) Z_fnl = Z_short/(Z_short@scale_factor) copy_VarCoords(Z_short,Z_fnl) Z_fnl@long_name = Z_short@hdf_name Z_fnl@units = Z_short@units ; printMinMax(Z_fnl,True) ; :::::::::::::::::::: ; 2nd Variable ;:::::::::::::::::::::: ; Raintype_short = hdf4_file->rainType ; printVarSummary(Raintype_short) Raintype_flt = tofloat(Raintype_short) ; printVarSummary(Raintype_flt) Raintype_flt@_FillValue = -9999 Raintype_float = where(Raintype_flt .le. -88, Raintype_flt@_FillValue,Raintype_flt) Raintype_float@_FillValue = -9999 ; printVarSummary(Raintype_float) ; print(Raintype_float(3162,27)) ; print("HEY") Raintype_float!0 = "longitude" Raintype_float!1 = "latitude" Raintype_float@lon2d = longitude Raintype_float@lat2d = latitude ; Plot(1) = gsn_csm_contour_map(WKS_FIG4,Raintype_float,Rtyperes) ; plot of Rain type ; ;-:::::::::::::::::::::::::::::::::::::::::::::::::: ; Find out the max value of Z amongst 80 levels ;:::::::::::::::::::::::::::::::::::::::::::::::::::: ; maxZ = dim_max_n(Z_fnl,2) ; Z_fnl is 3 dim. It s finding max Z at each level corresponding to each pixel ; printVarSummary(maxZ) ; print(maxZ(2915,28)) maxZ!0 = "longitude" maxZ!1 = "latitude" maxZ@lon2d = longitude maxZ@lat2d = latitude ; Plot(0) = gsn_csm_contour_map(WKS_FIG4,maxZ,Rmaxz) ; Plot of maximum Z in vertical ; ; CHK for THE RES (Rmaxz) for this PLOT UP. ;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; Now filter max_Z acording to the Raintype >300 and <170 are to be removed ;::::::::::::::::::::::::::::::::::::::::::::::::::::::: maxZ1 = where(Raintype_float .ge.(300),maxZ@_FillValue,maxZ) ; max_Z1 is when raintype >300 put fill value ; printVarSummary(maxZ1) ; printMinMax(maxZ1,True) maxZ2 = where(Raintype_float .le.(170),maxZ@_FillValue,maxZ1) ; Less than 170 data removed ; printVarSummary(maxZ2) ; printMinMax(maxZ2,True) maxZ_Conv = where(maxZ.lt.(17),maxZ@_FillValue,maxZ2) ; less than 17dBz removed ; printVarSummary(maxZ_Conv) ; printMinMax(maxZ_Conv,True) ;-------| maxZ_Conv is the final variable in which ge 300, le 170 and 17 dBz has been removed | ;-------|-------------------------------------------------------------------------------------| ; maxZ_Conv!0 = "longitude" maxZ_Conv!1 = "latitude" maxZ_Conv@lon2d = longitude maxZ_Conv@lat2d = latitude Plot(2) = gsn_csm_contour_map(WKS_FIG4,maxZ_Conv,Resconv) ; plot of only convective pixels ; ------Check Mary 's mail on 27.12.2016---------- ; ---------------- mkres = res mkres@gsMarkerSizeF = 2 mkres@gsMarkerIndex = 16 mkres@gsnCoordsNonMissingColor = "black" mkres@gsnCoordsMissingColor = "red" gsn_coordinates(WKS_FIG4,Plot(2),maxZ_Conv,mkres) print(Z_fnl(6352,17,:)) delete(Z_fnl) delete(Raintype_float) delete(Raintype_flt) delete(maxZ1) delete(maxZ2) NSstart = 6239 NSend = 6669 NPstart = 0 NPend = Total_Npix-1 do NS = NSstart, NSend do NP = NPstart, NPend NSS = NS NPP = NP Lat = latitude(NS,NP) Lon = longitude(NS,NP) if (.not.(ismissing(maxZ_Conv(NS,NP)))) then print("NSS = "+NSS+" NPP = "+NPP+" Lat = "+Lat+" Lon = "+Lon+" maxZ_Conv = "+maxZ_Conv(NS,NP)) end if end do end do end