以度为单位的像元面积转换为平方千米（matlab实现）

第一步：

function cellwt = GetAreaWeight(varargin)

if nargin == 1

    deg = varargin{1};

    rows=180/deg;

    cols=360/deg;

    halfSize = deg/2;

    tmp = ones(rows,cols,'single');

    R = makerefmat(-180+halfSize, 90-halfSize, deg, -deg);

    EP = almanac('earth','wgs84','kilometers');

    [~,wt] = areamat((tmp>-1),R,EP);

    cellwt = repmat(wt,[1 cols]);

elseif nargin == 2

    rows=varargin{1};

    cols=varargin{2};

    tmp = ones(rows,cols,'single');

    R = makerefmat(-180+(360/cols)/2, 90-(180/rows)/2, 360/cols, -180/rows);

    EP = almanac('earth','wgs84','kilometers');

    [~,wt] = areamat((tmp>-1),R,EP);

    cellwt = repmat(wt,[1 cols]);

else

    error('not support!');

end

end

注释：利用自定义的area=GetAreaWeight(0.1)函数得到0.1度像元的面积权重。

第二步：

进行仿射变换：R = georasterref('RasterSize', size(area), 'Latlim', [-90 90], 'Lonlim', [-180 180],'ColumnsStartFrom', 'south')。

第三部：

输出为tiff格式：

geotiffwrite('****.tif',area, R)；自定义路径。

第四步：

一个0.08°的像元面积≈(40075.7/360.0*0.08)^2=79.311691(Km2)，然而实际情况并非每个像元面积都相等，像元面积从中纬度到高纬度，表现出由大到小。所以需要将变量与上述转换后的面积相乘，得到实际情况。否则会使得总和的估计偏高。

感谢ZhuZaichun博士提供的代码！