python gdal实现栅格数据的固定步长重分类

这里以DEM为例，分成间隔为50米的类别。

from osgeo import gdal,osr
import numpy as np
import pandas as pd
import bisect

demfile = 'G:/data/dem/dem.tif' #栅格数据文件路径
dem = gdal.Open(demfile) #GDAL库打开栅格文件
projection = dem.GetProjection() #得到栅格数据的投影
geotransform = dem.GetGeoTransform() #得到仿射变换，计算坐标偏移
im_width = dem.RasterXSize #栅格矩阵的列数
im_height = dem.RasterYSize #栅格矩阵的行数
dem_array = dem.GetRasterBand(1).ReadAsArray().astype(int) #读取栅格的第一波段，读取后以数组形式存储，类型为整型
dem_mask=np.ma.masked_where(dem_array==-2147483647.0,dem_array)#将背景值(NoData)为-2147483647.0的掩膜掉
rows, columns = dem_mask.shape
dem_reclass = np.full_like(dem_mask, 0.00, dtype=float)#存放重分类值
dem_max = np.max(dem_mask)
dem_min = np.min(dem_mask)
demrec_point = []#存放分界点
for ii in range(dem_min,dem_max+1): #构建分类临界点
    if ii% 50 == 0: #分类区间为50米
        demrec_point.append(ii)
        print(ii)
demrec_point.append(dem_max+1)
demrec_point.insert(0,dem_min-1)        
for rr in range(rows):  #执行重分类
    for cc in range(columns):
        reclass1 = bisect.bisect_left(demrec_point,dem_mask[rr,cc])
        dem_reclass[rr,cc] = reclass1
        
outtifpath1 = 'G:/data/dem/dem_reclass.tif'
driver1 = gdal.GetDriverByName("GTiff")
dataset1 = driver1.Create(outtifpath1, im_width, im_height, 1, gdal.GDT_Float64)
dataset1.SetGeoTransform(geotransform)
dataset1.SetProjection(projection)
dataset1.GetRasterBand(1).WriteArray(dem_reclass)

该代码可以根据自己需求修改。