cdtgradientdocumentation
cdtgradientcalculates the spatial gradient of gridded data equally spaced in geographic coordinates.
Back to Climate Data Tools Contents.
Contents
Syntax
[FX,FY] = cdtgradient(lat,lon,F) [FX,FY] = cdtgradient(lat,lon,F,'km')
Description
[FX,FY] = cdtgradient(lat,lon,F)for the gridded variable F and corresponding geographic coordinateslatandlon,cdtgradientcalculatesFX, the spatial rate of west-to-east change inFper meter along the Earth's surface, andFY, the south-to-north change inFper meter. This function assumes an ellipsoidal Earth as modeled by thecdtdimandearth_radius. A positive value ofFXindicates thatF增加在网格单元,从西到东a positive value ofFYindicatesFincreases from south to north.Fcan be a 2D or 3D matrix whose first two dimensions must correspond tolatandlon. IfFis 3D, outputsFXandFYwill also be 3D, with each grid along the third dimension calculated separately.
[FX,FY] = cdtgradient(lat,lon,F,'km')returns gradients per kilometer rather than the default meters.
Example 1: Theory
Here's a global grid of some variableF, which is only a function of latitude and increases from south to north at a rate of 1 unitFper degree of latitude. Usecdtgridto create the grid and then defineFas the value of the latitude plus 100:
% Create a quarter-degree grid:[lat,lon] = cdtgrid(0.25);% Define F:F = lat + 100;% Plot F on a globe:figure globepcolor(lat,lon,F) globeborders('color',rgb('gray')) axistightcb = colorbar; ylabel(cb,'data F')
In the figure above we see that the range ofFvalues go from about 10 near the South Pole to about 190 near the North Pole. That's becauseFis just the latitude of each grid cell plus 100.
With this dataset we know that at any given latitude,Fis the same for all longitudes. That meansFxshould be zero everywhere on the globe, becauseFnever changes from west to east.
However,Fdoes change from south to north, at a rate of 1 unitFper degree latitude. It's handy to know that degrees of latitude are separated by about 111 km (from the original definition of the meter, which said the distance from the Equator to the North Pole is 10 million meters). So if each degree of latitude is separated by 111 km, andFincreases by 1 unit for each degree of latitude, thenFyshould be 1/111 = 0.009 everwhere on the globe. Let's see:
[Fx,Fy] = cdtgradient(lat,lon,F,'km'); figure subplot(1,2,1) pcolor(lon,lat,Fx) borders('countries','color',rgb('gray')) shadinginterpaxisoff% removes tickscb1=colorbar('location','southoutside'); xlabel(cb1,'\partialF/\partialx (units of F)/km') subplot(1,2,2) pcolor(lon,lat,Fy) shadinginterp边界('countries','color',rgb('gray')) axisoff% removes tickscb2=colorbar('location','southoutside'); xlabel(cb2,'\partialF/\partialy (units of F)/km')
Above, we see that our predictions fared pretty well.Fxis zero everywhere in the world, as we predicted, and the values ofFycenter on our predicted value of 9e-3. There is some latitudinal variation inFy, however, because the Earth is not a perfect sphere, but an ellipsoid whose degrees of latitude are not perfectly spaced by 111,111 m, as the original definition of the meter might have you believe.
Example 2: Reality
For this example, we'll use the example surface pressure data that comes with CDT. Start by loading it:
文件名='ERA_Interim_2017.nc'; sp = ncread(filename,'sp');% surface pressurelat = double(ncread(filename,'latitude')); lon = double(ncread(filename,'longitude')); [Lat,Lon] = meshgrid(lat,lon);
The absolute surface pressure isn't particularly interesting because it mostly tracks surface elevation. But surface pressureanomaliesare quite a bit more interesting because they give us an idea of how far out of balance the system is at a given point in time.
For convenience, let's calculate the surface pressure anomaly for January of 2017 as the surface pressure field for that month, minus the mean surface pressure from all of 2017:
% Surface pressure "anomaly":温泉= sp(:: 1)意味着(sp, 3);
Here's the surface pressure anomaly we want to analyze:
figure globepcolor(Lat,Lon,spa); globeborders('color',rgb('gray')) axistightcmocean('delta','pivot')% sets colormapcb = colorbar; ylabel(cb,'surface pressure anomaly (Pa)') view(125,5)% sets viewing angle
Calculate the zonal and meridional surface pressure gradients in Pa/km like this
[Sx,Sy] = cdtgradient(Lat,Lon,spa,'km'); figure subplot(1,2,1) pcolor(Lon,Lat,Sx) borders('countries','color',rgb('gray')) shadinginterpaxisoff% removes tickscb1=colorbar('location','southoutside'); xlabel(cb1,'zonal pressure gradient Pa/km') caxis([-1 1]) cmoceandiffsubplot(1,2,2) pcolor(Lon,Lat,Sy) shadinginterp边界('countries','color',rgb('gray')) axisoff% removes tickscb2=colorbar('location','southoutside'); xlabel(cb2,'meridional pressure gradient /km') caxis([-1 1]) cmoceandiff
The figure above may be somewhat difficult to interpret, so we'll take a different approach. Instead of plottingSxandSyas color, plot them as vectors atop the surface pressure anomaly map:
figure pcolor(Lon,Lat,spa) shadinginterpholdon边界('countries','color',rgb('gray'),'center'，180) cmocean('delta','pivot')% sets colormapcb = colorbar; ylabel(cb,'surface pressure anomaly (Pa)') quiversc(Lon,Lat,Sx,Sy,'k','density'，100)
In the map above, we see that gradient vectors always point from low values to high values. Naturally, that's the opposite direction of the pressure gradient force.
Author Info
This function is part of theClimate Data Toolbox for Matlab. The function and supporting documentation were written by Chad A. Greene of the University of Texas at Austin.