Re: gradient computation

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue, 07 Jul 2009 14:58:18 -0600

I just took a look at this.
I modified an existing test script.

[1] The gradients of a gaussian grid were computed
     via spherical harmonics [gradsg]. This was considered "truth".

[2] I computed the gradients using 'center_finite_diff'
     as outlined below.

Result: The root-mean-square difference of [2] from that
     calculated by [1] was 3.94e-06 for both.

     For all intents and purposes ... That is machine accuracy.

Test script attached.

Cheers
D

Mateus Teixeira wrote:
> Hi Dennis,
>
> I think making the relationship 1deg -> 111.123e3 m for all latitudes
> could result in errors away from 45 deg lat.
> I made some tests comparing the meridional gradient of temperature with
> two sets of data, both from reanalysis. Please, see below:
>
> ; read in data
> f = addfile( "air.1979.nc <http://air.1979.nc>", "r"
> ) ; data with dlat fixed
> ;f = addfile( "air.2m.gauss.1979.nc <http://air.2m.gauss.1979.nc>", "r"
> ) ; data with dlat not fixed
> T = short2flt( f->air(0,0,:,:) ) ;
> temperature at lowest level and for the first time
> lat = f->lat
> lon = f->lon
>
> ; deg to rad
> g2r = acos(-1.)/180
>
> ; for data with dlon and dlat fixed, equal to 2.5 degrees
> ; for data with dlat not fixed, the latitude difference used can
> introduces errors
> ; first way to compute the gradient => directly by center_finite_diff
> dT_dy = center_finite_diff( T(lon|:,lat|:), (lat(2)-lat(1))*g2r, False,
> 0 ) * 6.37e-6
>
> ; second way to compute gradient => computing variations of temperature
> and latitude in y direction
> dTy = center_finite_diff( T(lon|:,lat|:), 1., False, 0 ) ;
> variation of temperature in y direction
> dy = center_finite_diff( lat, 1., False, 0 ) *g2r ;
> variation of latitude in y direction (in radians)
> dyND = conform( T(lon|:,lat|:), dy, 1 ) ;
> same dimensions of dTy
> dTy_dyND = (dTy / dyND) * 6.37e-6 ; the
> GRADIENT
>
> ; copy coodinates
> copy_VarCoords( T(lon|:,lat|:), dTy_dyND )
> copy_VarCoords( T(lon|:,lat|:), dT_dy )
>
> Ploting 'dT_dy' and 'dTy_dyND' should result in the same graphics, since
> dlat = 2.5, it's fixed.
>
> But, when dlat isn't fixed, the graphics will have very little
> differences (this differences can be seen by ploting individual grid
> point values), but I think if the gradient will be used in other
> computations, the error can be propagated. So, for non-fixed dlat, the
> second way to compute gradient will be better.
>
> Best regards,
>
> Mateus
>
>
>
>
>
>
>
>
>
> 2009/7/6 Dennis Shea <shea_at_ucar.edu <mailto:shea_at_ucar.edu>>
>
> I must leave but this might be more clear
>
> ; compute K/deg_of_lat
> dT_dLat = center_finite_diff( T( lon|:, lat|: ), lat, False, 0 )
>
> ; convert to K/m: one deg lat @ 45N is approx 111.123e3 meters
> dT_dLat = dT_dLat/111.123e3
> dT_dLat_at_units = "K/m"
>
> Of course, this could be done in one statement.
>
> If this works, please post to the ncl-talk question.
>
> THX
> D
> >
>
> Mateus Teixeira wrote:
>
> Dennis,
>
> Couldn't it be as below
>
> dy = center_finite_diff( lat, 1., False, 0 ) * 6.37e6 * acos(
> -1. )/180.
> dyND = conform( T, dy, 1 )
> dTy = center_finite_diff( T( lon|:, lat|: ), 1., False, 0 )
>
> dTy_dyND = dTy / dyND
>
>
> Mateus
>
>
>
> 2009/7/6 Dennis Shea <shea_at_ucar.edu <mailto:shea_at_ucar.edu>
> <mailto:shea_at_ucar.edu <mailto:shea_at_ucar.edu>>>
>
>
> What I sent you was incorrect.
> I will have to think about this.
>
> D
>
> Mateus Teixeira wrote:
>
> Dear NCL users,
>
> I have a temperature field that is not equally spaced in
> y, but
> it is in x. I'm not sure how to compute its gradient in y
> direction. I'm wondering if I can apply this
>
> dphi = center_finite_diff( lat, 1., False, 0 ) * 6.37e6 *
> acos(
> -1. )/180. dT_dy = center_finite_diff( T( lon|:, lat|: ),
> dphi, False, 0 )
>
> or I just use
>
> dT_dy = center_finite_diff( T( lon|:, lat|: ), lat *
> acos(-1.)/180., False, 0 ) / 6.37e6
>
> Thanks.
>
>
> Best regards,
>
> -- Mateus da Silva Teixeira
>
> Registered Linux User #466740
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
>
> --
> Mateus da Silva Teixeira
>
> Registered Linux User #466740
>
>
>
>
> --
> Mateus da Silva Teixeira
>
> Registered Linux User #466740

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"
begin
;************************************************
; Read gaussian array
; ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/
;************************************************
  diri = "/Users/shea/Data/CDC/"
  fili = "air.2m.gauss.1979.nc"
  f = addfile(diri+fili,"r")

  TG = short2flt( f->air(0,:,:) ) ; only one time for test
  TG_at_long_name = "T-GausGrid" ; plot labeling clarity
  printVarSummary(TG) ; (94,192) gaussian [N-S]
  printMinMax(TG,True)

  dimTG= dimsizes(TG)
  nlat = dimTG(0) ; 94
  mlon = dimTG(1) ; 192
  
