Commit 5bc0a9d4 authored by Carl Schreck's avatar Carl Schreck
Browse files

Automated Nightly Commit - Sat Oct 4 00:01:31 EDT 2014

parent cd199657
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; meanvar.ncl
; Carl Schreck (cjschrec@ncsu.edu)
; September 2014
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Compare maps of variance or mean
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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/lib.cjs_graphics.ncl"
load "$CJS_NCL_LIB/lib.array.ncl"
load "$CJS_NCL_LIB/print_clock.ncl"
load "/home/carl/projects/olr_cdr/sub.calc_meanvar.ncl"
begin ; main
print_clock( "Here we go!" )
; These are some parameters that could be useful to have up top
plotType = "png"
plotDpi = 200
meanNotVar = True
varName = "olr"
if( meanNotVar ) then
fileName = "olr.total.nc"
plotName = "mean"
else
fileName = "olr.anom.nc"
; fileName = "olr.std.waves.nc"
plotName = "variance"
end if
; plotName = varName
dataList = NewList( "lifo" )
minTotal = -999.
minTotal@_FillValue = minTotal
maxTotal = minTotal
sensor = "20th"
minYear = 1911
maxYear = 1940
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
ancientData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
ancientData@long_name = long_name
minTotal = min( (/ minTotal, min( ancientData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( ancientData({-60:60},:) ) /) )
ListPush( dataList, ancientData )
sensor = "20th"
minYear = 1981
maxYear = 2010
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
modernData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
modernData@long_name = long_name
minTotal = min( (/ minTotal, min( modernData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( modernData({-60:60},:) ) /) )
ListPush( dataList, modernData )
sensor = "avhrr"
minYear = 1981
maxYear = 2010
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
avhrrData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
avhrrData@long_name = long_name
minTotal = min( (/ minTotal, min( avhrrData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( avhrrData({-60:60},:) ) /) )
ListPush( dataList, avhrrData )
sensor = "hirs"
minYear = 1981
maxYear = 2010
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
hirsData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
hirsData@long_name = long_name
minTotal = min( (/ minTotal, min( hirsData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( hirsData({-60:60},:) ) /) )
ListPush( dataList, hirsData )
nData = ListCount( dataList )
print( dataList )
; Customize base plot
totalRes = True
totalRes@cnFillPalette = "CBR_wet"
totalRes@cnLevelSelectionMode = "ManualLevels"
if( meanNotVar ) then
totalRes@gsnSpreadColorStart = 10
totalRes@gsnSpreadColorEnd = 0
end if
mnmxintvl = nice_mnmxintvl( minTotal, maxTotal, 10, False )
totalRes@cnMinLevelValF = mnmxintvl(0)
totalRes@cnMaxLevelValF = mnmxintvl(1)
totalRes@cnLevelSpacingF = mnmxintvl(2)
totalRes@pmLabelBarOrthogonalPosF = 0.30
diffRes = True
diffRes@cnFillPalette = "CBR_drywet"
diffRes@cnLevelSelectionMode = "ManualLevels"
diffRes@pmLabelBarOrthogonalPosF = 0.30
; Customize panel
panRes = True
panRes@gsnPanelRowSpec = True
panRes@gsnPanelCenter = False
print_clock( "Drawing the plot" )
; ...allows png or gif to work
if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
plotTypeLocal = "eps"
else
plotTypeLocal = plotType
end if
; ...open the workstation
wks = gsn_open_wks( plotTypeLocal, plotName )
; 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 )
diffRes@gsnLeftString = inttochar(97+plotCounter) + ") " \
+ yData@long_name + " minus " + xData@long_name
print_clock( diffRes@gsnLeftString )
currMax = 0.9 * max( abs(diffData({-60:60},:)) )
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, ispan( 1, nData, 1 ), panRes )
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
......@@ -27,46 +27,53 @@ begin ; main
plotType = "png"
plotDpi = 200
meanNotVar = True
varName = "olr"
meanNotVar = False
if( .not.isvar("varName") ) then
varName = "olr"
end if
if( meanNotVar ) then
fileName = "olr.total.nc"
plotName = "mean"
plotName = "figures/map.mean"
else
fileName = "olr.anom.nc"
; fileName = "olr.std.waves.nc"
plotName = "variance"
if( varName.eq."olr" ) then
fileName = "olr.anom.nc"
plotName = "figures/map.variance"
else
fileName = "olr.std.waves.nc"
plotName = "figures/map." + varName
end if
end if
; plotName = varName
dataList = NewList( "lifo" )
minTotal = -999.
minTotal@_FillValue = minTotal
maxTotal = minTotal
sensor = "20th"
minYear = 1911
maxYear = 1940
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
ancientData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
ancientData@long_name = long_name
minTotal = min( (/ minTotal, min( ancientData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( ancientData({-60:60},:) ) /) )
ListPush( dataList, ancientData )
sensor = "20th"
minYear = 1981
maxYear = 2010
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
modernData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
modernData@long_name = long_name
minTotal = min( (/ minTotal, min( modernData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( modernData({-60:60},:) ) /) )
ListPush( dataList, modernData )
if( False ) then
sensor = "20th"
minYear = 1911
maxYear = 1940
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
ancientData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
ancientData@long_name = long_name
minTotal = min( (/ minTotal, min( ancientData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( ancientData({-60:60},:) ) /) )
ListPush( dataList, ancientData )
sensor = "20th"
minYear = 1981
maxYear = 2010
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
long_name = sensor + " (" + minYear + " - " + maxYear + ")"
print_clock( "Calculating " + long_name )
modernData := calc_meanvar( inPath, varName, minYear, maxYear, meanNotVar )
modernData@long_name = long_name
minTotal = min( (/ minTotal, min( modernData({-60:60},:) ) /) )
maxTotal = max( (/ maxTotal, max( modernData({-60:60},:) ) /) )
ListPush( dataList, modernData )
end if
sensor = "avhrr"
minYear = 1981
......@@ -97,6 +104,7 @@ begin ; main
; Customize base plot
totalRes = True
totalRes@lbOrientation = "Vertical"
totalRes@cnFillPalette = "CBR_wet"
totalRes@cnLevelSelectionMode = "ManualLevels"
if( meanNotVar ) then
......@@ -108,16 +116,20 @@ begin ; main
totalRes@cnMaxLevelValF = mnmxintvl(1)
totalRes@cnLevelSpacingF = mnmxintvl(2)
totalRes@pmLabelBarOrthogonalPosF = 0.30
totalRes@tmXBLabelsOn = False
diffRes = True
diffRes@lbOrientation = "Vertical"
diffRes@cnFillPalette = "CBR_drywet"
diffRes@cnLevelSelectionMode = "ManualLevels"
diffRes@pmLabelBarOrthogonalPosF = 0.30
diffRes@tmXBLabelsOn = True
; Customize panel
panRes = True
panRes@gsnPanelRowSpec = True
panRes@gsnPanelCenter = False
panRes@gsnPanelRowSpec = False
panRes@gsnPanelCenter = True
panRes@gsnPanelYWhiteSpacePercent = 5
print_clock( "Drawing the plot" )
......@@ -168,7 +180,7 @@ begin ; main
end do
end do
gsn_panel( wks, plots, ispan( 1, nData, 1 ), panRes )
gsn_panel( wks, plots, (/ nPlots, 1 /), panRes )
print_clock( "Convert the image, if necessary" )
delete(wks)
......
#!/bin/bash
NCL_SCRIPT=time_series
echo $JOB_NAME `date`
NCL_DIR=`pwd`
mkdir log
for PREFIX in olr mjo er kelvin low
do
for VAR_NAME in total anom
do
NCL_OPTION="varName=\"$VAR_NAME\" prefix=\"$PREFIX\""
JOB_NAME="$PREFIX"."$VAR_NAME"
qsub \
-N $JOB_NAME \
-o $NCL_DIR/log/$JOB_NAME.log -j oe \
-q san \
-l walltime=12:00:00,nodes=1:ppn=1 \
-v NCL_DIR=$NCL_DIR,NCL_SCRIPT=$NCL_SCRIPT,NCL_OPTION="$NCL_OPTION" \
$NCL_DIR/run_ncl.pbs
done
done
#!/bin/bash
NCL_SCRIPT=meanvar
echo $JOB_NAME `date`
NCL_DIR=`pwd`
mkdir log
for VAR_NAME in olr mjo er kelvin low
do
NCL_OPTION="varName=\"$VAR_NAME\""
JOB_NAME="map.$VAR_NAME"
qsub \
-N $JOB_NAME \
-o $NCL_DIR/log/$JOB_NAME.log -j oe \
-q san \
-l walltime=12:00:00,nodes=1:ppn=1 \
-v NCL_DIR=$NCL_DIR,NCL_SCRIPT=$NCL_SCRIPT,NCL_OPTION="$NCL_OPTION" \
$NCL_DIR/run_ncl.pbs
done
......@@ -9,11 +9,12 @@ mkdir log
for VAR_NAME in olr mjo er kelvin low
do
for SENSOR in 20th avhrr hirs
for SENSOR in avhrr hirs 20th
do
NCL_OPTION="varName=\"$VAR_NAME\" sensor=\"$SENSOR\""
JOB_NAME="$SENSOR"."$VAR_NAME"
qsub \
-N $JOB_NAME \
-o $NCL_DIR/log/$JOB_NAME.log -j oe \
......
......@@ -31,7 +31,8 @@ function draw_series( \
)
begin ; draw_series
; These are some parameters that could be useful to have up top
nSmooth = 365*5
shortSmooth = 365
longSmooth = 365*5
inPath = "/home/carl/data/olr/compare/time_series/" + i_prefix + "." \
+ i_sensor + ".20S-20N.nc"
res = i_res
......@@ -70,7 +71,8 @@ begin ; draw_series
end do
print_clock("Smoothing")
smoothData = runave_n_Wrap( yData, nSmooth, 0, 1 )
yData = runave_n_Wrap( yData, shortSmooth, 0, 1 )
smoothData = runave_n_Wrap( yData, longSmooth, 0, 1 )
print_clock("Plotting")
resTick = True
......@@ -79,11 +81,13 @@ begin ; draw_series
month_axis_labels( i_minTime, 60, 12, 1, res, resTick )
res@xyLineThicknessF = 1
retVal = cjs_draw_xy( io_wks, xData, yData, res )
rawPlot = cjs_draw_xy( io_wks, xData, yData, res )
res@xyLineThicknessF = 3
smoothPlot = cjs_draw_xy( io_wks, xData, smoothData, res )
overlay( retVal, smoothPlot )
overlay( rawPlot, smoothPlot )
retVal = rawPlot
return(retVal)
end; draw_series
......
......@@ -23,9 +23,13 @@ begin ; main
print_clock( "Here we go!" )
; These are some parameters that could be useful to have up top
prefix = "low"
sensor = (/ "hirs", "avhrr", "20th" /)
varName = "total"
if(.not.isvar("prefix") )then
prefix = "olr"
end if
if(.not.isvar("varName") )then
varName = "anom"
end if
sensor = (/ "hirs", "avhrr" /)
plotType = "png"
plotName = "figures/" + prefix + "." + varName
......
......@@ -42,6 +42,7 @@ begin ; main
maxLat = 20
inPath = "/home/carl/data/olr/compare/" + sensor + "/" + fileName
rawPath = "/home/carl/data/olr/compare/" + sensor + "/olr.total.nc"
outPath = "/home/carl/data/olr/compare/time_series/" + varName + "." \
+ sensor + "." + lat_string(minLat) + "-" + lat_string(maxLat) + ".nc"
......@@ -51,20 +52,29 @@ begin ; main
print_clock( "Reading" )
rawFile = addfile( rawPath, "r" )
rawData = lon_subset( rawFile->olr(:,{minLat:maxLat},:), \
minLon, maxLon )
inFile = addfile( inPath, "r" )
inData = lon_subset( inFile->$varName$(:,{minLat:maxLat},:), minLon, maxLon )
time = inFile->time
print_clock( "Calculating series" )
inData = inData^2
inData@units = "(" + inData@units + ")^2"
totalSeries = dim_avg_n_Wrap( inData, (/ 1, 2 /) )
totalSereis = sqrt( totalSeries )
; totalSereis = sqrt( totalSeries )
missSeries = dim_num_n( ismissing(rawData), (/ 1, 2 /) )
totalSeries = where( missSeries.ne.0, totalSeries@_FillValue, totalSeries )
print_clock( "Standardizing series" )
mean = avg( totalSeries({climStart:climEnd}) )
std = stddev( totalSeries({climStart:climEnd}) )
anomSeries = ( totalSeries - mean ) / std
copy_VarMeta( totalSeries, anomSeries )
anomSeries = where( ismissing(anomSeries), 0, anomSeries )
anomSeries@units = "Std. Dev."
format = "%8.2f"
print( cd_string(time,"") + " " + sprintf( format, totalSeries ) \
......
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