; ; ; ============================================================================ ; $RCSfile$ ; $Source$ ; $Revision$ ; $Date$ ; $Author$ ; ============================================================================ ; ; Default Variable Initialization - DO NOT change! ========================== ; ============================================================================ ; External libraries 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"; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"; print(""); print("=============================================================================="); ; User specified configuration parameters ==================================== ; THREDDS AVHRR url ;satUrl = "http://tashtego.marine.rutgers.edu:8080/thredds/dodsC/cool/avhrr/bigbight/2010"; satUrl = "http://tds.maracoos.org/thredds/dodsC/SST-Agg.nc"; ; HF Radar url hfrUrl = "http://tashtego.marine.rutgers.edu:8080/thredds/dodsC/cool/codar/totals/macoora6km_clone"; ;hfrUrl = "http://tashtego.marine.rutgers.edu:8080/thredds/dodsC/cool/codar/totals/macoora6km"; ; Image destination defaultImgDir = "/www/home/kerfoot/public_html/beta/codar/img/maracoos/25hroverlays"; ; Multiplier to convert hours to seconds HOURS2SECONDS = 3600; ; Multiplier to convert cm/s to m/s CMS2MS = 0.01; ; Number of hours in a day (converting codar timestamps) DAYS2HOURS = 24; ; Domain plotting bounds - comment out to select default values from the ; dataset defaultSLat = 35.0; defaultNLat = 42.25 defaultWLon = -77.0; defaultELon = -68.0; ; Averaging/compositing window defaultWindow = 25.0; ; SST compositing method defaultCompMethod = "warmest"; ; Default workstation type (available: x11, ps, eps, epsi, pdf, png) defaultWksType = "ps"; ; Lower minimum value for color scaling SST imagery defaultMinTemp = 1.0; ; ============================================================================ ; SET UP all variables ======================================================= ; Set a default image destination if none was specified if (.not. isvar("imgDir") .and. .not. isvar("imgRoot")) then imgDir = ""; wksType = "x11"; else wksType = defaultWksType; end if print(wksType); ; If the variable window does NOT exist, use the value stored in defaultWindow if (.not. isvar("window")) then halfWindow = (defaultWindow/2)/DAYS2HOURS; window = defaultWindow; else if (isinteger(window)) then ; convert to float if string tmpWindow = window; delete(window); window = int2flt(tmpWindow); end if halfWindow = (window/2)/DAYS2HOURS; end if ;print(window); ;print(halfWindow); ;exit; ; Set the cloud mask if (.not. isvar("maskLevel")) then maskLevel = 0; end if ; Set up the color limits if specified if (isvar("cLim")) then print(cLim); exit; end if ; Set up NSEW bounds: ismissing(PARAM) returns true if the variable was not ; set explicitly with a float. If this is the case, use the default value ; from the dataset. if (.not. isvar("nLat")) then nLat = defaultNLat; end if if (.not. isvar("sLat")) then sLat = defaultSLat; end if if (.not. isvar("wLon")) then wLon = defaultWLon; end if if (.not. isvar("eLon")) then eLon = defaultELon; end if print("Plotting bounds:"); print("N: " + nLat); print("S: " + sLat); print("E: " + eLon); print("W: " + wLon); ; ============================================================================ ; Process HF Radar data ====================================================== print("Retrieving HF Radar data..."); ; Open up the url for reading hfrFid = addfile(hfrUrl, "r"); ;printVarSummary(hfrFid); if (.not. isvar("ts")) then print("Retrieving last " + window + " hours of data."); ; Get the number of records tDims = filevardimsizes(hfrFid, "time"); ; Select the last timestamp hft1 = hfrFid->time(tDims - 1); hft0 = hfrFid->time(tDims - floattointeger(window)); else ; Assign attributes to the timestamp for date conversions ts@units = "seconds since 1970-01-01 00:00:00"; ts@calendar = "gregorian"; print(ts); hfrTs = ut_convert(ts, "days since 2001-01-01 00:00:00"); ; Selected timeslices in codar date units hft0 = hfrTs - halfWindow; hft1 = hfrTs + halfWindow; hft0@units = "days since 2001-01-01 00:00:00"; hft0@calendar = "gregorian"; hft1@units = "days since 2001-01-01 00:00:00"; hft1@calendar = "gregorian"; end if ;print(hft0); ;print(hft1); ; Select the hourly files uvTs = hfrFid->time({hft0:hft1}); ;print(uvTs); hfCenterTs = dim_avg_Wrap(uvTs); hfCenterTsPieces = ut_calendar(hfCenterTs,0); ; Create center timestamp for filename (yyyymmddTHHMMSS) tsFilename = sprinti("%0.4i", floattointeger(hfCenterTsPieces(:,0))) + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,1))) + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,2))) + \ "T" + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,3))) + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,4))) + \ "00"; ;print(tsFilename); hfCenterTsLabel = sprinti("%0.4i", floattointeger(hfCenterTsPieces(:,0))) + \ "-" + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,1))) + \ "-" + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,2))) + \ " " + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,3))) + \ ":" + \ sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,4))) + \ ":" + \ "00" + \ " UTC"; ; Format the timestamps uvTsPieces = ut_calendar(uvTs,0); uvYear = floattointeger(uvTsPieces(:,0)); uvMonth = floattointeger(uvTsPieces(:,1)); uvDay = floattointeger(uvTsPieces(:,2)); uvHour = floattointeger(uvTsPieces(:,3)); uvMinute = floattointeger(uvTsPieces(:,4)); uvSeconds = floattointeger(uvTsPieces(:,5)); uvDates = sprinti("%0.4i", uvYear) + \ "-" + \ sprinti("%0.2i", uvMonth) + \ "-" + \ sprinti("%0.2i", uvDay) + \ " " + \ sprinti("%0.2i", uvHour) + \ ":" + \ sprinti("%0.2i", uvMinute) + \ ":" + \ sprinti("%0.2i", uvSeconds) + \ " UTC"; print(uvDates); ; Select the time slice u = hfrFid->u({hft0:hft1},{sLat:nLat},{wLon:eLon}); ; values are shorts ;printVarSummary(u); ;exit; ; Convert from cm/s to m/s ui = u * CMS2MS; ; Copy all variable metadata from u to ui copy_VarMeta(u, ui); ui@units = "m/s"; ;printVarSummary(ui); uMax = max(abs(ui)); ;print(uMax); ; Select the most recent time slice v = hfrFid->v({hft0:hft1},{sLat:nLat},{wLon:eLon}); ; values are shorts ;printVarSummary(v); ; Apply scale factor to convert to m/s (short -> float) - HFRNet vi = v * CMS2MS; ; Copy all variable metadata from v to vi copy_VarMeta(v, vi); vi@units = "m/s"; vMax = max(abs(vi)); ;print(vMax); ; Average the ui and vi grids over the time dimension print("Averaging U through the time dimension..."); U = dim_avg_n_Wrap(ui,0); print("Done."); print("Averaging V through the time dimension..."); V = dim_avg_n_Wrap(vi,0); print("Done."); ; Display the min and max values of u and v print("------------------------------------------------------------------------------"); print("Max u velocity: " + uMax + " " + ui@units); print("Max v velocity: " + vMax + " " + vi@units); print("------------------------------------------------------------------------------"); ; END HF Radar processing ==================================================== ; Process SST first ========================================================== print("Retrieving SST data..."); ; Open up the satUrl for reading satFid = addfile(satUrl, "r"); ;printVarSummary(satFid); ;exit; ; Set the compositing method compMethod = defaultCompMethod; if (isvar("method")) then compMethod = method; end if ; Convert codar dates (hft0 and hft1) to unix time for selecting satellite ; passes t0 = doubletointeger(ut_convert(hft0, "seconds since 1970-01-01 00:00:00")); t1 = doubletointeger(ut_convert(hft1, "seconds since 1970-01-01 00:00:00")); ; Add attributes t0@units = "seconds since 1970-01-01 00:00:00"; t0@calendar = "gregorian"; t1@units = "seconds since 1970-01-01 00:00:00"; t1@calendar = "gregorian"; ;print(t0); ;print(t1); passTs = satFid->time({t0:t1}); ;print(passTs); centerTime = dim_avg_Wrap(passTs); ;exit; centerTsPieces = ut_calendar(centerTime,0); passCenterTs = sprinti("%0.4i", floattointeger(centerTsPieces(:,0))) + \ "-" + \ sprinti("%0.2i", floattointeger(centerTsPieces(:,1))) + \ "-" + \ sprinti("%0.2i", floattointeger(centerTsPieces(:,2))) + \ " " + \ sprinti("%0.2i", floattointeger(centerTsPieces(:,3))) + \ ":" + \ sprinti("%0.2i", floattointeger(centerTsPieces(:,4))) + \ ":" + \ sprinti("%0.2i", floattointeger(centerTsPieces(:,5))) + \ " UTC"; print("SST Composite Center Time: " + \ passCenterTs + \ " (+/- " + \ floattoint(window/2) + \ " hours)"); ; Format the pass timestamps tsPieces = ut_calendar(passTs,0); tsYear = floattointeger(tsPieces(:,0)); tsMonth = floattointeger(tsPieces(:,1)); tsDay = floattointeger(tsPieces(:,2)); tsHour = floattointeger(tsPieces(:,3)); tsMinute = floattointeger(tsPieces(:,4)); tsSeconds = floattointeger(tsPieces(:,5)); passTimes = sprinti("%0.4i", tsYear) + \ "-" + \ sprinti("%0.2i", tsMonth) + \ "-" + \ sprinti("%0.2i", tsDay) + \ " " + \ sprinti("%0.2i", tsHour) + \ ":" + \ sprinti("%0.2i", tsMinute) + \ ":" + \ sprinti("%0.2i", tsSeconds) + \ " UTC"; ; Display the formatted pass times print(passTimes); ;exit; ; Select the SST from the file for the specified time mcsst = satFid->mcsst({t0:t1},{sLat:nLat},{wLon:eLon}); ;printVarSummary(mcsst); ;exit; ; Get and display the min/max SST values before applying cloud mask minSST = min(mcsst); maxSST = max(mcsst); print("------------------------------------------------------------------------------"); print("Original Temperature Minimum: " + minSST + " " + mcsst@units); print("Original Temperature Maximum: " + maxSST + " " + mcsst@units); print("------------------------------------------------------------------------------"); ; Apply cloud mask if cloudMask > 0 if (maskLevel .gt. 0) then print("Masking values < " + maskLevel + "..."); cloudMask = satFid->cloud_land_mask({t0:t1},{sLat:nLat},{wLon:eLon}); ; printVarSummary(cloudMask); mcsst = mask(mcsst, (cloudMask.lt.maskLevel), True); minSST = min(mcsst); maxSST = max(mcsst); print("------------------------------------------------------------------------------"); print("Masked Temperature Minimum : " + minSST + " " + mcsst@units); print("Masked Temperature Maximum : " + maxSST + " " + mcsst@units); print("------------------------------------------------------------------------------"); end if ;printVarSummary(mcsst); ; Warmest pixel composite if more than one pass was selected mcsstDims = dimsizes(mcsst); if (mcsstDims(0) .gt. 0) then if (compMethod .eq. "average") then print("Compositing method: Pixel Average..."); sst = dim_avg_n_Wrap(mcsst,0); else if (compMethod .eq. "warmest") then print("Compositing method: Warmest Pixel..."); sst = dim_max_n(mcsst,0); copy_VarMeta(mcsst(lat | :, lon | :, time | :),sst); end if end if else sst = mcsst; end if ; Don't need the original mcsst variable anymore delete(mcsst); ;printVarSummary(sst); ;exit; ; Unless the user has specified a min/max temperature for the colorbar, we're ; going to autoscale; however, limit the lower bounds of the autoscale to the ; value stored in sstAvg = avg(sst); ;print(sstAvg); sstDev = stddev(sst); ;print(sstDev); ;exit; if (.not. isvar("minTemp")) then minTemp = sstAvg - (2*sstDev); Take 2 standard deviations print("Autoscaling minimum temperature: " + minTemp); if (minTemp .lt. defaultMinTemp) then print("Minimum temperature (" + \ defaultMinTemp + \ ") exceeded."); print("Setting default minimum temperature: " + \ defaultMinTemp); minTemp = floattoint(defaultMinTemp); end if end if if (.not. isvar("maxTemp")) then maxTemp = floattoint(sstAvg + (2*sstDev)); print("Autoscaling maximum temperature: " + maxTemp); end if ;landMask = satFid->cloud_land_mask(0,{sLat:nLat},{wLon:eLon}); ;sst = mask(sst, (landMask.eq.99), True); ; END satellite SST processing =============================================== ; Set up filename imgName = "maracoos-overlay_" + \ tsFilename + \ "_" + \ window + \ "hroverlay"; img = imgName; ; Create the image destination. There are 2 possibilities: ; if imgDir exists, it is the absolute path to the write directory ; if imgDir does not exist, but imgRoot does, it is the root directory to ; we will write (and create if necessary) the imgRoot/yyyy/mm directory to ; write the images if (.not. isvar("imgDir")) then imgDir = ""; if (isvar("imgRoot")) then yyyy = sprinti("%0.4i", floattointeger(hfCenterTsPieces(:,0))); month = sprinti("%0.2i", floattointeger(hfCenterTsPieces(:,1))); imgDir = imgRoot + "/" + yyyy + "/" + month; success = systemfunc("/bin/mkdir -pm 755 " + imgDir + " && echo 0"); ;print(success); if (ismissing(success)) then print("Can't create output directory: " + imgDir); exit; end if end if end if img = imgDir + "/" + imgName; print("Writing " + img); ; Open up a work station and create the image wks = gsn_open_wks(wksType, img); ; Set up SST drawing resources sstRes = True; ; High-level gsn sstResources =================================================== sstRes@gsnMaximize = True; ; Maximize image size on paper sstRes@gsnAddCyclic = False; ; ============================================================================ ; Display best tick mark labeling sstRes@pmTickMarkDisplayMode = "Always"; ; Map settings =============================================================== sstRes@mpFillDrawOrder = "PostDraw"; draw the land mask after plotting contours sstRes@mpLandFillColor = "gray"; Map fill color is gray sstRes@mpDataBaseVersion = "HighRes"; High resolution map ;sstRes@mpDataSetName = "Earth..2"; ;sstRes@mpOutlineBoundarySets = "USStates"; State lines sstRes@mpProjection = "Mercator"; Mercator projection ; Map domains set by Latititude and Longitude bounds sstRes@mpLimitMode = "LatLon"; sstRes@mpMinLatF = sLat; sstRes@mpMaxLatF = nLat; sstRes@mpMinLonF = wLon; sstRes@mpMaxLonF = eLon; ; Turn on grid sstRes@mpGridAndLimbOn = True; ; Set gridline intervals sstRes@mpGridSpacingF = 0.5; ; 1/2 degree intervals sstRes@mpGridLineDashPattern = 2; ; dashed line ; Turn on map grid sstRes@mpGridAndLimbOn = True; ; Set gridline intervals sstRes@mpGridSpacingF = 0.5; ; 1/2 degree intervals sstRes@mpGridLineDashPattern = 2; ; dashed line sstRes@mpGridLineThicknessF = 0.5; ; ============================================================================ ; Set up colormap ============================================================ ;gsn_define_colormap(wks, "rainbow+gray"); gsn_define_colormap(wks, "sst"); sstRes@gsnSpreadColors = True; Use the entire colormap sstRes@gsnSpreadColorStart = 16; sstRes@gsnSpreadColorEnd = -4; ; ============================================================================ ; Temperature limits ========================================================= sstRes@cnLevelSelectionMode = "ManualLevels"; sstRes@cnMinLevelValF = minTemp; sstRes@cnMaxLevelValF = maxTemp; ; ============================================================================ ; LabelBar/Colorbar sstResources (lb) =========================================== ; Set to true to make sure colorbar labels don't overlap ;sstRes@lbLabelAutoStride = True; sstRes@lbLabelStride = 5; ; No lines on colorbar intervals sstRes@lbBoxLinesOn = False; ; Place the label bar on the right side of the plot sstRes@lbOrientation = "Vertical"; ; faster drawing of label bar sstRes@lbRasterFillOn = True; sstRes@lbLabelFontHeightF = .015; Cut fontsize in 1/2 from default sstRes@lbTitleString = "Sea Surface Temperature (" + sst@units + ")"; sstRes@lbTitleFontHeightF = 0.015; sstRes@lbTitlePosition = "Right"; sstRes@lbTitleDirection = "Across"; sstRes@lbTitleAngleF = 90.0; sstRes@lbLeftMarginF = -0.25; sstRes@lbRightMarginF = 0.35; ; ============================================================================ ; Contouring sstResources ======================================================= sstRes@cnFillOn = True; color-filled contours sstRes@cnFillMode = "RasterFill"; Faster plotting and no plotting on land mask sstRes@cnLevelSpacingF = 0.2; ; Spacing between contour levels sstRes@cnMonoLevelFlag = True; Allow me to turn off lines and labels sstRes@cnLevelFlag = "NoLine"; turn off lines and labels ;sstRes@cnMissingValFillPattern = -1; ;sstRes@cnMissingValFillColor = -1; ; ============================================================================ ; Plot annotation sstResources ================================================== sstRes@tiMainString = ""; sstRes@gsnLeftString = ""; sstRes@gsnRightString = ""; sstRes@gsnCenterStringFontHeightF = 0.015; sstRes@gsnCenterString = "MARACOOS SST/Current Field Overlay: " + \ hfCenterTsLabel + \ " (+/- " + \ floattoint(window/2) + \ " hours)" ; ; ============================================================================ ; Tickmarks ================================================================== sstRes@tmXBLabelFontHeightF = 0.012; sstRes@tmYLLabelFontHeightF = 0.012; ; ============================================================================ ; Do not draw contour immediately (wait for the codar plot) sstRes@gsnDraw = False; ; Do not advance the frame (wait for the codar plot) sstRes@gsnFrame = False; ; Draw sstPlot = gsn_csm_contour_map(wks, sst, sstRes); ;; Load and plot state boundaries ;pres = True; ;pres@tfPolyDrawOrder = "Draw"; ;slFile = "/www/home/kerfoot/public_html/shapefiles/usgs/state_boundaries/sblines.dat"; ;; Number of rows in the file ;rows = numAsciiRow(slFile); ;;; Read the data into a 2-column array ;gps = asciiread(slFile, (/rows, 2/), "float"); ;;lons = gps(:,0); ;;lats = gps(:,1); ;;;print(lons); ;;;print(lats); ;;;exit; ;pRes = True; ;pRes@tfPolyDrawOrder = "PostDraw"; ;boundaries = gsn_add_polyline(wks, sstPlot, gps(:,0), gps(:,1), pRes); ; Vector map plot resource settings uvres = True; uvres@gsnDraw = False; uvres@gsnFrame = False; ; Vector styling uvres@vcGlyphStyle = "CurlyVector"; uvres@vcMinDistanceF = 0.005; spacing between vectors (larger value = less vectors) uvres@vcRefLengthF = 0.08; vector scaling factor ; (NDC units: http://www.ncl.ucar.edu/Document/Graphics/ndc.shtml) ; larger value = longer vectors uvres@vcRefMagnitudeF = 0.5 ; Length (in uv units) of reference vector. Any ; vector of this size will be rendered at the ; length of vcRefLength uvres@vcLineArrowThicknessF = 0.5; ; Y-dir position of reference vector uvres@vcRefAnnoOrthogonalPosF = -1.0; ; X-dir position of reference vector uvres@vcRefAnnoParallelPosF = 0.18; uvres@vcMagnitudeFormat = "@0.2f"; ; '%02.2f' uvres@vcRefAnnoString1 = "$VMG$ " + u@units; uvres@vcRefAnnoString2 = ""; uvres@vcLineArrowColor = 1; ; Vector color ; Plot annotation resources uvres@tiMainString = ""; uvres@gsnLeftString = ""; uvres@gsnRightString = ""; uvres@gsnCenterString = ""; print("Plotting surface current field..."); ; Plot the vectors vPlot = gsn_csm_vector(wks, U, V, uvres); print("Done."); ; overlay the plots overlay(sstPlot, vPlot); maximize_output(wks, True); ; If an image root (imgRoot) has been specified, the script is being run to ; generate imagery in the imgRoot/yyyy/mm directories (ie: cron). Convert ; created .ps files to png and then remove them if (isvar("imgRoot")) then system("/home/kerfoot/bin/bash/img/codarPs2png.sh " + \ imgDir + \ "/*.ps"); system("rm " + \ imgDir + \ "/*.ps"); end if