;************************************************
; Calculate gradients via spherical harmonics
;************************************************
  TG = TG(::-1,:) ; spherical harmonics want S-N

  TG_gradsg_lon = new( dimTG, typeof(TG), getFillValue(TG) )
  TG_gradsg_lat = new( dimTG, typeof(TG), getFillValue(TG) )
  gradsg (TG, TG_gradsg_lon, TG_gradsg_lat)
                                         ; add meta data
  copy_VarCoords(TG, TG_gradsg_lon)
  copy_VarCoords(TG, TG_gradsg_lat)
  TG_gradsg_lon_at_long_name = TG_at_long_name + ": Zonal Gradient"
  TG_gradsg_lon_at_units = "K/m"
  TG_gradsg_lat_at_long_name = TG_at_long_name + ": Meridional Gradient"
  TG_gradsg_lat_at_units = "K/m"

  printMinMax(TG_gradsg_lat, True)
  printMinMax(TG_gradsg_lon,False)
  
;************************************************
; Miscellaneous
;************************************************
  rad = 4.*atan(1.)/180.
  lat = TG&lat
  TT = TG(lon|:,lat|:) ; for convenience

  re = 6.37122e6 ; spherical earth
  con = re*rad ; one deg lat = 111198.8 meters
                               
  scale = 1e6 ; scale => nicer plots
;************************************************
; Use 'center_finite_diff' for meridional gradient
;************************************************

; METHOD 1 [simple]
  work = center_finite_diff( TT, lat, False, 0 ) ; K/deg_of_lat
  work = work/con ; K/m
  copy_VarCoords( TT, work) ; coord variables

  TG_cfd1 = work(lat|:,lon|:) ; reorder for plot
  TG_cfd1_at_long_name = "TG_cfd1: meridional"
  TG_cfd1_at_units = "K/m"
  printVarSummary(TG_cfd1)
  printMinMax(TG_cfd1,False)
  delete(work)

; METHOD 2 [Mateus]
  dlat = center_finite_diff( lat, 1., False, 0 )*rad
  dTT = center_finite_diff( TT, 1., False, 0 )
  work = (dTT/conform(dTT, dlat, 1))/re

  copy_VarCoords( TT, work) ; meta
  TG_cfd2 = work(lat|:,lon|:) ; reorder for plot
  TG_cfd2_at_long_name = "TG_cfd2: meridional"
  TG_cfd2_at_units = "K/m"
  printVarSummary(TG_cfd2)
  printMinMax(TG_cfd2,False)
  delete(work)

