Amphan_水稻受损面积说明

1、生成的文件及结果说明

30m空间分辨率的 mDVDI图：

2020_mDVCI_11_win_BGD_19_21_7_rice.tif

细像元（30m)的受损率分布图（每个像素点为该像素点的受损率（0-1) \(P_{fine}\))：

P_fine.tif

水稻像素点（细像元 30m)的受损率分布图(每个像素点的值大于 0表示该像素点为水稻，且值为水稻的受损率（0-1))：

2020_rice_influence_by_fragility_bgd.tif

水稻像素点（细像元 30m)的受损率分布图经过Albers等面积投影投影变换后生成的图：

2020_projected_rice_influence_by_fragility_bgd.tif

统计图鉴上受影响的几个地区的矢量边界图，并经过了Albers等面积投影投影变换

gadm41_BGD_2_projected.shp

统计结果汇总表：

summary.xlsx

统计结果

水稻总面积对比汇总结果如下：

散点图绘制使用的是第二行和第三行的数据，左边这张图是去掉 Jhalokathi 地区的散点图（R2 为 0.91)，右边是去掉 Jhalokathi 和 Pabna 这两个地区的散点图（R2 为 0.92)

（注意这里未使用 Jhalokathi 地区的数据，因为该地区属于沿海地区，分类效果不太好，评估出来的面积与官方数据相差较大；取消 Pabna 是因为官方的统计损失数据感觉有问题，太少了)

说明，

第一行是统计年鉴中4.2.3节的图表中的 Area Under Crops(Acre)这一项的数据（即下面绿色的这一列)

第二行是统计年鉴第3 章预估产量和面积我汇总后的数据

第三行使用2020_BGD_PRED_rice_CUT_mix.tif经过投影变换后统计出来的各地区的水稻种植面积

注意目前看上去第一行和第二行的数据有较大的差异，且第一行的数据看起来问题很大，所以选用第二行的官方统计面积数据进行评价。

水稻受损面积汇总结果：

左边这张图是去掉 Jhalokathi 地区的散点图（R2 为 0.7965)，右边是去掉 Jhalokathi 和 Pabna 这两个地区的散点图（R2 为 0.8622)

说明，

第一行是统计年鉴中4.2.3节的图表中的 Partially damage(Acre)这一项的数据（即下面红色的这一列)

第二行是经过预估后的各地区受损面积（Acre)

2、计算过程

计算水稻受损面积

根据 DVDI 生成 30m mDVDI

输入是

空间分辨率

公式如下：

使用 python 代码进行计算的代码如下（使用 qgis 无法使用 exp 函数，所以用 python 代码处理了)：

紧接着计算高分辨率像素的受损率 \(P_{fine}\)

高分辨率像素的受损率 \(P_{fine}\)的计算公式如下：

生成粗像元农田盖度图

在地理坐标系下粗像素元的分辨率定为0.00898315284230767272度（大概一个像素点是1km)

农田盖度是指农田中植物（如作物）的覆盖面积占总地表面积的比例。它是一个重要的农业生态指标，用于评估农田的生态效益、土壤保持能力以及作物生长状况。

单位

农田盖度的单位通常是百分比（%）。例如，如果一个农田区域的盖度为70%，这意味着该区域的70%被植物覆盖，其余30%为裸地或其他非植被覆盖区域。

输入图层为 30m 左右的水稻分布图：

空间分辨率是：0.0002694945851924635（30m)

其中 1 表示水稻（白色像素)，0 表示非水稻（黑色像素)： 1733987253611

生成 1km 粗像元水稻盖度图只需要直接重投影即可即可，并且设置重采样方法设置为平均即可（取平均的话，每个粗像元像素的值恰好就是农作物的盖度)，并设置好输出分辨率：

生成的（水稻)农田盖度图如下：

栅格直方图统计结果如下：

生成细像元的农田盖度图

因为细像元的空间分辨率就是0.0002694945852404068143度（30 m)

所以直接使用即可（无需重新生成)

计算细像元的受损率 \(P_{fine}\)

公式为： 1734089716820

本来这里还需要再除以一个粗像元的作物盖度，但是发现前面计算\(P_{coarse}\)的时候没乘这个，所以这里看起来不除似乎也不会有问题

主要代码如下（主要思想是先计算出每个粗像元的 mDVDI 加权和和农田盖度总和，即\(P_{fine}\)公式中两个累加的部分，然后计算细像元的受损率\(P_{fine}\))：

生成的细像元的受损率结果图P_fine.tif：

0 为黑色的，1 为白色，可能看不太清楚，下面有 0 为白色的图，1 为黑色的图

对比下粗像元受损率的图(fragility_curve_bgd_1km.tif)如下：

0 为白色的图，1 为黑色的图

计算水稻损失栅格（细像元)

因为已经有了水稻的分布图层（细像元)，1 表示有水稻，0 表示非水稻，同时也有了细像元的受损概率栅格图（细像元)，值为 0-1 的数，表示该像元的受损概率，对应像素相乘就是水稻的损失栅格图

