搜 索

基于Globeland统计2020年东亚的耕地面积(上)

  • 768阅读
  • 2021年12月28日
  • 0评论
首页 / 学习区 / 正文

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

获取“东亚”范围的 shapefile

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

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

east_asia_area

分幅统计耕地面积

因为各幅遥感影像已经投影在所处的 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 投影带:

N43

最右端 N55 投影带:

N55计算结果如下:

calculation

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

反思与不足

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

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

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

评论区
暂无评论
avatar