Starting new all-in-one version of pouch_lag_hov

pouches/pouch_lag_hov.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; storm_lag_map.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; February 2014
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Draw a lagged composite for storms in a given lag of a wave
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+load "$CJS_NCL_LIB/lib.ibtracs.ncl"
+load "$CJS_NCL_LIB/lib.composites.ncl"
+load "$CJS_NCL_LIB/lib.cjs_graphics.ncl"
+load "/home/carl/projects/kelvin/lib/sub.match_storm_crest.ncl"
+
+
+begin ; main
+  print_clock( "Here we go!" )
+
+  ; These are some parameters that could be useful to have up top
+  if( .not.isvar( "level") ) then
+    level  = 850
+  end if
+  if( .not.isvar( "shdVarPre") ) then
+    shdVarPre  = "ow"
+  end if
+  if( .not.isvar("targetBasin") ) then
+    targetBasin = 0
+  end if
+  if( .not.isvar("stormMinLag") ) then
+    stormMinLag = 2
+  end if
+  if( .not.isvar("waveSpeed") ) then
+    waveSpeed = -0
+  else
+    waveSpeed = -waveSpeed
+  end if
+  print( waveSpeed )
+  if( .not.isvar("stormMaxLag") ) then
+    stormMaxLag = stormMinLag + 0.75
+  end if
+  if( .not.isvar("stormType") ) then
+    stormType = "genesis"
+  end if
+  shdVar = shdVarPre + level
+
+  plotType = "png"
+  plotDpi  = 300
+  fontHeightF = 0.02
+
+  waveName = "kelvin"
+  crestThresh = 1.0
+
+  txString = ibt_basin_name( targetBasin, True ) + "_" \
+      + stormType + "_" + stormMinLag + "_" + waveSpeed + "_" + shdVar
+  plotName = "figures/" + txString
+
+  stormMinLon = 000
+  stormMaxLon = 360
+  if( targetBasin.eq.-1 ) then
+    stormMinLat = -30
+    stormMaxLat = 30
+    minLat = -30
+    maxLat =  30
+  else
+    if( any( targetBasin.eq.(/ -3, 1, 4, 6, -4, -6, -9, 9 /) ) ) then
+      stormMinLat = -20
+      stormMaxLat = 00
+      minLat = -25
+      maxLat =  00
+      northHem = False
+    else
+      stormMinLat = 00
+      stormMaxLat = 20
+      minLat = -00
+      maxLat =  25
+      northHem = True
+    end if
+  end if
+
+;  lag = ispan( stormMinLag-3, stormMinLag, 1 )
+; lag = (/ stormMinLag-2, stormMinLag /)
+  lag = stormMinLag + ispan( -4, 0, 2 )
+  twoColumns = False
+  nTests  = 1
+  pThresh = 0.0
+
+  print_clock( txString )
+  scratchPath = "~/data/scratch/" + txString + ".nc"
+  newScratch = .not.isfilepresent( scratchPath ) 
+  if( newScratch ) then
+    scratchFile = addfile( scratchPath, "c" )
+  else
+    scratchFile = addfile( scratchPath, "r" )
+  end if
+
+
+  print_clock( "Getting the wave data for the storms" )
+  if( newScratch ) then
+    stormWaveData = match_storm_crest( stormType, targetBasin, \
+        stormMinLat, stormMaxLat, stormMinLon, stormMaxLon, waveName, \
+        crestThresh )
+    scratchFile->stormWaveData = stormWaveData
+  else
+    stormWaveData = scratchFile->stormWaveData
+  end if
+  crestInd = ind( ( stormWaveData@allCrestLag.ge.stormMinLag ) \
+      .and.( stormWaveData@allCrestLag.le.stormMaxLag ) )
+  stormInd = stormWaveData@allCrestStormInd(crestInd)
+  compLag = stormWaveData@allCrestLag(crestInd)
+  compLon = stormWaveData@lon(stormInd)
+  compLat = stormWaveData@lat(stormInd)
+  compDates = stormWaveData@time(stormInd) - compLag
+  compDates@units = stormWaveData@timeUnits
+
+  print( "Number of dates: " + dimsizes(compDates) )
+; print( cd_string( compDates, "" ) )
+  
+  centerLon = dim_median( compLon )
+  minLon = centerLon - 30
+  maxLon = centerLon + 30
+  print( minLon + " " + centerLon + " " + maxLon )
+
+  miss = new( 1, float )
+  miss = miss@_FillValue
+
+  if( newScratch ) then
+; if( False ) then
+    print_clock( "Compositing for u-wind" )
+    uVar = "u" + level
+    uPath = "/home/carl/data/merra/total/" + uVar + ".total.nc"
+    inData = read_xyt( uPath, uVar, miss, minLat, maxLat, miss, miss )
+    uData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( uData, True )
+
+
+    print_clock( "Compositing for v-wind" )
+    vVar = "v" + level
+    vPath = "/home/carl/data/merra/total/" + vVar + ".total.nc"
+    inData = read_xyt( vPath, vVar, miss, minLat, maxLat, miss, miss )
+    vData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( vData, True )
+
+
+    print_clock( "Compositing for OKW" )
+    owVar = "ow" + level
+    owPath = "/home/carl/data/merra/total/" + owVar + ".total.nc"
+    inData = read_xyt( owPath, owVar, miss, minLat, maxLat, miss, miss )
+    owData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( owData, True )
+
+    smthOwData = smth9_Wrap( owData, 0.5, 0.25, True )
+    printMinMax( smthOwData, True )
+
+    print_clock( "Compositing for vorticity" )
+    vortVar = "vort" + level
+    vortPath = "/home/carl/data/merra/total/" + vortVar + ".total.nc"
+    inData = read_xyt( vortPath, vortVar, miss, minLat, maxLat, miss, miss )
+    vortData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( vortData, True )
+
+  ; vortData = smth9_Wrap( vortData, 0.5, -0.25, True )
+  ; printMinMax( vortData, True )
+
+
+    print_clock( "Compositing for Kelvin wind" )
+    inPath = "/home/carl/data/merra/waves/" + uVar + ".waves.nc"
+    inData = read_xyt( inPath, waveName, miss, minLat, maxLat, miss, miss )
+    kelWindData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( kelWindData, True )
+
+
+    print_clock( "Compositing for Kelvin rain" )
+    inPath = "/home/carl/data/trmm/notc500/trmm3b42.waves.notc500.nc"
+    inData = read_xyt( inPath, waveName, miss, minLat, maxLat, miss, miss )
+    kelRainData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+        nTests, pThresh )
+    delete(inData)
+    printMinMax( kelRainData, True )
+
+
+    if( shdVarPre.eq."u" ) then
+      shdData = uData
+    else if( shdVarPre.eq."vort" ) then
+      shdData = vortData
+    else if( shdVarPre.eq."ow" ) then
+      shdData = owData
+    else
+      print_clock( "Compositing for shading" )
+      shdPath = "/home/carl/data/merra/total/" + shdVar + ".total.nc"
+      inData = read_xyt( shdPath, shdVar, miss, minLat, maxLat, miss, miss )
+      shdData = rotate_comp_xy_preread( compDates, compLon, lag, inData, \
+      nTests, pThresh )
+      delete(inData)
+    end if
+    end if
+    end if
+    printMinMax( shdData, True )
+
+    scratchFile->uData = uData
+    scratchFile->vData = vData
+    scratchFile->owData = owData
+    scratchFile->vortData = vortData
+    scratchFile->kelWindData = kelWindData
+    scratchFile->kelRainData = kelRainData
+    scratchFile->shdData = shdData
+
+  else
+    uData = scratchFile->uData
+    vData = scratchFile->vData
+    owData = scratchFile->owData
+    vortData = scratchFile->vortData
+    kelWindData = scratchFile->kelWindData
+    kelRainData = scratchFile->kelRainData
+    shdData = scratchFile->shdData
+
+  end if
+
+  uData = uData - waveSpeed
+  printMinMax( uData, True )
+
+  if( shdVarPre.eq."u" ) then
+    shdData = shdData - waveSpeed
+  end if
+
+
+  trfData = vData
+  if( northHem ) then
+    trfData = where( vortData.le.4, trfData@_FillValue, trfData )
+  else
+    trfData = where( vortData.ge.-4, trfData@_FillValue, trfData )
+  end if
+
+
+  ; Customize base plot
+  res = True
+  res@lbLabelBarOn = False
+
+  if( shdVarPre.eq."u" ) then
+    res@cnFillPalette = "colorbrewer_coldhot"
+    res@cnLevels = ispan(-5,5,1) * 2
+  end if
+  if( shdVarPre.eq."ow" ) then
+    res@cnFillPalette = "colorbrewer_coldhot"
+    res@cnLevels = ispan(-5,5,1) * 1.0
+  end if
+  if( shdVarPre.eq."epv" ) then
+    res@cnFillPalette = "colorbrewer_coldhot"
+    res@cnLevels = ispan(-5,5,1) * 0.05
+  end if
+  if( shdVarPre.eq."vort" ) then
+    res@cnFillPalette = "colorbrewer_coldhot"
+    res@cnLevels = ispan(-5,5,1) * 3.0
+  end if
+
+  ; ...set the bounds of a map plot
+  res@mpMinLatF    = minLat
+  res@mpMaxLatF    = maxLat
+  res@mpMinLonF    = minLon
+  res@mpMaxLonF    = maxLon
+
+  ; Customize overlaid contours
+  rainRes = True
+  rainRes@cnInfoLabelOn = False
+  rainRes@cnMonoLineColor   = False
+  rainRes@gsnContourNegLineDashPattern = 0
+  rainRes@cnLevels = (/-1,1/) * 0.1
+  rainRes@cnLineColors = where( rainRes@cnLevels.lt.0, "sienna1", "green3" )
+
+  windRes = True
+  windRes@cnInfoLabelOn = False
+  windRes@cnMonoLineColor   = False
+  windRes@gsnContourNegLineDashPattern = 0
+  windRes@cnLevels = (/-1,1/) * 0.5
+  windRes@cnLineColors = where( windRes@cnLevels.lt.0, "blue", "red" )
+
+  trfRes = True
+  trfRes@cnInfoLabelOn = False
+  trfRes@cnMonoLineColor   = True
+  trfRes@cnLevels = 0
+  trfRes@cnLineColor = "black"
+  trfRes@gsnContourZeroLineThicknessF = 5
+
+  ; Customize overlaid vectors
+  strRes = True
+  strRes@gsnLeftString     = ""
+  strRes@gsnRightString    = ""
+  strRes@gsnDraw = False
+  strRes@gsnFrame = False
+  strRes@stArrowLengthF = 0.01
+  strRes@stMinLineSpacingF = 0.020
+  strRes@gsnAddCyclic = True
+
+  txRes = True
+  txRes@txFontHeightF = fontHeightF * 2.0
+  txRes@txFontColor = "magenta"
+  txRes@txFontThicknessF = 2
+
+  ; Customize panel
+  panRes = True
+  panRes@txString       = txString \
+      + " (" + min(shdData@nDates) + " storms)"
+  panRes@gsnPanelBottom = 0.05
+  panRes@gsnPanelYWhiteSpacePercent = 10
+
+; if( twoColumns.or.( dimsizes(lag).eq.1 ) ) then
+; if( twoColumns ) then
+  if( True ) then
+    panRes@lbTitlePosition = "Right"
+    panRes@lbTitleDirection = "Across"
+    panRes@lbOrientation = "Horizontal"
+    panRes@pmLabelBarOrthogonalPosF = -.010
+    panRes@pmLabelBarHeightF        = 0.04
+  else
+    panRes@lbTitlePosition = "Top"
+    panRes@lbTitleDirection = "Across"
+    panRes@lbOrientation = "Vertical"
+  end if
+  
+  panRes@gsnPanelLabelBar   = True
+  panRes@lbTitleString      = shdData@units
+  panRes@lbLabelAutoStride  = True
+  panRes@lbTitleFontHeightF = fontHeightF * 0.8
+  panRes@lbLabelFontHeightF = fontHeightF * 0.8
+
+
+  print_clock( "Drawing the plot" )
+
+  ; ...allows png or gif to work
+  if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then
+    plotTypeLocal = "eps"
+  else
+    plotTypeLocal = plotType 
+  end if
+
+  ; ...open the workstation
+  wks = gsn_open_wks( plotTypeLocal, plotName )
+
+  plot = new( dimsizes(lag), graphic )
+
+  do iLag = 0, dimsizes(lag)-1
+
+;   res@gsnLeftString = inttochar(97+iLag) + ") Lag: " + lag(iLag) + " days"
+    res@gsnLeftString = inttochar(97+iLag) + ") "  \
+        + ( lag(iLag) - stormMinLag ) + " days from genesis"
+
+    if( twoColumns ) then
+      if( iLag.eq.dimsizes(lag)-1 ) then
+        res@tmXBLabelsOn    = True
+      else
+        if( iLag.eq.dimsizes(lag)/2-1 ) then
+          res@tmXBLabelsOn    = True
+        else
+          if( iLag.gt.dimsizes(lag)/2-1) then
+            res@tmYLLabelsOn = False
+          end if
+          res@tmXBLabelsOn    = False
+        end if
+      end if
+     
+    else
+      if( iLag.eq.dimsizes(lag)-1 ) then
+        res@tmXBLabelsOn    = True
+      else
+        res@tmXBLabelsOn    = False
+      end if
+    end if
+
+    plot(iLag) = cjs_draw_shaded_map( wks, shdData(iLag,:,:), res )
+    strmPlot = gsn_csm_streamline( wks, uData(iLag,:,:), vData(iLag,:,:),  \
+        strRes )
+    overlay( plot(iLag), strmPlot )
+
+    cjs_add_contours( wks, plot(iLag), kelRainData(iLag,:,:), rainRes )
+    cjs_add_contours( wks, plot(iLag), kelWindData(iLag,:,:), windRes )
+;   cjs_add_contours( wks, plot(iLag), trfData(iLag,:,:), trfRes )
+
+
+    waveLon = dim_median(compLon)  \
+        + ( lag(iLag) - avg((/stormMinLag,stormMaxLag/)) ) \
+        * waveSpeed * 86400 / 111000
+
+    if( lag(iLag).lt.stormMinLag ) then
+      genSymb = gsn_add_text( wks, plot(iLag), "~F37~m", waveLon, \
+          dim_median(compLat), txRes ) 
+    else
+      if( dim_median(compLat).gt.0 ) then
+        genSymb = gsn_add_text( wks, plot(iLag), "~F37~p", \
+            waveLon, dim_median(compLat), txRes ) 
+      else
+        genSymb = gsn_add_text( wks, plot(iLag), "~F37~s", \
+            waveLon, dim_median(compLat), txRes ) 
+      end if
+    end if
+  end do
+
+  nPlots = dimsizes(lag)
+  if( .not.twoColumns ) then
+    gsn_panel( wks, plot, (/ nPlots, 1 /), panRes )
+  else
+    plotInd = ispan(0,nPlots-1,1)
+    plotInd(0::2) = ispan(0,nPlots/2-1,1)
+    plotInd(1::2) = ispan(nPlots/2,nPlots-1,1)
+    gsn_panel( wks, plot(plotInd), (/ nPlots/2, 2 /), panRes )
+  end if
+
+
+
+  print_clock( "Convert the image, if necessary" )
+  delete(wks)
+  if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
+    system( "convert -trim -border 5x5 -bordercolor white " \
+        + "+repage -density " + plotDpi + " " \\
+        + plotName + ".eps " + plotName + "." + plotType )
+    if( .not.isStrSubset( plotType, "e" ) ) then
+      system( "rm -f " + plotName + ".eps" )
+    end if
+  end if
+
+
+  print_clock( "Thank you, come again." )
+end; main
+