Commit 3e09a29e authored by Carl Schreck's avatar Carl Schreck
Browse files

Automated Nightly Commit - Fri Oct 17 00:02:35 EDT 2014

parent 3a19b1c9
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; map.ncl
; Carl Schreck (carl@cicsnc.org)
; October 2011
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Draw a quick-look map of data in a file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Carl Schreck (cjschrec@ncsu.edu)
; September 2014
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Compare maps for a day
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
......@@ -14,134 +14,109 @@ load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.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/lib.cjs_graphics.ncl"
load "$CJS_NCL_LIB/lib.array.ncl"
load "$CJS_NCL_LIB/lib.percentiles.ncl"
load "$CJS_NCL_LIB/print_clock.ncl"
load "/home/carl/projects/olr_cdr/sub.calc_meanvar.ncl"
begin
begin ; main
print_clock( "Here we go!" )
; These are some parameters that could be useful to have up top
sensor = "avhrr_25"
; sensor = "hirs_25"
timeUnits = "days since 1800-01-01 00:00:00"
; plotTime = cd_inv_calendar( 1999, 09, 15, 12, 0, 0, timeUnits, 0 )
; plotTime = cd_inv_calendar( 2011, 11, 25, 12, 0, 0, timeUnits, 0 )
plotTime = cd_inv_calendar( 2008, 02, 25, 12, 0, 0, timeUnits, 0 )
plotType = "png"
plotName = "figures/map_" + sensor + cd_string( plotTime, "_%Y-%N" )
plotDpi = 200
fontHeightF = 0.02
; shdLevels = ispan(-7,7,1) * 10
shdLevels = ispan(100,330,10)
; shdLevels = (/ 0.01, 0.1, 0.2, 0.4, 1.0, 2.0, 4.0 /)
; shdLevels = -999
minLon = 000
maxLon = 360
minLat = -90
maxLat = 90
centerLon = 180
print_clock( "Reading the data" )
var = "olr"
path = "/home/carl/data/olr/cdr/" + sensor + ".total.nc"
; var = "rain"
; path = "/home/carl/data/trmm/trmm3b42.nc"
; var = "brtmp"
; path = "/home/carl/data/claus/claus.1deg.fill.nc"
plotDpi = 200
inFile = addfile( path, "r" )
plotTime = cd_convert( plotTime, inFile->time@units )
timeUnits = plotTime@units
data = inFile->$var$({plotTime},{minLat:maxLat},:)
plotTime = inFile->time({plotTime})
totalNotAnom = True
if( .not.isvar("varName") ) then
varName = "olr"
end if
if( totalNotAnom ) then
fileName = "olr.total.nc"
plotName = "figures/map.total"
else
fileName = "olr.anom.nc"
plotName = "figures/map.anom"
end if
minLat = -60
maxLat = 60
printMinMax( data, True )
timeUnits = "days since 1800-01-01 00:00:00"
plotTime = cd_inv_calendar( 2011, 11, 27, 00, 0, 0, timeUnits, 0 )
plotTime = cd_inv_calendar( 1987, 05, 15, 00, 0, 0, timeUnits, 0 )
plotName = plotName + cd_string( plotTime, ".%Y-%N-%D" )
dataList = NewList( "lifo" )
minTotal = -999.
minTotal@_FillValue = minTotal
maxTotal = minTotal
sensor = "avhrr"
print_clock( "Calculating " + sensor )
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
inFile = addfile( inPath, "r" )
avhrrData := inFile->$varName$({plotTime},:,:)
avhrrData@long_name = sensor
minTotal = min( (/ minTotal, percent_to_value( avhrrData({minLat:maxLat},:), 01 ) /) )
maxTotal = max( (/ maxTotal, percent_to_value( avhrrData({minLat:maxLat},:), 99 ) /) )
ListPush( dataList, avhrrData )
sensor = "hirs"
print_clock( "Calculating " + sensor )
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
inFile = addfile( inPath, "r" )
hirsData := inFile->$varName$({plotTime},:,:)
hirsData@long_name = sensor
minTotal = min( (/ minTotal, percent_to_value( hirsData({minLat:maxLat},:), 01 ) /) )
maxTotal = max( (/ maxTotal, percent_to_value( hirsData({minLat:maxLat},:), 99 ) /) )
ListPush( dataList, hirsData )
nData = ListCount( dataList )
print( dataList )
; Customize base plot
res = True
res@gsnDraw = False
res@gsnFrame = False
; res@gsnLeftString = ""
res@gsnRightString = cd_string(plotTime,"")
; ...standardize the font sizes
res@cnInfoLabelFontHeightF = fontHeightF
res@gsnLeftStringFontHeightF = fontHeightF
res@gsnRightStringFontHeightF = fontHeightF
res@lbTitleFontHeightF = fontHeightF
res@lbLabelFontHeightF = fontHeightF
res@tiMainFontHeightF = fontHeightF * 1.3
res@tiYAxisFontHeightF = fontHeightF
res@tiXAxisFontHeightF = fontHeightF
res@tmXBLabelFontHeightF = fontHeightF
res@tmYLLabelFontHeightF = fontHeightF
; ...shading resources
res@cnFillOn = True
res@cnFillMode = "CellFill"
res@cnRasterSmoothingOn = True
res@cnMonoFillColor = False
res@cnInfoLabelOn = False
res@cnLineLabelsOn = False
res@cnLinesOn = False
; res@cnLineColor = "gray"
res@cnMissingValFillPattern = "SolidFill"
; res@cnMissingValFillColor = "gray"
; res@cnFillPalette = "colorbrewer_wet"
; ...set the contour levels
if( shdLevels(0).eq.-999 ) then
res@cnMaxLevelCount = 16
res@lbOrientation = "Horizontal"
res@cnLevelSelectionMode = "ManualLevels"
res@gsnMajorLonSpacing = 45
res@gsnMinorLonSpacing = 15
res@pmLabelBarOrthogonalPosF = 0.3
res@gsnRightString = cd_string( plotTime, "%d %c %Y" )
totalRes = res
if( totalNotAnom ) then
totalRes@cnFillPalette = "CBR_wet"
totalRes@gsnSpreadColorStart = 10
totalRes@gsnSpreadColorEnd = 0
mnmxintvl = nice_mnmxintvl( minTotal, maxTotal, 10, False )
totalRes@cnMinLevelValF = mnmxintvl(0)
totalRes@cnMaxLevelValF = mnmxintvl(1)
totalRes@cnLevelSpacingF = mnmxintvl(2)
else
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = shdLevels
totalRes@cnFillPalette = "CBR_drywet"
totalRes@gsnSpreadColorStart = 10
totalRes@gsnSpreadColorEnd = 0
maxTotal = max( abs( (/ minTotal, maxTotal /) ) )
mnmxintvl = nice_mnmxintvl( -maxTotal, maxTotal, 10, False )
totalRes@cnMinLevelValF = mnmxintvl(0)
totalRes@cnMaxLevelValF = mnmxintvl(1)
totalRes@cnLevelSpacingF = mnmxintvl(2)
end if
; ...labelbar resources
res@lbLabelBarOn = True
res@lbOrientation = "Horizontal"
res@lbTitlePosition = "Right"
res@lbTitleDirection = "Across"
res@lbTitleString = data@units
res@pmLabelBarOrthogonalPosF = fontHeightF * 10
res@pmLabelBarHeightF = fontHeightF * 2
; ...map gridlines
res@mpGridLatSpacingF = 10
res@mpGridLonSpacingF = 10
res@mpGridMaskMode = "MaskLand"
; res@mpGridLineColor = "gray50"
; ...set the bounds of a map plot
res@mpMinLatF = minLat
res@mpMaxLatF = maxLat
res@mpMinLonF = minLon
res@mpMaxLonF = maxLon
res@mpCenterLonF = centerLon
res@gsnAddCyclic = True
; ...set map resources
res@mpFillOn = False
res@mpGeophysicalLineThicknessF = 2
; res@mpGeophysicalLineColor = "gray50"
; res@mpNationalLineColor = "gray50"
; res@mpUSStateLineColor = "gray50"
res@mpOutlineBoundarySets = "AllBoundaries"
; res@mpDataBaseVersion = "MediumRes"
; res@mpDataSetName = "Earth..1"
diffRes = res
diffRes@cnFillPalette = "CBR_coldhot"
; Customize panel
panRes = True
panRes@gsnPanelRowSpec = False
panRes@gsnPanelCenter = True
panRes@gsnPanelYWhiteSpacePercent = 5
print_clock( "Drawing the plot" )
; ...allows png or gif to work
if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then
if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
plotTypeLocal = "eps"
else
plotTypeLocal = plotType
......@@ -149,20 +124,60 @@ begin
; ...open the workstation
wks = gsn_open_wks( plotTypeLocal, plotName )
plot = gsn_csm_contour_map( wks, data, res )
draw(plot)
frame(wks)
; gsn_define_colormap( wks, "default" )
nPlots = nData * ( nData + 1 ) / 2
plots = new( nPlots, graphic )
meanNotDiff = new( nPlots, logical )
plotCounter = 0
do xCounter = 0, nData-1
xData := dataList[xCounter]
totalRes@gsnLeftString = inttochar(97+plotCounter) + ") " \
+ xData@long_name
print_clock( totalRes@gsnLeftString )
plots(plotCounter) = cjs_draw_shaded_map( wks, xData, totalRes )
plotCounter = plotCounter + 1
do yCounter = 0, xCounter-1
yData := dataList[yCounter]
yRegrid := copy_Gridding( xData, yData )
xRegrid := copy_Gridding( yData, xData )
diffData := yData - xRegrid
copy_VarMeta( yData, diffData )
; diffData := yRegrid - xData
; copy_VarMeta( yRegrid, diffData )
diffRes@gsnLeftString = inttochar(97+plotCounter) + ") " \
+ yData@long_name + " minus " + xData@long_name
print_clock( diffRes@gsnLeftString )
currMax = percent_to_value( diffData, 99 )
mnmxintvl := nice_mnmxintvl( -currMax, currMax, 10, False )
diffRes@cnMinLevelValF = mnmxintvl(0)
diffRes@cnMaxLevelValF = mnmxintvl(1)
diffRes@cnLevelSpacingF = mnmxintvl(2)
plots(plotCounter) = cjs_draw_shaded_map( wks, diffData, diffRes )
plotCounter = plotCounter + 1
end do
end do
gsn_panel( wks, plots, (/ nPlots, 1 /), panRes )
print_clock( "Convert the image, if necessary" )
delete(wks)
if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then
if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
system( "convert -trim -border 5x5 -bordercolor white " \
+ "+repage -density " + plotDpi + " " \\
+ plotName + ".eps " + plotName + "." + plotType )
system( "rm -f " + plotName + ".eps" )
if( .not.isStrSubset( plotType, "e" ) ) then
system( "rm -f " + plotName + ".eps" )
end if
end if
print_clock( "Thank you, come again." )
end
end; main
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment