;---------------------------------------------- ; CALIPSO VFM FIGURE LOOP ;---------------------------------------------- ; Copyright (C) 2014-2017 by The HDF Group ; All Rights Reserved ; ; This example code illustrates how to access and visualize LaRC CALIPSO Lidar ; Level 2 Vertical Feature Mask Version 4.10 HDF4 file in NCL. ; ; If you have any questions, suggestions, comments on this example, please use ; the HDF-EOS Forum (http://hdfeos.org/forums). ; ; If you would like to see an example of any other NASA HDF/HDF-EOS data ; product that is not listed in the HDF-EOS Comprehensive Examples page ; (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org or ; post it at the HDF-EOS Forum (http://hdfeos.org/forums). ; ; Tested under: NCL 6.4.0 ; Last updated: 2017-05-05 ;---------------------------------------------- ; LOAD LIBRARIES ;---------------------------------------------- load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl" ;---------------------------------------------- ; START OF CODE ;---------------------------------------------- begin ;---------------------------------------------- ; SET DIRECTORIES ;---------------------------------------------- file_hdf = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/2007_to_2017_January_PalmerLTERDATA_for_Overlay_and_VFM/2010_January_PalmerLTERDATA_Overlay_Data/" ; input directory file_name = systemfunc("cd "+file_hdf+" ; ls CAL*.hdf") printVarSummary(file_name) nfile = dimsizes(file_name) print(file_hdf) print("---") print("nfile="+nfile) print("---") ;---------------------------------------------- ; FIGURE DESTINATION ;---------------------------------------------- fig_diri = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/2007_to_2017_January_PalmerLTER_VFM_Figures/2010_January_PalmerLTER_VFM_Figures/" ;directory that image gets saved in ;---------------------------------------------- ; LOOP START ;---------------------------------------------- do nf = 0,nfile-1 print("==============================================") hdf4_file = addfile(file_hdf+file_name(nf), "r") dat = str_get_cols(file_name(nf),31,40) ; to obtain the date string print(dat) ;---------------------------------------------- ; RENAME AND SAVE FIGURE IN FIGURE DIRECTORY ;---------------------------------------------- fig_save = fig_diri + file_name(nf) ;---------------------------------------------- ; OPEN WORKSPACE ;---------------------------------------------- wks = gsn_open_wks("png", fig_save + ".v.ncl") ; open workstation ;---------------------------------------------- ; NOTES from NCL Site ;---------------------------------------------- ; Print information about the file to know what variables and attributes are ; available for plotting. ; print(hdf4_file) ; Users need to understand the layout of the Feature_Classification_Flag ; dataset. ; ; The Feature_Classification_Flag values are stored as a sequence of 5515 ; element arrays (i.e., as an N x 5515 matrix, where N is the number of ; separate records in the file). In this file, N is 4224. ; ; Each array represents a 5 km "chunk" of data, ; and each array element contains the feature classification information for a ; single range resolution element in the Level 0 lidar data downlinked from ; the satellite. As shown in the summary below, the vertical and horizontal ; resolution of the CALIPSO data varies as a function of altitude above mean ; sea level (see Hunt et al., 2009). ; ; Here's the summary of number of profiles per 5 km. ; ; 3 profiles for 20.2km (base) to 30.1km (top) @ 180m ; (index 1-165 / 55 samples per profile) ; ; 5 profiles for 8.2km (base) to 20.2km (top) @ 60m ; (index 166 - 1165 / 200 samples per profile) ; ; 15 profiles for -0.5km (base) to 8.2km (top) @ 30m ; (index 1166 - 5515 / 290 samples per profile) ; ; 3 profiles mean horizontal resolution is 1667m because 3 * 1667m = 5km. ; 55 samples mean vertical resolution is 180m because 55 * 180m = 9.9km (= ; 30.1km - 20.2km). ; ; In short, profile size determines horizontal resolution and sample size ; determines the vertical resolution. ; ; Each record is an unsigned 16 bit integer. See [1] for details. ; Bits | Description ; ---------------- ; 1-3 | Feature Type ; 4-5 | Feature Type QA ; ... |... ; 14-16 | Horizontal averaging ; ; ; In this example, we'll focus only on Featrue Type (bits 1-3). ; ; There are many possibilites to plot this data. ; ; Here, we'll subset -05km to 8.2km (e.g., 1165:5514) ; over latitude approx. 40N to 62N (e.g., 3500:3999) ; and plot altitude vs. latitude. ;---------------------------------------------- ; NOTES END ;---------------------------------------------- ;---------------------------------------------- ; CLASSIFY FEAUTRE CLASSIFICATION FLAGS ;---------------------------------------------- fcf = hdf4_file->Feature_Classification_Flags ;---------------------------------------------- ; What does this mean? V ; Subset dataset that latitude value increases monotonically. ; Monotonicity is required by NCL. ;---------------------------------------------- ; lat = hdf4_file->Latitude(400:800,0) ;---------------------------------------------- ; Why identify lat above ^ and again below? ;---------------------------------------------- lat = hdf4_file->Latitude(:,0) tmp = dimsizes(lat) print(tmp) ;---------------------------------------------- ; What is going on below?? v ;---------------------------------------------- if (tmp .le. 150) then print("failed to plot" + dat) delete(hdf4_file) delete(fcf) delete(lat) delete(tmp) continue end if delete(lat) lat = hdf4_file->Latitude(50:150,0) ; WHY IDENTIFY LAT HERE? ;---------------------------------------------- ; What is going on above??? ^ ;---------------------------------------------- ; You can visualize other blocks by changing subset parameters. ; profile = fcf(470:695, 0:164) ; 20.2km to 30.1km ; profile = fcf(470:695,165:1164) ; 8.2km to 20.2km profile = fcf(50:150,1165:5514) ; -0.5km to 8.2km ;---------------------------------------------- ; Errors above ^ if chunk of code is edited ;---------------------------------------------- ;---------------------------------------------- ; GET TIME FIELDS TO GET A TIME STRING, WHICH WILL BE USED TO STRIDE X-AXIS LATITUDE AND LONGITUDE ; MUST LOAD IN UTC LIBRARY IN THE BEGGINING OF SCRIPT ;---------------------------------------------- time := (/hdf4_file->Profile_Time/) ; overwrite syntax because each new file in loop may have a different size time@units = "seconds since 1993-01-01 00:00" tstring := ut_string(time(:,0), "%Y-%N-%D %H:%M:%S") stride = dimsizes(tstring)/15 size = dimsizes(profile) ; Select the first 3 bit for the feature type. See [2, 3] for dim_gbits. data2d = dim_gbits(profile, 4, 3, 13, size(1)) ; Reshape 2-D data to 3-D for n profiles * y samples. ; You can visualize other blocks by changing subset parameters. ; data3d = reshape(data2d, (/size(0), 3, 55/)) ; data3d = reshape(data2d, (/size(0), 5, 200/)) data3d = reshape(data2d, (/size(0), 15, 290/)) ; data3d = reshape(data2d, (/size(0), 290, 15/)) ; printVarSummary(data3d) ; Subset horizontally. Variation along longitude (i.e., profile) is very ; small. We'll just pick the first location from each profile. data_raw = data3d(:,0,:) ; data_raw = data3d(:,:,0) ;0 = not determined ;1 = clean marine ;2 = dust ;3 = polluted continental ;4 = clean continental ;5 = polluted dust ;6 = smoke ;7 = other levels = (/0.5,1.5,2.5,3.5,4.5,5.5,6.5/) ;---------------------------------------------- ; IDENTIFY COLORMAP ;---------------------------------------------- colormap= (/"white","blue","orange","green","brown","black","pink","grey"/) gsn_define_colormap(wks,colormap) cmap = gsn_retrieve_colormap(wks) ;toggling different aspects of the figure below res = True ; plot mods desired res@cnFillOn = True ; enable contour fill res@gsnMaximize = True; make plot large res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False; turn off contour line labels res@gsnSpreadColors = True ; use the entire color spectrum res@cnFillMode = "RasterFill" ; faster res@cnMissingValFillPattern = 0 ; missing value pattern is set to "SolidFill" res@cnMissingValFillColor = 0; white color for missing values res@lbLabelAutoStride = True ; ensure no label overlap res@cnFillColors = cmap res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = levels; (/0.5,1.5,2.5,3.5,4.5,5.5,6.5/) res@lbLabelAlignment = "BoxCenters" res@lbLabelStrings = (/"Not Determined","Clean Marine","Dust","Polluted Continental","Clean Continental","Polluted Dust","Smoke","Other"/) res@gsnDraw = False res@gsnFrame = False res@lbOrientation = "Vertical" res@lbLabelPosition = "Right" res@tiMainFontColor = 5 res@tiXAxisFontColor = 5 res@tiYAxisFontColor = 5 res@tmBorderLineColor = 5 res@tmXBMajorLineColor = 5 res@tmXBLabelFontColor = 5 res@tmXBMinorLineColor = 5 res@tmXTLabelFontColor = 5 res@tmXTMajorLineColor = 5 res@tmXTMinorLineColor = 5 res@tmYLLabelFontColor = 5 res@tmYLMajorLineColor = 5 res@tmYLMinorLineColor = 5 res@tmYRLabelFontColor = 5 res@tmYRMajorLineColor = 5 res@tmYRMinorLineColor = 5 ;0 = not determined res@tiMainString = "Vertical Feature Mask WAP" res@gsnStringFontColor = 5 res@tiMainFontHeightF = 0.025 ;1 = clean marine ;2 = dust ;3 = polluted continental ;4 = clean continental ;5 = polluted dust ;6 = smoke ;7 = other ; Create altitude dataset manually like NCL CALIPSO example [4]. ; You can visualize other blocks by changing subset parameters. ; ; hgt = ispan(20200, 30020, 180) ; hgt = ispan(8200, 20140, 60) hgti = ispan(-500, 8170, 30) ; hgti = ispan(8170, -500, -30) hgt = tofloat(hgti) / 1000.0 hgt = hgt(::-1) hgt@long_name = "Altitude (km)" hgt@units = "km" hgt!0 = "hgt" lat!0 = "lat" lat@long_name = "Latitude" lat@units = "degrees_north" ; Select "Cloud" feature type. data = todouble(data_raw) ;data@long_name = "Feature Type (Bits 10-12) in Feature Classification Flag" data!0 = "lat" data!1 = "hgt" data&lat = lat ;print(lat) data&hgt = hgt ;---------------------------------------------- ; SET UP PANEL PLOTS ;---------------------------------------------- plot = new(2,"graphic") ;---------------------------------------------- ; PLOT VFM ;---------------------------------------------- plot(0) = gsn_csm_contour(wks, transpose(data(::-1,::-1)),res) ;Transpose data so that x-axis becomes latitude and y-axis becomes altitude. ;---------------------------------------------- ; MAP RESOURCES ;---------------------------------------------- res2 = True ;LOOK UP mpres = True mpres@gsnMaximize = True ; Maximize plot in window. mpres@gsnDraw = False mpres@gsnFrame = False mpres@tmXBLabelFontHeightF = 0.02 mpres@tmYLLabelFontHeightF = 0.02 mpres@pmTickMarkDisplayMode = "Always" mpres@tmBorderLineColor = 5 mpres@tmXBLabelFontColor = 5 mpres@tmXBMajorLineColor = 5 mpres@tmXBLabelFontColor = 5 mpres@tmXBMinorLineColor = 5 mpres@tmXTLabelFontColor = 5 mpres@tmXTMajorLineColor = 5 mpres@tmXTMinorLineColor = 5 mpres@tmYLLabelFontColor = 5 mpres@tmYLMajorLineColor = 5 mpres@tmYLMinorLineColor = 5 mpres@tmYRLabelFontColor = 5 mpres@tmYRMajorLineColor = 5 mpres@tmYRMinorLineColor = 5 mpres@tiMainFontHeightF = 0.02 mpres@tiDeltaF = 0.02 mpres@tiMainFontColor = 5 mpres@tiMainString = "UTC: "+tstring(0)+" to "+tstring(dimsizes(tstring)-1) ;---------------------------------------------- ; LATITUDE AND LONGITUDE BELOW. LIMIT MAP AREA ;---------------------------------------------- mpres@mpLimitMode = "LatLon" mpres@mpMinLatF = -71 mpres@mpMaxLatF = -58 mpres@mpMinLonF = -80 mpres@mpMaxLonF = -58 mpres@gsnDraw = False ; Look up mpres@gsnFrame = False ; Look up ;---------------------------------------------- ; POLYLINE (SWATH) RESOURCES ;---------------------------------------------- plres = True plres@gsLineColor = "red" latitude = hdf4_file->Latitude longitude = hdf4_file->Longitude ;---------------------------------------------- ; PLOT MAP AND SWATH ;---------------------------------------------- plot(1) = gsn_csm_map(wks,mpres) lnid = gsn_add_polyline(wks,plot(1),longitude(:,0),latitude(:,0),plres) resP = True resP@gsnPanelCenter = True resP@tmBorderLineColor = 5 gsn_panel(wks,plot,(/2,1/),resP) ;---------------------------------------------- ; CLEAN UP RESOURCES ;---------------------------------------------- delete(plot) delete(wks) delete(data) delete(res) delete(hdf4_file) delete(fcf) delete(lat) delete(latitude) delete(longitude) delete(profile) ;---------------------------------------------- ; END LOOP ;---------------------------------------------- end do end ;---------------------------------------------- ; REFERENCES ;---------------------------------------------- ; [1] http://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/index.php ; [2] https://www.ncl.ucar.edu/Document/Functions/Built-in/dim_gbits.shtml ; [3] http://www.ncl.ucar.edu/Applications/Scripts/hdf4eos_5.ncl ; [4] http://www.ncl.ucar.edu/Applications/Scripts/calipso_3.ncl