MATLAB计算地理网格面积

等经纬网格对于地学初学者来说是一种极具欺骗性的网格，这是因为地球作为一个球体，纬度越高经圈长度越小，网格面积也就越小。下面提供两种计算方案供参考，欢迎大家讨论

1.不使用WGS84进行计算

这里直接指路Matlab处理气象数据（六）数据的面积加权

2.使用WGS84进行计算

在网格跨度不大的情况下，将网格近似为一个长方形（或者梯形），取网格的中心纬度计算d1（长），取网格的西边界计算d2（高），即有：

function ll_area=m_areat(lat,lon)
    %%M_AREAT 根据输入的lat和lon矩阵计算每个网格的面积
    %当输入数据为等经纬网格时:
    %lat: 1维数组，为每个网格的起始纬度
    %lon: 1维数据，为每个网格的起始经度，0-180代表东经，180-360代表西经
    %ll_area:length(lat)*length(lon)，程序会自动向lat和lon数组末尾按照网格间距补充最后一个网格终止经纬度
    %暂不支持输入数据为曲面网格
    %需要安装Mapping Toolbox 
    % % % % % % % % % % % % % %

    nlat=length(lat);
    res_lat=lat(2)-lat(1);
    lat(nlat+1)=lat(end)+res_lat;
    nlon=length(lon);
    res_lon=lon(2)-lon(1);
    lon(nlon+1)=lon(end)+res_lon;
    ll_area=zeros(nlat,nlon);
    wgs84 = wgs84Ellipsoid("m");
    for i=1:nlat
        for j=1:nlon
            d1=distance((lat(i)+lat(i+1))/2,lon(j),(lat(i)+lat(i+1))/2,lon(j+1),wgs84);       %units:m
            d2=distance(lat(i),lon(j),lat(i+1),lon(j),wgs84);
            ll_area(i,j)=d1*d2/1e6;             %units: km^2
        end
    end
end

经过比较，在低纬地区方案1计算得到的结果较大，在30-40°N地带相差逐渐减小，尔后随着纬度的增大方案2计算得到的结果较大（附比较程序如下）

diff=NaN(80,1);         %save compare result
for k=1:80  
    lat=0+k:0.25:1+k;       %testing from 5 degree N to 65 degree N 
    lon=140:0.25:141;
    ll_area=m_areat(lat,lon);
    csdn_compare=优快云_area(lat+0.25,lon+0.25,0.25);
    diff(k,1)=mean(csdn_compare(:,1)-ll_area(:,1));
end
plot(0:79,diff,'LineWidth',1.2);grid on; 
ylabel('mean diff (km^2)','FontSize',15);
xlabel('latitude \circN','FontSize',15);
yline(0,'color','r','LineWidth',1);

%% function adaped from https://blog.youkuaiyun.com/ymyz1229/article/details/106039683
function output=优快云_area(dLat1,dLon2,res)
    R =  6371.007181;
    dLat2=dLat1-res;
    dLon1=dLon2-res;
    for i=1:length(dLat1)
	    for j=1:length(dLon2)
		    output(i,j) = R*R*abs(sin(dLat1(i)/180.0*pi) - sin(dLat2(i)/180.0*pi)) * abs(dLon1(j) - dLon2(j))*(pi/180.0);
	    end
    end
end