IDL+code

=;Here we (together) build the code for IDL to do the homework.=

;Data Import

 * I recommend [[file:cdf2idl.pro]] for quick grabs of all the data in a file!
 * Caution: data are packed as short ints and may need to be unpacked.
 * hlp.pro is useful for seeing the dimensions of arrays and the max and min values. [[file:hlp.pro]]
 * I use make_key.pro to make a fast simple color bar. [[file:make_key.pro]]
 * I used map_world.pro to overlay continents on my plots. [[file:map_world.pro]]


 * Land.pro and psout.pro are used for plotting. [[file:land.pro]] [[file:psout.pro]]


 * ;;;;Code to complete MPO581 HW 1 **
 * ;;;;ER - 2/3/2011 **
 * ;;;;NOTE: The codes land.pro, psout.pro, hlp.pro, make_key.pro , and map_world.pro were used in this code. **
 * ;;;; They can be found on the IDL code page. **

pro go

land ;;; This uses land.pro device, file = 'HW1_basicstats.ps'

;;; If you don't want to use land.pro the code is: ; set_plot,'ps' ; device,/color,bits_per_pixel=8 ; loadct, 39 ; Rainbow + white ; device,/landscape

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; 1. Download and Read the data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; I just downloaded the data to my computer

;;;; Convert netcdf to idl - I do this at command line level ;;;; make sure the path to .nc is correct ; cdf2idl, '/Users/mapesgroup/MPO581/Data/air.mon.ltm.nc' ;air T ; cdf2idl, '/Users/mapesgroup/MPO581/Data/olr.mon.ltm.nc' ;olr ; cdf2idl, '/Users/mapesgroup/MPO581/Data/prate.mon.ltm.nc' ;rain gauge ; cdf2idl, '/Users/mapesgroup/MPO581/Data/precip.mon.ltm.nc' ;model precip ; cdf2idl, '/Users/mapesgroup/MPO581/Data/uwnd.mon.ltm.nc' ;u wind ; cdf2idl, '/Users/mapesgroup/MPO581/Data/vwnd.mon.ltm.nc' ;v wind

;;;Example: ; **IDL>** cdf2idl, '/Users/mapesgroup/MPO581/Data/prate.mon.ltm.nc' ;rain gauge ; You may now type @prate.mon.ltm.idl at the IDL prompt to input data from /Users/mapesgroup/MPO581/Data/prate.mon.ltm.nc.

;;;; Let's acess the data ;;;;Check the dimensions of each array, use ncdump -h in terminal or ;;;;hlp from command line ; prate and precip have diff dims from the others, so I want to load ; them first and rename their lat & lon dimensions, so they ; don't get overwritten later. @prate.mon.ltm.idl ;; [192, 94, 12] lat_prate = lat lon_prate = lon

@precip.mon.ltm.idl ;; [144, 72, 12] lat_precip = lat lon_precip = lon

;; air, olr, and uwnd, vwnd have some dimensions, will leave their lat & ;; lon alone @air.mon.ltm.idl ;; [144, 73, 12] @olr.mon.ltm.idl @uwnd.mon.ltm.idl @vwnd.mon.ltm.idl

;;;;;; See what I have: ;;; **IDL>** help ;;; % At $MAIN$ ;;; AIR FLOAT = Array[144, 73, 12] ;;; LAT FLOAT = Array[73] ;;; LAT_PRATE FLOAT = Array[94] ;;; LAT_PRECIP FLOAT = Array[72] ;;; LON FLOAT = Array[144] ;;; LON_PRATE FLOAT = Array[192] ;;; LON_PRECIP FLOAT = Array[144] ;;; NCID LONG = 12 ;;; OLR INT = Array[144, 73, 12] ;;; PRATE FLOAT = Array[192, 94, 12] ;;; PRECIP FLOAT = Array[144, 72, 12] ;;; TIME DOUBLE = Array[12] ;;; UWND FLOAT = Array[144, 73, 12] ;;; VWND FLOAT = Array[144, 73, 12]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; 2. Averaging p and v over longitude ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

s = size(precip) ;;; gives number of dims, size of each dim, dim type, and total elements

;;; Average precip and vwnd over longitude precip_xbar = total(precip, 1)/n_elements(lon) ;;; or could put s[1] as denominator vwnd_xbar = total(vwnd, 1)/n_elements(lon) ;;; note n_elements(lon) is a long int, but vwnd is a float, so no worries about improper dividing

;;;;;;;;;;;;;;;;;;;;;; ;;; Plot it up

;;; create time axis array time_axis = indgen(12) + 1

;;; Plot precip_xbar, without making it pretty contour, precip_xbar, lat_precip, time_axis, /fill, nlev = 33

;;; Oh-but Brian wants time on the x-axis and lat on the y, better ;;; rotate the array, and let's add title and axis title and label the contours. spacing = (max(precip_xbar) - min(precip_xbar))/10. levels = min(precip_xbar) + findgen(11)*spacing ;;; create levels for color shading. these levels are exact, but not nice

levels = findgen(10) ;;; much nicer levels

contour, rotate(precip_xbar, 4), time_axis, lat_precip, /fill, levels = levels, $ title = 'Precipitation Rate', xtitle = 'Month' , ytitle = 'latitude' contour, rotate(precip_xbar, 4), time_axis, lat_precip, /over, levels = levels, $ c_labels = intarr(10)+1 ;;; label the contours

;;; Plot vwnd_xbar levels = -4 + findgen(17)*.5 ;;this looks nice, but you can play with amounts of levels and spacing