使用 栅格计算器 工具，输入表达式 ("rice_layer" * "damage_probability_layer")，生成一个新的水稻损失栅格图层。这个图层的每个像素值表示该像素点的水稻损失概率。计算结果如下（2020_rice_influence_by_fragility_bgd.tif)：

裁剪受影响的几个地区的栅格边界数据

使用 栅格 > 提取 > 裁剪栅格 工具，将水稻分布栅格和作物受损概率栅格裁剪到各省份的边界范围内

步骤

使用Albers等面积投影进行投影变换

选择适当的等面积投影坐标系（这里选Albers等面积投影）

算了这里直接将整个孟加拉国的水稻损失栅格都进行投影变换了

然后生成要使用的地区的矢量文件并做投影变换

计算每个省份水稻损失总面积：

步骤如下，思路主要是先统计选中的各个地区的矢量边界图层中的各个地区的所有像素点的和（即总的水稻损失概率)，然后乘以每个像素点的面积即为受损的水稻面积

各省份水稻受损率总和（注意还没乘面积)：

Jhalokati:3290.155333611402 Rajbari:26800.093430639397 Chuadanga:48730.667508211176 Jhenaidah: 46816.21448709527 Kustia:72339.49813155603 Narail:22613.452791931308 Pabna: 66233.13581311288 Dinajpur:49931.87775322891

在Albers等面积投影下，一个像素点的面积是 27.8376mx30.0661m,而 1 英亩=4046.8648平方米

水稻受损面积计算结果:

Jhalokati:

3290.155333611402 x27.837x30.0661/4046.8648=680.4515908273128

Rajbari:

26800.093430639397x27.837x30.0661/4046.8648=5542.645972639374

Chuadanga:

48730.667508211176 x27.837x30.0661/4046.8648=10078.20508937573

Jhenaidah:

46816.21448709527x27.837x30.0661/4046.8648=9682.268584349815

Kustia:

72339.49813155603x27.837x30.0661/4046.8648=14960.851872375602

1734078970465 Narail:

22613.452791931308x27.837x30.0661/4046.8648=4676.788286915996

1734079032451 Pabna:

66233.13581311288x27.837x30.0661/4046.8648=13697.968047011705

Dinajpur:

49931.87775322891x27.837x30.0661/4046.8648=10326.632698185073

水稻受损面积分析

首先年鉴中的受损面积如下，这里选择第 3 列作为数据来源，但是去掉了 Pabna 地区的数据（感觉不太真实) 1734328075454

汇总的原始受损数据和预估的受损面积数据如下：

统计水稻种植面积

还是使用分区统计工具进行计算

矢量边界选择要统计几个地区的矢量图形并进行投影变换

栅格图层使用水稻分布图（1 表示水稻，0 表示非水稻)，并记得投影

还是使用总和作为统计信息

计算结果如下：

与前面一样以英亩为单位进行转化：

年鉴中的数据如下：

这个是年鉴中的作物面积和预估的作物面积（好像有点对不上，Pabna先不管）

数据校验

1、水稻品种问题

突然发现孟加拉国的水稻是有三个品种的，不过这三个品种的水稻的种植期是不一样的

1734512348388 1734184652087

3、再仔细统计一下2020 年年检中这几个地区的数据：

先水稻有这三个品种，且三个品种还有两到三个品种

这三个品种的日期如下：

在农业统计中，像“2020-21”这样的时间段通常表示一个跨年度的农业年度或财政年度。具体到孟加拉国的数据，“2020-21”通常指的是从2020年7月1日到2021年6月30日的时间段

所以我们选用 2019-2020 的数据

Jhalokati:

Aus:26600英亩

Aman：118050 英亩（Aman 应该用不到)

Boro：22542 英亩

Rajbari:

Aus：4493 英亩

Boro：32406 英亩

Chuadanga:

Aus：96922 英亩

Boro：77565 英亩

Jhenaidah:

Aus：75694 英亩

Boro：200233 英亩

Kustia:

Aus：66222 英亩

Boro：84391 英亩

Narail:

Aus：15868 英亩

Boro：121723 英亩

Pabna:

Aus: 46487英亩

Boro：128094 英亩

Dinajpur:

Aus：17023 英亩

Boro：426619 英亩

整理数据如下：

突然发现代码写错了，计算 P_fine的时候括号放的地方放错了

就下面红框这

完蛋，得出来的受损率好像偏高了（0 是黑的，越亮表示受损率越高）

对比下 P_coarse,大致趋势倒是还是有的

现在来看看具体细节

这个是水稻的分布，黑色代表水稻

P_fine如下：

仔细看下大部分地方是无水稻的地方变成了黑色，而有水稻的地方变为了浅灰色，但是还是有很多地方异常值

上面有的像元是斜的，这是 mDVDI 引入的，DVDI 的计算使用的 modis 数据集，重采样到 30m 的时候需要加个影像投影和重采样的过程，这个是我之前写的代码https://code.earthengine.google.com/e870be706993fa650f75670aa70d13ce，参考的是https://blog.csdn.net/qq_35591253/article/details/127528788

先来分析下异常的地方