01 引言：
栅格计算面积相对于矢量而言，不似矢量那么方便，后者可以直接在arcgis等软件，属性表计算几何获取面积。而栅格计算面积的思路通常是从地理坐标系投影到投影坐标系，将栅格由度转为米等单位，最后统计栅格数量乘以单位栅格面积，以获取整个栅格的面积。
现在记录分享下gdal投影栅格并计算面积，记录在此，分享给有需要的同学。
原始栅格数据如下图所示，为2000地理坐标系的sentinel数据。
02 投影代码：
import numpy as npfrom osgeo import gdal,osrimport timedef 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 outtifif __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 arrarr = tif2arr(outtif)= 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]]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 libimport numpy as npfrom osgeo import gdal,osrimport timedef 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 outtifdef tif2arr(tif):ds = gdal.Open(tif)arr = ds.ReadAsArray()return arrif __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}s'.format(t2-t1))
戳这里关注我
请点赞、在看、关注，你们的支持是我更新的动力。
