Re: gradient computation

Thanks a lot for the support, Dennis.

Best regards,

Mateus

2009/7/7 Dennis Shea <shea_at_ucar.edu>

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

-- 
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
Received on Tue Jul 07 2009 - 15:21:10 MDT