【更新】根据Google地球测算中国国土面积

点击文件名下载附件



附程序

function [area_result]=area_sphere(point_list_raw)

% Calculate the area inside a loop defined by the 'point_list_raw' on a elliptical
% sphere. The 'point_list_raw' contains the theta and phi coordinate.

N=2000;

earth_radius_long=6378.1370;
earth_radius_short=6356.7523;

theta=90-point_list_raw(:,2);
phi=point_list_raw(:,1);

[row,col]=size(point_list_raw);

theta=theta/180*pi;
phi=phi/180*pi;

theta_min=min(theta);
theta_max=max(theta);
theta_span=theta_max-theta_min;
theta_dalta=theta_span/N;

phi_min=min(phi);
phi_max=max(phi);
phi_span=phi_max-phi_min;
phi_dalta=phi_span/N;

area_result=0.;
for i=1:N
    theta_current=theta_min+(i-1+0.5)*theta_dalta;
    sum_theta_current=0.;
    for j=1:(row-1)
        cross_tag=(theta(j)-theta_current)*(theta(j+1)-theta_current);
        dif_theta=abs(theta(j)-theta(j+1));
        if (cross_tag<0. & dif_theta<10/180*pi)
            phi_current=(phi(j)+phi(j+1))/2;
            if (theta(j)<theta_current)
                sum_theta_current=sum_theta_current+phi_current;
            else
                sum_theta_current=sum_theta_current-phi_current;
            end
        end
    end
    ds=(earth_radius_long*sin(theta_current))^2+(earth_radius_short*cos(theta_current))^2;
    area_result=area_result+sin(theta_current)*sum_theta_current*theta_dalta*ds;
end




之前的计算把纬度和theta方向弄混了，他们相差90度角。新的计算并且考虑了椭圆球面的影响。从新的网站可以下载更精确的边界。

最近大家为了中国国土大小争吵不休，我使用google地球估算了一下中国国土面积。

方法如下：

一、 获取中国边界

首先从网站下载边界的数据文件。我使用的是 http://www.gadm.org/countryres 网站的信息。获得各国的边界的kml文件（见附件world_Countries.kml）。在这个网站中，中国主要由3部分组成：中国大陆，海南岛，其他岛屿。其他岛屿太小忽略不计。另外台湾列为了独立国家。这样需要提取中国大陆，海南岛，台湾岛的边界（见附件Mainland.kml，Hainan.kml，Taiwan.kml）。 把kml文件中的坐标复制出来放入一个文本文件中（我使用dat做为后缀）（mainland.dat，hainan.dat，taiwan.dat）。

二、 计算面积

Google地图提供的是经纬度坐标，相当于球坐标中的theta和phi分量。利用球坐标的积分公式，可以计算出边界包围的面积。我使用Matlab计算（程序为area_sphere.m）。把刚才得到的dat文件load进matlab就可以调用area_sphere.m求解出椭圆圆球上围成的面积。地球的长短半轴分别取6378.1370和6356.7523公里。

三、 结果

根据这一计算，中国大陆（含内陆湖泊，不含藏南，阿克赛钦）的面积是933.81万平方千米，海南和台湾分别是3.4198和3.6096万平方千米。

四、 精度

从图中可以看到这个网站提供的边界和google提供的边界比较接近，中国大陆的边界比较精细但是海南和台湾的边界比较粗糙。另外可以计算中国的边界长度（包括海洋和陆地）是3.2409e+004千米，如果认为google地图边界误差不超过10千米，地球半径平方误差为<1/150，可以根据此估计误差大小为40万平方公里。如果各位有更精细的地图文件，可以使用这一程序进行技术。

同时计算了美国的面积（含五大湖）为973.43万平方千米。

中国国土面积计算.rar(943 KB, 下载次数: 62)2011-10-19 12:15 上传点击文件名下载附件



附程序

function [area_result]=area_sphere(point_list_raw)

% Calculate the area inside a loop defined by the 'point_list_raw' on a elliptical
% sphere. The 'point_list_raw' contains the theta and phi coordinate.

N=2000;

earth_radius_long=6378.1370;
earth_radius_short=6356.7523;

theta=90-point_list_raw(:,2);
phi=point_list_raw(:,1);

[row,col]=size(point_list_raw);

theta=theta/180*pi;
phi=phi/180*pi;

theta_min=min(theta);
theta_max=max(theta);
theta_span=theta_max-theta_min;
theta_dalta=theta_span/N;

phi_min=min(phi);
phi_max=max(phi);
phi_span=phi_max-phi_min;
phi_dalta=phi_span/N;

area_result=0.;
for i=1:N
    theta_current=theta_min+(i-1+0.5)*theta_dalta;
    sum_theta_current=0.;
    for j=1:(row-1)
        cross_tag=(theta(j)-theta_current)*(theta(j+1)-theta_current);
        dif_theta=abs(theta(j)-theta(j+1));
        if (cross_tag<0. & dif_theta<10/180*pi)
            phi_current=(phi(j)+phi(j+1))/2;
            if (theta(j)<theta_current)
                sum_theta_current=sum_theta_current+phi_current;
            else
                sum_theta_current=sum_theta_current-phi_current;
            end
        end
    end
    ds=(earth_radius_long*sin(theta_current))^2+(earth_radius_short*cos(theta_current))^2;
    area_result=area_result+sin(theta_current)*sum_theta_current*theta_dalta*ds;
end