;========================================================================= ; I've tried sum 12 hour precipitation to obtain 24 hour precipitation using attached script. ; However it seems doesn't work correctly! ; I had to sum tp(0,:,:)+tp(1,:,:), also tp(2,:,;)+tp(3,;,:), and so on. ; I uploaded "2_2_ERAI_Prec_2015_Jan_00_12_12.nc" in ftp.cgd.ucar.edu. ; I'll be thankful if I hear from you. ;========================================================================= ; About summation of 12 hours precipitation to obtain daily precipitation, ; as suggested me I used "cdo" and it seems work correctly: ; "cdo -b F32 -f nc4c -z zip timselsum,2 2_2_ERAI_Prec_2015_Jan_00_12_12.nc 2_2_ERAI_Prec_2015_Jan_24_cdo.nc4" ; However I like to know how writing nc files with "Geo2D". ;========================================================================= ;Email: Email reply from: Dennis Shea ;Email: I am going to add a few more items to a previous response: ;Email: ;Email: The Golden Rules of data processing are: ;Email: ;Email: (a) Look at your data! ;Email: (b) Know your data! ;Email: (c) Make sure *you* know what a function or tool is doing to your data ;Email: --- ;Email: ;Email: The source netcDF '2_2_ERAI_Prec_2015_Jan_00_12_ 12.nc' was created from a GRIB file: ;Email: ;Email: history attribute: ;Email: 2018-05-16 09:15:25 GMT by grib_to_netcdf-2.7.3: grib_to_netcdf /data/data05/scratch/3c/71/_ma rs-atls00-a562cefde8a29a7288fa 0b8b7f9413f7-b4mPzK.grib -o /data/data04/scratch/41/14/_gr ib2netcdf-atls17-a562cefde8a29 a7288fa0b8b7f9413f7-PZaKeY.nc -utime" ; ;Email: ;Email: --- ;Email: There are 62 time steps on the source ERAI file: ;Email: ;Email: yyyymmddhh tp value ;Email: (0) ERAI: 2015010112 0.000330344 <=== only one value for this calendar day ;Email: ;Email: (1) ERAI: 2015010200 0.0051588 ...... two time steps per calendar day ... ;Email: (2) ERAI: 2015010212 0.00401913 ;Email: (3) ERAI: 2015010300 0 ;Email: (4) ERAI: 2015010312 0 ;Email: (5) ERAI: 2015010400 0 ;Email: (6) ERAI: 2015010412 0.000478998 ;Email: (7) ERAI: 2015010500 0.00344653 ;Email: [SNIP] ;Email: (55) ERAI: 2015012900 0 ;Email: (56) ERAI: 2015012912 0 ;Email: (57) ERAI: 2015013000 9.35942e-05 ;Email: (58) ERAI: 2015013012 7.15703e-05 ;Email: (59) ERAI: 2015013100 0.00390901 ;Email: (60) ERAI: 2015013112 0 ;Email: ;Email: (61) ERAI: 2015020100 0.00183889 <=== only one value for this calendar day ;Email: --- ;Email: To what 'time' does the source ERAI 'time stamp apply"? ;Email: That is ***USER RESPONSIBILITY*** to know or find the answer. ;Email: ;Email: Three possibilities: ;Email: ;Email: [1] beginning ;Email: [2] end ;Email: [3] 'center-of-mass' ... mid temporal time? ;Email: ;Email: --- ;Email: NCL's 'calculate_daily_values' with 'ncrit=2' requires two values on the same calendar day ;Email: https://www.ncl.ucar.edu/ Document/Functions/ Contributed/calculate_daily_ values.shtml ;Email: ;Email: With ncrit=2 the function returns 30 values ... not 31 or 32. Why? ;Email: Because it requires 2 values per calendar date. It ignores the 1st and last times because there is only one time step per calendar day ;Email: ;Email: tp2_ncl ;Email: NCL_tag : calculate_daily_values: arith=sum ;Email: ;Email: ---- ;Email: The function 'tags' it with the 1st time step of each *complete* calendar day ;Email: ;Email: (0) NCL: 2015010200 0.00917792 ;Email: (1) NCL: 2015010300 0 ;Email: (2) NCL: 2015010400 0.000478998 ;Email: (3) NCL: 2015010500 0.00879252 ;Email: [SNIP] ;Email: (28) NCL: 2015013000 0.000165164 ;Email: (29) NCL: 2015013100 0.00390901 ;Email: ;Email: --- ;Email: The CDO takes the first two values [0,1], then the next two [2,3], etc. ;Email: It 'tags' the time as the 2nd time step ;Email: ;Email: (0) CDO: 2015010200 0.00548913 ;Email: (1) CDO: 2015010300 0.00401912 ;Email: (2) CDO: 2015010400 -1.08611e-10 ;Email: [SNIP] ;Email: (29) CDO: 2015013100 0.00398058 ;Email: (30) CDO: 2015020100 0.00183889 ;Email: ;Email: Even though the time 'tags' are the same, they refer to different sums. ;Email: --- ;Email: ;Email: NCL can be made to handle the time and data any way the user wants. ;Email: ;Email: I'll let you look at the attached script and output ;========================================================================= ;E:diri = "/Users/shea/Data/ECMWF/ERAI/" diri = "/home/taghizadeh/97_1_17/5_Era-Interim/" fili = "2_ERAI_Prec_2015_Jan_00_12_12.nc" pthi = diri+fili f = addfile(pthi, "r") tp = short2flt(f->tp) printVarSummary(tp) ; [time | 62] x [latitude | 241] x [longitude | 480] printMinMax(tp,0) print("--------") ymdh_erai = cd_calendar(tp&time, -3) LAT = 10 ; arbitrary grid point LON = 145 print("ERAI: "+ymdh_erai+" "+tp(:,{LAT},{LON})) print("==================================") ; NCL: calculate_daily_values ; default mode: Regardless of the number of values per calendar day: calculate a sum ; There are 32 calendar days represented on the file opt = False tp_ncl = calculate_daily_values(tp, "sum", 0, opt) printVarSummary(tp_ncl) ; [time | 32] x [latitude | 241] x [longitude | 480] printMinMax(tp_ncl,0) print("----") ymdh_ncl = cd_calendar(tp_ncl&time, -3) print("NCL: "+ymdh_ncl+" "+tp_ncl(:,{LAT},{LON})) print("==================================") ; NCL: calculate_daily_values ; set 'nval_crit': require 2-values ****per calendar day**** opt = True opt@nval_crit = 2 ; require 2 or more values per calendar date tp2_ncl = calculate_daily_values(tp, "sum", 0, opt) printVarSummary(tp2_ncl) printMinMax(tp2_ncl,0) print("----") ymdh2_ncl = cd_calendar(tp2_ncl&time, -3) print("NCL: "+ymdh2_ncl+" "+tp2_ncl(:,{LAT},{LON})) print("==================================") ; CDO: Two successive values ;E:F = addfile(diri+"2_2_ERAI_Prec_2015_Jan_24_cdo.nc4", "r") F = addfile(diri+"/cdo_sum/3_ERAI_Prec_2015_Jan_00_12_24_cdo.nc4", "r") tp_cdo = F->tp printVarSummary(tp_cdo) printMinMax(tp_cdo,0) ymdh_cdo = cd_calendar(tp_cdo&time, -3) print("CDO: "+ymdh_cdo+" "+tp_cdo(:,{LAT},{LON})) print("----") ; NCL: Manually program loops ; Two successive values: start at subscript 0 [span two calendar days] dim_tp = dimsizes( tp ) ; ERAI source variable ntim = dim_tp(0) nlat = dim_tp(1) mlon = dim_tp(2) ntim3 = ntim/2 ; 31 tp3_ncl= new( (/ntim3,nlat,mlon/), typeof(tp), getVarFillValue(tp)) ntJump = 2 ; just to be explicit ntStrt = 0 ; " " ntLast = ntim-1 ; " " do nt=ntStrt,ntLast,ntJump tp3_ncl(nt/2,:,:) = tp(nt,:,:) + tp(nt+1,:,:) end do tp3_ncl@long_name = "NCL sum of two successive values regardless of time stamp" tp3_ncl@units = tp@units copy_VarCoords(tp(::2,:,:), tp3_ncl) printVarSummary(tp3_ncl) printMinMax(tp3_ncl,0) print("----") ymdh3_ncl = cd_calendar(tp3_ncl&time, -3) print("NCL_3: "+ymdh3_ncl+" "+tp3_ncl(:,{LAT},{LON})) print("==================================") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;Start: writing tp3_ncl in nc file ; setfileoption("nc", "FileStructure", "Advanced") setfileoption("nc", "Format", "NetCDF4") ;fn = "ncl_wrt_opaque.nc" fn = "3_ERAI_Prec_2017_Jan_00_12_24_ncl.nc4" system("/bin/rm -f " + fn) fo = addfile(fn, "c") fAtt = True ; assign file attributes fAtt@title = "NCL generated netCDF file with enum" fAtt@source_file = fn fAtt@Conventions = "None" fAtt@creation_date = systemfunc("date") ;print(fAtt) fileattdef(fo, fAtt) ;=================================================================== ;nvals = nv(0) ;dimNames = (/"opaque_base_size"/) dimNames = (/"time", "latitude", "longitude"/) ;dimSizes = (/ nvals /) ;dimUnlim = (/ False /) dimUnlim = (/ True , False , False /) ;filedimdef(fo, dimNames, dimSizes, dimUnlim) filedimdef(fo, dimNames, dim_tp, dimUnlim) ;var_name = "opaque_variable" var_name = "tp" ;fileopaquedef(fo, "opaque_type", var_name, var_size, dimNames) fo->$var_name$ = (/tp3_ncl/) ;;;End: writing tp3_ncl in nc file ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; NCL: Two successive values: start at subscript 1 ntim4 = ntim/2-1 tp4_ncl= new( (/ntim4,nlat,mlon/), typeof(tp), getVarFillValue(tp)) ntJump = 2 ; just to be explicit ntStrt = 1 ; " " ntLast = ntim-3 ; " " do nt=ntStrt,ntLast,ntJump tp4_ncl(nt/2,:,:) = tp(nt,:,:) + tp(nt+1,:,:) end do tp4_ncl@long_name = "NCL sum of two successive values regardless of time stamp" tp4_ncl@units = tp@units copy_VarCoords(tp(1:ntim-3:2,:,:), tp4_ncl) printVarSummary(tp4_ncl) printMinMax(tp4_ncl,0) print("----") ymdh4_ncl = cd_calendar(tp4_ncl&time, -3) print("NCL_4: "+ymdh4_ncl+" "+tp4_ncl(:,{LAT},{LON})) print("==================================")