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:
>
> 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
>
> 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
>
> ; 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.
>
> 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

begin
;************************************************
; ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface_gauss/
;************************************************
diri = "/Users/shea/Data/CDC/"
fili = "air.2m.gauss.1979.nc"

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) )

;************************************************
; Miscellaneous
;************************************************
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_cfd1 = TG_cfd1*scale
TG_cfd2 = TG_cfd2*scale
rmse1 = rmse1*scale
rmse2 = rmse2*scale

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

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_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