;First define the file name to be opened. you will need to change this to the file location in your directory
args = command_line_args()
print, args[0]
file1='command_line_args(a)'
;;;;;;idl opens the file
id=NCDF_OPEN(file1)
NCDF_VARGET,id,'trackid',trackno_xbylr
NCDF_VARGET,id,'variables',vars_xbylr
NCDF_VARGET,id,'data', data_xbylr
NCDF_VARGET,id,'points', points_xbylr
NCDF_CLOSE,id
num_tracks=n_elements(trackno_xbylr)
print, num_tracks
num_points=n_elements(points_xbylr)
print, num_points
;Define variables and create arrays
vort=fltarr(num_tracks,100) ; vorticity at each point on each track
vortcolour=strarr(num_tracks,100) ; colour for vorticity value at each point on each track
;;;;If you would like to check it's working so far you can remove the comments here and see the vorticity values it prints
;sanity check? Should show vort vals for first 20 points on first 20 tracks
;print, data_xbylr(0:20,0:20,2)
;;;;For missing values
data_xbylr(where(data_xbylr eq 0.0))=!VALUES.F_NAN
; open postscript file - Andy Heaps syntax used from here---from the IDL guide
;
; Adjust file name to save image in your directory
PSOPEN, file='command_line_args(b)'
;;;;Map.pro gives you options on projections
;default is MAP,/NH, can reduce the size if give lat/long, hires = high resolution coastlines
MAP,/cylindrical, OCEAN=10,LAND=11,/draw,latmin=30,latmax=80,lonmin=-90,lonmax=10 ;;can adjust the lat/long to reduce plot size
;MAP,/sector,OCEAN=10,LAND=11,/draw ;;;I modified sector in map.pro which i would not recommend doing. See if you can adjust the lat/long like in other projections
;Colorscale default is 1 but here we are specifying the number of colours on scale 4 we would like to use to show intensity in vorticity
CS,SCALE=4,NCOLS=254
;;;;;;;;;;Loop thru tracks to assign vorticity into intensity bins
for i=0, num_tracks-1 do begin
; loop through track number
; Max jvalue is 2 less than total points - as loops onto next point when plotting
for j=0, 98 do begin
; loop through points on track
; create array for the vorticity values, which will define the track colouring
vort(i,j) = data_xbylr(i,j,2)
;;;;;;;;;;;;; AMEND THE BIN VALUES DEPENDING ON THE INTENSITY DISTRIBUTIONS
if (vort(i,j) gt 0) and (vort(i,j) le 2) then begin
vortcolour(i,j) = 2
endif
if (vort(i,j) gt 2) and (vort(i,j) le 4) then begin
vortcolour(i,j) = 50
endif
endif
if (vort(i,j) gt 4) and (vort(i,j) le 6) then begin
vortcolour(i,j) = 100
endif
if (vort(i,j) gt 6) and (vort(i,j) le 8) then begin
vortcolour(i,j) = 150
endif
if (vort(i,j) gt 8) and (vort(i,j) le 10) then begin
vortcolour(i,j) = 200
endif
if (vort(i,j) gt 10) then begin
vortcolour(i,j) = 250
endif
;VOR='!9'+SCROP(BYTE(120))+'!X'
intensity=['VOR<2','2<VOR<4','4<VOR<6','6<VOR<8','8<VOR<10','10<VOR']
cols=[2,50,100,150,200,250]
;;;;;Plots tracks on map
GPLOT, X=[data_xbylr(i,j,0), data_xbylr(i,j+1,0)], Y=[data_xbylr(i,j,1), data_xbylr(i,j+1,1)], THICK=80, COL=vortcolour(i,j)
endfor
endfor
; plot x on first point of track to indicate direction
; we don't want x to indicate tracks outside our region of interest
for i=0, num_tracks-1 do begin
inds_lowerbound=where(data_xbylr(*,0,1) lt 30.) ; this will remove x marks below 30N
inds_upperbound=where(data_xbylr(*,0,1) gt 70.) ; this will remove x marks above 70N
inds_westbound=where(data_xbylr(*,0,0) lt 270.) ; this will show x marks for 0 to 90 W, but anything east of 0 up to 90W has no x
;inds_eastbound=where(data_xbylr(*,0,0) gt 10. ) ;here if needed to adjust x marks in area of interest
new_data_xbylr=data_xbylr
new_data_xbylr(inds_westbound)=!values.f_nan
new_data_xbylr(inds_upperbound)=!values.f_nan
new_data_xbylr(inds_upperbound)=!values.f_nan
new_data_xbylr(inds_lowerbound)=!values.f_nan
GPLOT, X=new_data_xbylr(i,0,0), Y=new_data_xbylr(i,0,1), SYM=6, SIZE=50
endfor
;Can change value of LEGPOS to change legend position.
GLEGEND, LEGPOS=9, LEGXOFFSET=1500,COL=cols, LABELS=intensity, FILLCOL=0, TYPE=0 ; TYPE=1 has legend with thicker colour blocks and TYPE=2 has thicker blocks without border around each block, LEGXOFFSET=-5000 to move legend left
;;;Plots Title on image. make sure you change this to match which season you are plotting
GPLOT, X=14752,Y=19000, CHARSIZE=200,/DEVICE, TEXT='command_line_args(c) DJF season vorticity maximums'
;;can adjust labels on axes. i believe default is 30 and 30
AXES, XSTEP=15, YSTEP=10
PSCLOSE
end