;************************************************
; Calculate RMS difference
;************************************************
  gwt = latGauWgt(nlat, "lat", "gaussian weights", "")

  rmse1 = wgt_arearmse(TG_gradsg_lat, TG_cfd1, gwt, 1.0, 0)
  rmse2 = wgt_arearmse(TG_gradsg_lat, TG_cfd2, gwt, 1.0, 0)
  print(" ")
  print("rmse1="+rmse1+" rmse2="+rmse2)

;************************************************
; create default plot
;************************************************
  wks = gsn_open_wks("ps","tst_grad") ; open a ps file
  gsn_define_colormap (wks, "ViBlGrWhYeOrRe")
  plot = new(3,graphic) ; create a plot array
  
  res = True
  res_at_mpFillOn = False
  plot(0) = gsn_csm_contour_map_ce(wks,TG,res)

  res_at_gsnDraw = False ; don't draw
  res_at_gsnFrame = False ; don't advance frame
  res_at_gsnSpreadColors= True ; use full range of color map
  res_at_cnFillOn = True
  res_at_cnLinesOn = False
  res_at_cnLineLabelsOn = False
  res_at_lbLabelBarOn = False ; turn off individual lb's

 ;res_at_cnInfoLabelOrthogonalPosF = -0.07 ; move the label inside th plot
  res_at_lbOrientation = "vertical" ; vertical label bar

  TG_gradsg_lat = TG_gradsg_lat*scale
  TG_gradsg_lon = TG_gradsg_lon*scale
  TG_cfd1 = TG_cfd1*scale
  TG_cfd2 = TG_cfd2*scale
  rmse1 = rmse1*scale
  rmse2 = rmse2*scale

  symMinMaxPlt (TG_gradsg_lat,20,False,res)
  plot(0) = gsn_csm_contour_map_ce(wks,TG_gradsg_lat,res)

  symMinMaxPlt (TG_gradsg_lon,20,False,res)
  plot(1) = gsn_csm_contour_map_ce(wks,TG_gradsg_lon,res)

  resP = True
  resP_at_gsnMaximize = True
  gsn_panel(wks,plot(0:1),(/2,1/),resP)

  symMinMaxPlt (TG_gradsg_lat,20,False,res)
  plot(0) = gsn_csm_contour_map_ce(wks,TG_gradsg_lat,res)

  res_at_gsnCenterString = "rmse1: "+sprintf("%5.2f", rmse1) \
                      + " [*"+sprintf("%5.4g",(1/scale))+"]"
  plot(1) = gsn_csm_contour_map_ce(wks,TG_cfd1 ,res)

  res_at_gsnCenterString = "rmse2: "+sprintf("%7.2f", rmse2) \
                      + " [*"+sprintf("%5.4g",(1/scale))+"]"
  plot(2) = gsn_csm_contour_map_ce(wks,TG_cfd2 ,res)

 ;resP_at_txString = "Gradient Test"
  resP_at_gsnPanelLabelBar = True ; add common colorbar
  resP_at_lbLabelFontHeightF = 0.007 ; make labels smaller
  resP_at_lbLabelAutoStride = True
  resP_at_lbTitleOn = True ; turn on title
  resP_at_lbTitleString = "[*"+sprintf("%5.4g",(1/scale))+"]"
  resP_at_lbTitlePosition = "Right" ; title position
  resP_at_lbTitleFontHeightF = .0125 ; make title smaller
  resP_at_lbTitleDirection = "Across" ; title direction
  resP_at_pmLabelBarParallelPosF = 0.05

  gsn_panel(wks,plot,(/3,1/),resP)
end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Jul 07 2009 - 14:58:18 MDT

This archive was generated by hypermail 2.2.0 : Wed Jul 08 2009 - 14:48:16 MDT