;---------------------------------------------- ;CALIPSO Test Code for Total Backscatter at 532 nm ;---------------------------------------------- load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl" begin ;---------------------------------------------- ; OPEN CALIPSO FILE AND WORKSTATION ;---------------------------------------------- dir_name = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/Backscatter_Test_Data/" file_name = "CAL_LID_L2_05kmAPro-Prov-V3-01.2010-01-30T16-21-48ZD_Subset.hdf" ;dir_name = "./" ;file_name = "CAL_LID_L2_05kmAPro-Prov-V3-01.2006-06-13T01-37-01ZN.hdf" hdf4_file = addfile(dir_name+file_name, "r") fig_diri = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/Backscatter_Test_Images/" ;directory that image gets saved in ;fig_diri = "./" fig_save = fig_diri + "CALIPSO_Backscatter_Test_2010-01-30T16-21-48ZD"; will save this image in the directory specified (fig_diri) wks = gsn_open_wks("png", fig_save) ; open workstation & send graphics to PNG file. Saves graphics with file as indicated above. ;---------------------------------------------- ;Print information about the file to know what variables and attributes are ;available for plotting. ;---------------------------------------------- print(hdf4_file) ;---------------------------------------------- ; READ CALIPSO DATA ;---------------------------------------------- TBC_532 = hdf4_file->Total_Backscatter_Coefficient_532 ; [fakeDim70 | 3001] x [fakeDim71 | 399] printVarSummary(TBC_532) ; Info: ; Variable: TBC_532 TBC_532 = TBC_532(:,::-1) ; I most likely need to reverse this 2D array (altitude is also reversed when viewing text file) TBC_532@_FillValue = -9999.0 ;---------------------------------------------- ; GET LATITUDE AND LONGITUDE FOR THE X-AXIS AND FORMAT THE LABELING ;---------------------------------------------- Latitude = hdf4_file->Latitude Latitude@units = "degrees_north" ;printVarSummary(Latitude) Longitude = hdf4_file->Longitude Longitude@units = "degrees_east" ;printVarSummary(Longitude) xlabel = sprintf("%.2f",Latitude)+"~C~"+sprintf("%.2f",Longitude) ; format labeling on x-axis ;---------------------------------------------- ; GET TIME FIELDS TO GET A TIME STRING, WHICH WILL BE USED TO STRIDE X-AXIS LATITUDE AND LONGITUDE ;---------------------------------------------- ; stride might be time step? ask Dennis Shea for clarification time = (/hdf4_file->Profile_Time/) time@units = "seconds since 1993-01-01 00:00" tstring = ut_string(time(:,0), "%Y-%N-%D %H:%M:%S") stride = dimsizes(tstring)/10 ;printVarSummary(time) ;printVarSummary(tstring) ;---------------------------------------------- ; GET ALTITUDE VALUES; Figure out a way to extract metadata ;---------------------------------------------- ;asciiread function: If you don't know how many data values you have, you can use the special "-1" value for the dimension size. ;When you use -1, data values will be read from left-to-right, top-to-bottom, ;into a 1D array, until there are no values left. ;once converted into a text file, it should theoretically be accessed through: ;---------------------------------------------- ; set a directory above this and write the loop hgt = asciiread("/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/Backscatter_Test_Data/altitudes.txt",-1,"float") ; hgt = height ; hgt = asciiread("altitudes.txt",-1,"float") ; hgt = height ; [399] hgt = hgt(::-1) ; Reverse this array hgt!0 = "hgt" ; first dimension being referenced, which is altitude. hgt@long_name = "Altitude, km" hgt@units = "km" printVarSummary(hgt) ;Info: ;---------------------------------------------- ; ASSIGN NEW DIMENSIONS AND ATTRIBUTES TO COORDINATE CALIPSO Total Backscatter Coefficient WITH ALTITUDE ;---------------------------------------------- TBC_532!1 = "hgt" ; assigning hgt as first dimension TBC_532&hgt = hgt ; xcoord = ispan(0,dimsizes(Latitude(:,0)) - 1,1) ;ispan creates an array of equally-spaced values. ;dimsizes returns the dimension sizes of the input variable. Confirm if this command is correct. TBC_532!0 = "xcoord" TBC_532&xcoord = xcoord ;---------------------------------------------- ; PLOTTING RESOURCES ;---------------------------------------------- setvalues NhlGetWorkspaceObjectId() "wsMaximumSize": 300000000 end setvalues colorMap = read_colormap_file("BlAqGrYeOrRe") colorMap(0,:) = (/ 0., 0., 0., 1./) res = True res@cnFillOn = True ; color plot desired res@cnLinesOn = False res@cnFillPalette = colorMap res@cnRasterModeOn = True res@gsnAddCyclic = False res@gsnMaximize = True res@gsnLeftString = "UTC: "+tstring(0)+" to "+tstring(dimsizes(tstring)-1) res@tiMainString = "Total Backscatter Coefficient 532nm" res@tiMainFontHeightF = 0.025 res@tiYAxisFontHeightF = 0.02 res@sfYArray = hgt ;--------------------------------------------- ; HERE, SET SPECIAL LABEL HEIGHTS AND ORIENTATION ;--------------------------------------------- res@lbLabelFontHeightF = 0.01 res@lbLabelAngleF = 0.0 res@lbOrientation = "vertical" res@vpWidthF = 1.0 ;--------------------------------------------- ; HERE YOU CAN SET A MAXIMUM ALTITUDE ;--------------------------------------------- res@trYMaxF = 5 ; THIS IS IN KILOMETERS. Needed? Check once plot has been constructed ;--------------------------------------------- ; SET SPECIAL LEVELS ;--------------------------------------------- res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/ 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.0009, \ 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, 0.005, \ 0.0055, 0.006, 0.0065, 0.007, 0.0075, 0.008, \ 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 /) ;levels from example NCL script. Appears to be correct from Browse Images section in Subsetter Tool ;--------------------------------------------- ; HERE IS WHERE YOU STRIDE THE DATA AND PLOT LATITUDE AND LONGITUDE ON X-AXIS ;--------------------------------------------- ;tm = tickmark res@tmXTOn = False ; turns top tick marks off res@tmXBMode = "Explicit" ;uses the array resource tmXBValues to determine the coordinate positions of the tick marks res@tmXBValues = xcoord(::stride) ;(::stride) is used to subsample an array. subsamples xcoord. res@tmXBLabels = xlabel(::stride,0) res@tmXBLabelFontHeightF = 0.01 ; toggle this if you want ;--------------------------------------------- ; HERE, SET SPECIAL LABEL HEIGHTS AND ORIENTATION ;--------------------------------------------- res@lbLabelFontHeightF = 0.01 res@lbLabelAngleF = 0.0 res@lbOrientation = "vertical" res@vpWidthF = 1.0 ;--------------------------------------------- ; PLOT ;--------------------------------------------- plot = new(2,"graphic") ;Code below is for map and swath res@gsnDraw = True res@gsnFrame = False ;printVarSummary(TBC_532) ; [xcoord | 3001] x [hgt | 399] plot(0) = gsn_csm_contour(wks,TBC_532(hgt|:,xcoord|:),res) ; make 'hgt' the Y-Axis. Pipe (|:) required for re-ordering frame(wks) ;delete([/Latitude,Longitude, TBC_532, xcoord, xlabel/]) -- probably important once loop is constructed ;--------------------------------------------- ; MAP ;--------------------------------------------- ;n = dimsizes(file_name) ;calls out the number of files mpres = True mpres@gsnDraw = False ; False so that it doesn't create a new image everytime and name it something like ".00002" etc. mpres@gsnFrame = False ; mpres@mpLimitMode = "LatLon" ; mpres@mpMinLatF = -88 ; mpres@mpMinLonF = -100 ; mpres@mpMaxLatF = -40 ; mpres@mpMaxLonF = -20 latitude = hdf4_file->Latitude longitude = hdf4_file->Longitude plot(1) = gsn_csm_map(wks,mpres) plres = True plres@gsLineColor = "red" lnid = gsn_add_polyline(wks,plot(1),longitude(:,0),latitude(:,0),plres) resP = True resP@gsnPanelCenter = True gsn_panel(wks,plot,(/2,1/),resP) ;res2 = True end