;;; Perhaps I only want every other contour labeled labels = intarr(17) ;;; first define an array of zeros for i=0, 16 do $ if i mod 2 eq 0 then labels[i] = 1 ;;; now fill everyother element with 1

;;; Since we have both neg and pos winds let's make the neg ;;; dashed contours and thicken the zero contour contour_style = (levels lt 0)*2 thick_array = intarr(17) + 3 thick_array[where(levels eq 0)] = 5

contour, rotate(vwnd_xbar, 4), time_axis, lat, /fill, levels = levels, $ title = 'Meridional Wind', xtitle = 'Month' , ytitle = 'latitude' contour, rotate(vwnd_xbar,4), time_axis, lat, /over, levels = levels, $ c_labels = labels, c_linestyle = contour_style, $ c_thick = thick_array ;;; label the contours

;;; more axis labels values = indgen(19)*10-90 contour, rotate(vwnd_xbar, 4), time_axis, lat, /fill, levels = levels, $ title = 'Meridional Wind', xtitle = 'Month' , ytitle = 'latitude' , $ yticks = 20, ytickv = values, ytickn = str(values)

contour, rotate(vwnd_xbar,4), time_axis, lat, /over, levels = levels, $ c_labels = labels, c_linestyle = contour_style, $ c_thick = thick_array ;;; label the contours

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Make comparison plots for 90°W ;;; find array index where lon = 90°W lon90 = where(lon eq 90 + 180) vwnd90 = vwnd[lon90, *, *] vwnd90 = reform(vwnd90) ;;; eliminate single element dimensions

;;; find array index where lon_precip = 90°W ;;; lon_precip = where(lon eq 90) ;returns -1 cause index doesn't exist! ;;; So find what index is closest to 90°W dummy = min(abs(lon_precip - (90+180)), index) precip90 = reform(precip[index, *, *]) ;;; now I've done the reform all in one step

;;;;;;;;;;;;;;;;;;;;;; ;;; Plot it up

;;; Precip levels = findgen(17) ;;; again makes some color levels

contour, rotate(precip90, 4), time_axis, lat_precip, /fill, levels = levels, $ title = 'Precipitation Rate at 90!Uo!NW', xtitle = 'Month' , ytitle = 'latitude' contour, rotate(precip90, 4), time_axis, lat_precip, /over, levels = levels, $ c_labels = intarr(17)+1 ;;; label the contours

;;; V Wind levels = -5 + findgen(14) ;;; this looks nice, but you can play with amounts of levels and spacing labels = intarr(14) + 1 contour_style = (levels lt 0)*2 thick_array = intarr(14) + 3

contour, rotate(vwnd90, 4), time_axis, lat, /fill, levels = levels, $ title = 'Meridional Wind at 90!Uo!NW', xtitle = 'Month' , ytitle = 'latitude' contour, rotate(vwnd90,4), time_axis, lat, /over, levels = levels, $ c_labels = labels, c_linestyle = contour_style, $ c_thick = thick_array ;;; label the contours

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; 3. Air (temp) time series ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; First do lon average air_xbar = total(air, 1)/n_elements(lon)

;;; Now lat average -- but have to weight the latitudes coslat = cos(lat *!pi / 180.) weight = coslat#(fltarr(12)+1) ;;; make coslat same dimensions as air_xbar by doing matrix multiplication air_xybar = total(air_xbar*weight, 1) / total(coslat)

;;; I want to make these plots smaller, so I plot 2X2 figs per page !p.multi = [0, 2, 2]

plot, air_xybar, title = 'Temp Time Series', xtitle = 'Month' , $ ytitle = 'Temp (!Uo!NC)', thick = 5

;;; Make the y axis go from 0-16 plot, air_xybar, title = 'Temp Time Series', xtitle = 'Month' , $ ytitle = 'Temp (!Uo!NC)', yrange = [0,16], thick = 5

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; 4. Temporal stdev of Precipitation ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; First find the temporal average precip_bart = total(precip, 3)/n_elements(time)

;;; Find variance Var(X) = E[(X - mu)^2] = E(X^2) - mu^2 precip_var = total(precip^2, 3)/12 - precip_bart^2

;;; Find standard deviation sd = sqrt(precip_var)

;;; Get the answer mystat = sd/precip_bart

;;; Plot it up ;;; I no longer want 2X2 plots per page, but 1 plot per page !p.multi = 0

levels20 = findgen(21)/10. levels10 = findgen(11)/5. contour, mystat, lon_precip, lat_precip, /fill, levels = levels20, $ title = '"Monsoon Climate"', xtitle = 'longitude' , ytitle = 'latitude'

map_world, /cont , /over , c_thick = 7

;;; Add a color bar MAKE_KEY, COLORS= findgen(n_elements(levels20)) * !D.N_colors /n_elements(levels20), $ labels = round_any([min(levels20),(min(levels20)+max(levels20))/4.,(min(levels20)+max(levels20))/2., $ (min(levels20)+max(levels20))*3/4.], sig= 2), $ units = '%', /orientation, charsize = .75

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; 5. Space-time standard deviation of 'air' (T) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; air_bart = total(air, 3)/n_elements(time) air_bartx = total(air_bart, 1)/n_elements(lon) air_bartxy = total(air_bartx*coslat)/total(coslat)

air2_bart = total(air^2, 3)/n_elements(time) air2_bartx = total(air2_bart, 1)/n_elements(lon) air2_bartxy = total(air2_bartx*coslat)/total(coslat)

air_var = air2_bartxy - air_bartxy^2 air_sd = sqrt(air_var)

psout ;;; If you don't use psout.pro, the code is: ; device,/close ; x

stop end