aofd165.bio.mie-u.ac.jp
/work1/am/2016.PolarLow/AVHRR/Arctic.Polar.Pathfinder
Fri Jun 24 14:10:29 JST 2016
======================
avhrr.run.sh
======================
#!/bin/bash
exe=./runncl.sh
ncl=avhrr.ncl
indir='/work4/data/AVHRR/Polar.Pathfinder'
inlist="\
Polar-APP_v01r01_Nhem_0400_d20110118_c20150111.nc \
Polar-APP_v01r01_Nhem_0400_d20110119_c20150111.nc \
Polar-APP_v01r01_Nhem_0400_d20110120_c20150111.nc \
Polar-APP_v01r01_Nhem_0400_d20110121_c20150111.nc \
Polar-APP_v01r01_Nhem_0400_d20110122_c20150111.nc \
Polar-APP_v01r01_Nhem_1400_d20110118_c20150111.nc \
Polar-APP_v01r01_Nhem_1400_d20110119_c20150111.nc \
Polar-APP_v01r01_Nhem_1400_d20110120_c20150111.nc \
Polar-APP_v01r01_Nhem_1400_d20110121_c20150111.nc \
Polar-APP_v01r01_Nhem_1400_d20110122_c20150111.nc \
"
#inlist="\
#Polar-APP_v01r01_Nhem_0400_d20110118_c20150111.nc \
#"
for infle in $inlist; do
figfile=$(basename $infle .nc).png
$exe $ncl "$indir" "$infle" "$figfile"
done
exit 0
----------------------
End of avhrr.run.sh
----------------------
======================
runncl.sh
======================
#!/bin/bash
#
# Universal wrapper script for ncl.
# Pass arguments from the command line to environment variables
#
# version 0.1, Thierry Corti, C2SM ETH Zurich
#
E_BADARGS=65
if [ ! -n "$1" ]
then
echo "Usage: `basename $0` script.ncl argument1 argument2 etc."
exit $E_BADARGS
fi
# save number of arguments to environment variable NCL_N_ARG
export NCL_N_ARGS=$#
# save command line arguments to environment variable NCL_ARG_#
for ((index=1; index<=$#; index++))
do
eval export NCL_ARG_$index=\$$index
done
# run ncl
ncl -nQ $1
----------------------
End of runncl.sh
----------------------
======================
avhrr.ncl
======================
;*************************************************
; avhrr.ncl
;
; Concepts illustrated:
; - Contouring AVHRR data
; - Drawing raster contours
; - Manually creating lat/lon coordinate arrays
; - Converting "byte" data to "float"
;
;======================================================================
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"
;********************************************
;----------------------------------------------------------------------
;
; Subroutine
;
; This procedure attaches lat/lon labels to a masked lambert plot
;
; You will likely need to change lat_values and/or lon_values to
; contain the locations where you want lat/lon labels.
;----------------------------------------------------------------------
procedure add_lc_labels(wks,map,minlat,maxlat,minlon,maxlon)
local lat_values, nlat, lat1_ndc, lat2_ndc, lon1_ndc, lon2_ndc,slope,txres, \
lon_values, PI, RAD_TO_DEG, dum_lft, dum_rgt, dum_bot
begin
PI = 3.14159
RAD_TO_DEG = 180./PI
;---Pick some "nice" values for the latitude labels.
lat_values = ispan(toint(minlat),toint(maxlat),10) * 1.
nlat = dimsizes(lat_values)
;
; We need to get the slope of the left and right min/max longitude lines.
; Use NDC coordinates to do this.
;
lat1_ndc = new(1,float)
lon1_ndc = new(1,float)
lat2_ndc = new(1,float)
lon2_ndc = new(1,float)
datatondc(map,minlon,lat_values(0),lon1_ndc,lat1_ndc)
datatondc(map,minlon,lat_values(nlat-1),lon2_ndc,lat2_ndc)
slope_lft = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc)
datatondc(map,maxlon,lat_values(0),lon1_ndc,lat1_ndc)
datatondc(map,maxlon,lat_values(nlat-1),lon2_ndc,lat2_ndc)
slope_rgt = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc)
;---Set some text resources
txres = True
txres@txFontHeightF = 0.015
txres@txPosXF = 0.1
;
; Loop through lat values, and attach labels to the left and
; right edges of the masked LC plot. The labels will be
; rotated to fit the line better.
;
dum_lft = new(nlat,graphic) ; Dummy array to hold attached strings.
dum_rgt = new(nlat,graphic) ; Dummy array to hold attached strings.
do n=0,nlat-1
; Add extra white space to labels.
lat_label_rgt = " " + lat_values(n) + "~S~o~N~"
;---Check if North, South, or Zero
if(lat_values(n).lt.0) then
lat_label_lft = lat_values(n) + "~S~o~N~S "
lat_label_rgt = lat_label_rgt + "S"
end if
if(lat_values(n).gt.0) then
lat_label_lft = lat_values(n) + "~S~o~N~N "
lat_label_rgt = lat_label_rgt + "N"
end if
if(lat_values(n).eq.0) then
lat_label_lft = lat_values(n) + "~S~o~N~ "
end if
;---Left label
txres@txAngleF = RAD_TO_DEG * atan(slope_lft) - 90
dum_lft(n) = gsn_add_text(wks,map,lat_label_lft,minlon,lat_values(n),txres)
;---Right label
txres@txAngleF = RAD_TO_DEG * atan(slope_rgt) + 90
dum_rgt(n) = gsn_add_text(wks,map,lat_label_rgt,maxlon,lat_values(n),txres)
end do
;----------------------------------------------------------------------
; Now do longitude labels. These are harder because we're not
; adding them to a straight line.
;
; Loop through lon values, and attach labels to the bottom edge of the
; masked LC plot.
;
delete(txres@txPosXF)
txres@txPosYF = -5.0
;---Pick some "nice" values for the longitude labels.
lon_values = ispan(toint(minlon+10),toint(maxlon-10),10) * 1.
nlon = dimsizes(lon_values)
dum_bot = new(nlon,graphic) ; Dummy array to hold attached strings.
do n=0,nlon-1
;
; For each longitude label, we need to figure out how much to rotate
; it, so get the approximate slope at that point.
;
datatondc(map,lon_values(n)-0.25,minlat,lon1_ndc,lat1_ndc)
datatondc(map,lon_values(n)+0.25,minlat,lon2_ndc,lat2_ndc)
slope_bot = (lat1_ndc-lat2_ndc)/(lon1_ndc-lon2_ndc)
txres@txAngleF = atan(slope_bot) * RAD_TO_DEG
;
; Create longitude label. Add extra carriage returns to
; move label away from plot.
;
;---Check if East, West, or Zero
lon_label_bot = " ~C~ ~C~" + abs(lon_values(n)) + "~S~o~N~"
if(lon_values(n).lt.0) then
lon_label_bot = lon_label_bot + "W"
end if
if(lon_values(n).gt.0) then
lon_label_bot = lon_label_bot + "E"
end if
;---Attach to map.
dum_bot(n) = gsn_add_text(wks,map,lon_label_bot,lon_values(n),minlat,txres)
end do
end
;----------------------------------------------------------------------
; Main code.
;----------------------------------------------------------------------
begin
wallClock1 = systemfunc("date") ; retrieve wall clock time
scriptname = getenv("NCL_ARG_1")
; infle="Polar-APP_v01r01_Nhem_0400_d20110118_c20150111.nc"
; indir="/work4/data/AVHRR/Polar.Pathfinder"
; figfile="Polar-APP_v01r01_Nhem_0400_d20110118_c20150111.png"
indir = getenv("NCL_ARG_2")
infle = getenv("NCL_ARG_3")
figfile = getenv("NCL_ARG_4")
; print(indir)
print("Processing " + infle + "...")
input=indir + "/" + infle
f = addfile(input, "r")
; print(f)
bt1= f->fcdr_brightness_temperature_ch4
lon = f->longitude
lat = f->latitude
dsizes = new( (/3/),integer)
dsizes = dimsizes(bt1)
; print(dsizes)
im=dsizes(1)
jm=dsizes(2)
bt=new( (/im,jm/), float) ;typeof(bt1))
bt(:,:)=short2flt( bt1(0,:,:) )- 273.15 ;*bt1@scale_factor
opts = True
opts@gsnMaximize = False
opts@gsnDraw = False ; Don't draw yet
opts@gsnFrame = False ; Don't advance frame yet
opts@gsnMaskLambertConformal = True
opts@gsnAddCyclic = False
opts@sfXArray = lon
opts@sfYArray = lat
; Figure size
opts@vpXF = 0.05
opts@vpYF = 1
opts@vpHeightF = 0.8
opts@vpWidthF = 0.8
opts@cnFillOn = True
opts@cnLinesOn = False ; turn off contour lines
opts@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
opts@cnMinLevelValF = -60 ; set min contour level
opts@cnMaxLevelValF = 0 ; set max contour level
opts@cnLevelSpacingF = 5 ; set contour spacing
opts@pmLabelBarOrthogonalPosF = 0.0
;Map
opts@mpFillOn = False
opts@mpProjection = "LambertConformal"; choose projection
opts@mpLimitMode = "LatLon"
opts@pmTickMarkDisplayMode = "Always"
opts@mpDataBaseVersion = "MediumRes"
opts@mpGeophysicalLineColor = "Red" ;changes the outline line color.
opts@mpGeophysicalLineThicknessF = 3 ;changes the thickness of continental outlines.
minlon = -30
maxlon = 35
minlat = 60
maxlat = 80
opts@mpMinLonF = minlon
opts@mpMaxLonF = maxlon
opts@mpMinLatF = minlat
opts@mpMaxLatF = maxlat
type="png"
; figfile=infle + "." + type
; print(bt)
wks = gsn_open_wks(type, figfile)
gsn_define_colormap(wks, "MPL_Greys")
opts@pmLabelBarOrthogonalPosF = .05 ; move whole thing down
setvalues NhlGetWorkspaceObjectId()
"wsMaximumSize" : 400000000
end setvalues
plot = gsn_csm_contour_map(wks,bt,opts)
; plot = gsn_csm_map(wks,opts)
;---Attach latitude labels
add_lc_labels(wks,plot,minlat,maxlat,minlon,maxlon)
; Header and footer
;
; drawNDCGrid(wks)
txres = True ; text mods desired
txres@txFontHeightF = 0.02 ; font smaller. default big
gsn_text_ndc(wks,infle,0.45,0.98,txres)
txres@txFontHeightF = 0.01 ; font smaller. default big
text=systemfunc("pwd")
gsn_text_ndc(wks,text,0.2,0.010,txres)
text=scriptname
gsn_text_ndc(wks,text,0.5,0.03,txres)
;---Drawing the plot will also draw all the attached labels.
draw(plot)
frame(wks)
wallClock2 = systemfunc("date") ; retrieve wall clock time
print("Started at " + wallClock1)
print("Finished at " + wallClock2)
print("Input: " + infle)
print("Output: " + figfile)
end
----------------------
End of avhrr.ncl
----------------------