交换的学期结束了，自己也从香港飞回了上海。隔离期间自己在并行多线程处理许多数据，其中就包括计算亚洲2020年的耕地面积（数据的来源为 Globeland）。下面的部分是以东亚为例，提供了一个简单的介绍。

获取“东亚”范围的 shapefile

根据国际上的定义，“东亚”共包含若干个地区。GDAM 数据库包含了全球的行政边界，因此可以从 GADM 数据库内依次下载各地区的行政边界，而后在 Arcmap 中将其合并。需要注意的是，GADM 数据库中的部分地区与中国政府的国界主张存在冲突，因此在使用时需要格外注意。

合并后的行政边界的投影为默认的 WGS 84, 而 Globeland 中的数据则是按 UTM 投影的 6 度带分幅存储。从GlobeLand的官网上查看研究区域：

分幅统计耕地面积

因为各幅遥感影像已经投影在所处的 6 度带内，为谨慎起见，我不打算将它们合并，而是分幅统计计算。处理的思路如下：

1. 将 shapefile 文件投影到栅格所处的投影中；
2. 相同投影下，将两者叠加（即按掩膜提取）；
3. 统计叠加后栅格内耕地像元的数目；

其中前两步都很基础，代码如下：

import arcpy
import glob

result = 0
raster_list = []
temp_list = glob.glob(r"D:\Desktop\land_use\2020 global total landuse part2\n*.tif")

for i in temp_list:
    if int(i.split("\\")[-1][1:3]) in range(43,56,1):  # east asia extent
        raster_list.append(i)
    
for raster in raster_list:
    # get raster projection
    out_projection = arcpy.Describe(raster).spatialReference

    # 1. reproject shp file
    num = raster.split("\\")[-1][0:3]
    input_shp = r"D:\Desktop\land_use\land_sum\East_Asia.shp"
    output_shp = r"D:\Desktop\land_use\temp_ea/" + num + ".shp"
    arcpy.Project_management(input_shp, output_shp, out_projection)

    # 2. extract by mask and export raster
    result = arcpy.sa.ExtractByMask(raster, output_shp)
    result.save(r"D:\Desktop\land_use\temp_ea/" + num + "_clipped.tif")

而第三步值得进一步把玩，我在这里取了个巧：

读取栅格属性表字段

在执行完按掩膜提取后，裁剪后的栅格在当前投影下自行统计了本幅影像的像元个数，并将结果存储在了属性表中。因此，可以通过读取属性表中“耕地”一项值的个数，并将其乘以像元大小（也就是30*30），来快速计算该幅影像中耕地的面积。

据查阅资料发现，arcpy 中并没有直接、易用的适用于提取栅格属性表字段的函数。但是，栅格的属性表是存储在 .dbf文件中的，因此一个很好的迂回方法是绕过 arcpy，直接读取.dbf文件，获得字段并进行计算。这里需要导入dbfread库，而这个库是自带于 arcpy 包中的，因此可以直接使用：

代码可以类似这样写：

from dbfread import DBF
import pandas as pd

result = 0
for raster in raster_list:
    # some codes......
    table = DBF(r"D:\Desktop\land_use\temp_ea/" + num + "_clipped.tif.vat.dbf")
    df = pd.DataFrame(iter(table))
    arable_area = df.loc[1]['Count'] * 30 * 30  # important!!! count*size
    print(raster.split("\\")[-1] + ' area: ' + str(arable_area))
    result += arable_area

计算结果

根据上述代码，最终批量生成了各个栅格。

最左端 N43 投影带：

最右端 N55 投影带：

计算结果如下：

汇总后，求得的总面积为225 * 10^4 km^2.

反思与不足

在刚开始想到用这种方法计算的时候，我还很兴奋，觉得这是个绝佳的方法。不过回过神定心想了一下，在对比 10 年的数据，发现这个值其实是偏大的，原因在于虽然在特定投影下面积的计算是准确的，但两幅相邻的栅格在某些情况下存在重叠区域，而按上述方法计算则没有去除重叠的这些部分。

现在自己想到了一种迭代的方法逐个计算，将之前计算过的区域给剔除掉。后续再写。

另外，自己已经很久没有定下心写博客了。“我们走了很久，却不知为何出发”，一年到头，自己也该收拾下思绪，好好总结、复盘、忏悔这一年了。