【GDAL基础教程】Python投影栅格并计算栅格面积

【GDAL基础教程】Python投影栅格并计算栅格面积

01 引言：

栅格计算面积相对于矢量而言，不似矢量那么方便，后者可以直接在arcgis等软件，属性表计算几何获取面积。而栅格计算面积的思路通常是从地理坐标系投影到投影坐标系，将栅格由度转为米等单位，最后统计栅格数量乘以单位栅格面积，以获取整个栅格的面积。

现在记录分享下gdal投影栅格并计算面积，记录在此，分享给有需要的同学。

原始栅格数据如下图所示，为2000地理坐标系的sentinel数据。

02 投影代码：

import numpy as np
from osgeo import gdal,osr
import time
​
def projtif(intif,outtif,xres = 10,yres = 10):
    # 创建Albers等面积投影对象
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromProj4("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
    
    gdal.Warp(outtif , # 输出栅格路径
            intif,  # 输入栅格
            format='GTiff',  # 导出tif格式
            dstSRS = spatial_ref, # 投影
            xRes=xres, # 重采样后的x方向分辨率
            yRes=yres, # 重采样后的y方向分辨率
            resampleAlg=gdal.GRA_NearestNeighbour, #最临近重采样
            creationOptions=['COMPRESS=LZW']#lzw压缩栅格
            )
    return outtif​
​
if __name__ == '__main__':
    intif = r'D:\ForestMeteorology\FM230327\data\wheat.tif'
    outtif = r'D:\ForestMeteorology\FM230327\data\wheat_Alberts.tif'
    projtif(intif,outtif,10,10)

03 投影结果：

如下图所示，通过gdal warp成功将地理坐标系的sentinel数据投影为10m分辨率的等面积投影。

04 计算面积：

def tif2arr(tif):
    ds = gdal.Open(tif)
    arr = ds.ReadAsArray()
    return arr
    
arr = tif2arr(outtif)
value,count = np.unique(arr,return_counts=True)
print(arr)
print(value,count)
# 输出结果
[[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]]
[0 1] [341716483  31627211]

通过gdal读取栅格数组，借助np.unique(）获取，栅格唯一值以及唯一值的计数。通过栅格个数以及单位面积栅格，即可获取栅格总面积3162721100平方米。值得注意的是，需要排除栅格的无效值（上述栅格的无效值为0）。

05 完整代码：

# -*- encoding: utf-8 -*-
'''
@File    :   calculate_area.py
@Time    :   2023/03/27 20:50:23
@Author  :   HMX
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''
​
# here put the import lib
import numpy as np
from osgeo import gdal,osr
import time
​
def projtif(intif,outtif,xres = 10,yres = 10):
    # 创建Albers等面积投影对象
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromProj4("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
    
    gdal.Warp(outtif , # 输出栅格路径
            intif,  # 输入栅格
            format='GTiff',  # 导出tif格式
            dstSRS = spatial_ref, # 投影
            xRes=xres, # 重采样后的x方向分辨率
            yRes=yres, # 重采样后的y方向分辨率
            resampleAlg=gdal.GRA_NearestNeighbour, #最临近重采样
            creationOptions=['COMPRESS=LZW']#lzw压缩栅格
            )
    return outtif
​
def tif2arr(tif):
    ds = gdal.Open(tif)
    arr = ds.ReadAsArray()
    return arr
​
​
if __name__ == '__main__':
    t1 = time.time()
    intif = r'D:\ForestMeteorology\FM230327\data\wheat.tif'
    outtif = r'D:\ForestMeteorology\FM230327\data\wheat_Alberts.tif'
    projtif(intif,outtif,10,10)
    arr = tif2arr(outtif)
    value,count = np.unique(arr,return_counts=True)
    print(arr)
    print(value,count)
    areas= [i*10*10 for i in count]
    print(value,areas)
    t2 = time.time()
    print('共计用时{:.2f}'.format(t2-t1))

以上就是本期推文的全部内容了，如果对你有帮助的话，请‘点赞’、‘收藏’，‘关注’，你们的支持是我更新的动力。