source: trunk/src/basic_variable_evaluation_tpr_timeseries.pro

;+
;
; .. _basic_variable_evaluation_tpr_timeseries.pro:
;
; ============================================
; basic_variable_evaluation_tpr_timeseries.pro
; ============================================
;
; DESCRIPTION
; ===========
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
; :ref:`project_init.pro`
; :ref:`cm_project.pro`
;
; use :func:`read_basic_var`, :func:`x_site_location`, :func:`y_site_location`
;
; TODO
; ====
;
; make it work : missing extract_tpr_location.pro
;
; uncomment var parameter
;
; correct bad check on unavailable tpr data
;
; EXAMPLES
; ========
;
; .. code-block:: idl
;
;    date1 = 20000101L
;    date2 = 20091231L
;    basic_variable_evaluation_tpr_timeseries, date1, date2
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20111130T133857Z cratos (Linux)
;
;   * try to make it work on my account on cratos
;   * remove x_site_location and  x_site_location because already exist
;
;-
pro basic_variable_evaluation_tpr_timeseries, $
;                   var,    $ ;;  basic variable (sst, t2m, q2m, ws) to calculate the statistics
                    date1,  $ ;;  start date (in Julian date. eg. 20000101)
                    date2     ;;  end date (in Julian date eg. 20091231)
@common
@cm_project
;------------------------------------------------------------
reinitplt, /z,/invert
key_portrait = 0
openps, FILENAME = project_od_env+'basic_variable_evaluation_tpr_timeseries.ps'
;------------------------------------------------------------
;; part to change

min_obs=10.  ;; this will allow to calculate statistics at locations with more than 180. valid observation

;------------------------------------------------------------
;; Before running this program, you have to compile the following subroutines
;; .r read_basic_var

;; TPR locations.  This needs to be updated with time since more locations are added to the array.

sitelist=['5n165e','8s67e','12s55e', '8s55e', '8s80.5e', '1.5s80.5e', '0n80.5e', '1.5n80.5e', '1.5s90e', $
           '0n90e', '1.5n90e', '4n90e','8n90e','12n90e', '15n90e', '5s95e', $
           '8s165e', '8s180w',  '8s155w', '8s125w', '8s110w', '8s95w',  '5s156e', '5s165e', '5s180w', '5s170w', $
          '5s155w', '5s140w', '5s125w', '5s110w', '5s95w', '2s156e', '2s165e', '2s180w', '2s170w', '2s155w', '2s140w', $
          '2s125w', '2s110w', '2s95w', '0n147e', '0n156e', '0n165e', '0n180w', '0n170w', '0n155w', '0n140w', '0n125w', $
          '0n110w', '0n95w', '2n147e', '2n156e', '2n165e', '2n180w', '2n170w', '2n155w', '2n140w', '2n125w', '2n110w', $
          '2n95w', '5n147e', '5n156e', '5n170w', '5n155w', '5n140w', '5n125w', '5n110w', '5n95w', $
          '8n156e', '8n165e', '8n180w', '8n170w', '9n140w', '8n125w', '8n110w', '8n95w', $
          '0n0e', '0n10w', '0n23w', '0n35w', '10s10w', '12n23w', '12n38w', '14s32w', '15n38w', '19s34w', '20n38w', $
          '21n23w', '4n23w', '4n38w', '6s10w', '8n38w', '8s30w']

;------------------------------------------------------------------------------------------------------------------------
;;   This program will create the following text files with statistics of respective variables
;------------------------------------------------------------------------------------------------------------------------
close,/all

fi = project_id_env + 'basic_var_stat.txt'
openw,1,fi

printf,1, 'x     y      cor    bias     std ratio     rmsd'
;------------------------------------------------------------------------------------------------------------------------
;; First, this program reads the full TropFlux data and later extract it at specific TPR locations

file = project_id_env + "sst_tropflux_1d_1989_2010.nc"
initncdf, file
var=read_ncdf("sst", date1, date2, file=file,/nostr)
help, var

;------------------------------------------------------------------------------------------------------------------------
nn=n_elements(sitelist)
date1=date1
date2=date2
nsmooth=1. ;; this will return daily TPR values
mooring=0 & product=0

for n=0, nn-1 do begin

;; reading data from mooring

    site=sitelist(n) & csite=site
    print, csite
    x=x_site_location(site)
    y=y_site_location(site)
    if (y ge 0. and y le 30.) then y=y+360.
    dx=0.5 & dy=0.5 & box=[y-dy, y+dy, x-dx, x+dx]

    read_basic_var, csite,date1,date2,nsmooth, $
            at, rh, sst, wu, wv, ws
help, sst
print, sst
;; select the appropriate variables for evaluation (trp = sw or lw or sh or lh)

   tpr=sst

   ind=where(finite(tpr)) & no_valid=n_elements(ind)

   if (no_valid ge min_obs) then begin


       extract_tpr_location,var,box, $
       var_tpr
       var_tpr=reform(var_tpr)

       stats_5d, tpr,var_tpr, $  ;; tpr=TPR observation and var_tpr=gridded product extracted at TPR location
       cor, bias, std, rmsd

       printf, 1, x, y, cor, bias, std, rmsd, format='(f6.2, 3x, f6.2, 3x, f4.2, 3x, f7.2, 3x, f4.2, 3x, f5.2)'
       cstat=string(cor, bias, std, rmsd, format='(f4.2,3x,f7.2,3x,f4.2,3x,f5.2)')
       print, cstat

;;     PLOTTING THE TIME-SERIES
       array=[tpr, var_tpr] & mi=min(array,/nan) & ma=max(array,/nan) & int=(ma-mi)/3.
       pltt, ts_smooth(tpr,5,/nan), "t",/rempl, small=[1,3,1], lct=65, $
       title='Five day stats are shown below.  TPR (black) and Product (red) at'+csite+' ', charsize=1., $
       subtitle=cstat
       ind=where(finite(tpr,/nan)) & var_tpr(ind)=!Values.f_nan
       pltt, ts_smooth(var_tpr,5,/nan), "t",/ov1d, color=250
       erase
       mooring=[mooring,tpr] & product=[product,var_tpr]
    endif
endfor

close,/all
closeps

end