0.5°的nc数据，想要计算每个grid的面积，求讨论

sunzy025

我的数据是nc格式的，有0.5*0.5（全球720*360）的，也有其他分辨率的数据但都是以°为单位，我想要计算每一个格网的面积。

https://gis.stackexchange.com/qu ... ate-area-of-rasters      Eq.1
这个帖子有介绍如何计算以°为单位的raster的数据的面积的。


我的指导教师和我讨论后得到计算公式为：

grid_area(km2) = 111*111*degree*degree*cos(latitude*(3.14159274/180))     Eq.2

我有点理解不理解Eq.2的具体含义，还望大家伙帮帮忙，指点一二，感恩！

sunzy025

学习渣渣 发表于 2019-4-20 14:51
百度。曲面多边形面积。

嗯嗯！好主意，我后来试验了，确实精度还是可以接受的，感谢

sunzy025

文质彬 发表于 2019-4-17 15:47
每个1*1°格点的面积等于所在格点的dlat*dlon
dlat-------大约111km，固定值
dlon------会随纬度变化，dl ...

后来，我就这么算的了，感觉比111km这个准，但是是UTM投影

a = 6378137;
b = 6356752.3142;
e = sqrt(1 - (b/a)^2);
f1 = center_lat + pixel_y./2;
f0 = center_lat - pixel_y./2;
zm0 = 1 - e.*sind(f0);zm1 = 1 - e.*sind(f1);
zp0 = 1 + e.*sind(f0);zp1 = 1 + e.*sind(f1);
area0 = pi .* b.^2 .* (log(zp0./zm0) ./ (2.*e) + sind(f0) ./ (zp0.*zm0));
area1 = pi .* b.^2 .* (log(zp1./zm1) ./ (2.*e) + sind(f1) ./ (zp1.*zm1));
area = (area1 - area0) .* pixel_x ./ 360;
area = repmat( area , 1 , length(x) );
area = area';

sunzy025

虫儿飞 发表于 2019-4-17 10:49
我做过面积加权的格点平均降水，精度还可以的，总体的思路就是求梯形面积，如果考虑地形，那就考虑每一个格 ...

后来，我这么算了
a = 6378137;
b = 6356752.3142;
e = sqrt(1 - (b/a)^2);
f1 = center_lat + pixel_y./2;
f0 = center_lat - pixel_y./2;
zm0 = 1 - e.*sind(f0);zm1 = 1 - e.*sind(f1);
zp0 = 1 + e.*sind(f0);zp1 = 1 + e.*sind(f1);
area0 = pi .* b.^2 .* (log(zp0./zm0) ./ (2.*e) + sind(f0) ./ (zp0.*zm0));
area1 = pi .* b.^2 .* (log(zp1./zm1) ./ (2.*e) + sind(f1) ./ (zp1.*zm1));
area = (area1 - area0) .* pixel_x ./ 360;
area = repmat( area , 1 , length(x) );
area = area';

感觉还不错呢

学习渣渣

百度。曲面多边形面积。

文质彬

每个1*1°格点的面积等于所在格点的dlat*dlon
dlat-------大约111km，固定值
dlon------会随纬度变化，dlat=111*cos（（lat/180）*3.14) ：角度单位——度化π
那么0.5*0.5的面积就是1*1°的面积*0.5*0.5

虫儿飞

我做过面积加权的格点平均降水，精度还可以的，总体的思路就是求梯形面积，如果考虑地形，那就考虑每一个格点的地形包络，这个我没有研究过，但是思路都是类似的：在每一个格点内，利用高分辨率的高程和地形数据，把格点再次分割细小网格，用求解梯形面积的方法，求解每一个细小网格的面积