Python读取栅格数据、进行重分类、阈值分割、分区统计

其实这些Python的栅格操作总体上是以处理数组的思想去处理的，下面就放代码你细细品味

1. 背景介绍

以实例来复现这个过程：使用遥感的比值指数进行冰川识别，对不同高程范围的冰川进行统计分析，流程图如下：

2. 代码实操

①首先以gdal库去读取高程栅格

(GDAL库是处理栅格地理数据主要的Python库，作用非同一般，可以多查阅一下，还有就是处理矢量数据的ogr)

demfile=r’dem.tif’ #栅格数据文件路径
Dem=gdal.Open(demfile) #GDAL库打开栅格文件
proj=dem.GetProjection() #得到栅格数据的投影
geotrans=dem.GetGeoTransform() #得到仿射变换，计算坐标偏移
dem=dem.GetRasterBand(1).ReadAsArray().astype(float) #读取栅格的第一波段，读取后以数组形式存储，类型为浮点型
dem_mask=np.ma.masked_where(dem==32767,dem)#将背景值(NoData)为32767的掩膜掉

②次之对高程数据重分类

demrec=dem #将dem栅格矩阵重新赋值给另外一个变量
demrec[demrec<2000]=0 #高程小于2000的重分类为0
demrec[demrec==32767]=-1 #背景NoData值为32767的重分类为0
#其他范围高程值进行循环赋值
for i in range(14):
    demrec[((2000+i*500)<=demrec)&(demrec<(2000+(i+1)*500))]=i+1

③计算积雪比值指数并进行“阈值分割”提取积雪范围

#GDAL读取2015年遥感影像
LC815=r'image/landsat/lc8150930.tif'
LC815=gdal.Open(LC815)
LC815_b4=LC815.GetRasterBand(4).ReadAsArray().astype(float) #读取第四波段红光波段
LC815_b6=LC815.GetRasterBand(6).ReadAsArray().astype(float) #读取第六波段短波红外波段

red=LC815.GetRasterBand(4).ReadAsArray().astype(float)
green=LC815.GetRasterBand(3).ReadAsArray().astype(float)
blue=LC815.GetRasterBand(2).ReadAsArray().astype(float)
data2015=np.dstack((red,green,blue)) #红绿蓝三波段真彩色显示
plt.subplot(234)
plt.imshow(data2015)

#计算比值指数
LC815_b5 = np.ma.masked_where(LC815_b6== 0, LC815_b6)  
LC815_ratio = LC815_b4/LC815_b6 

#阈值分割
LC815ratio_rec=LC815_ratio #将
# 阈值分割操作，将符合某一范围的值赋值为1，其他赋值为0
LC815ratio_rec[LC815ratio_rec<1.4]=0 
LC815ratio_rec[LC815ratio_rec>=1.4]=1 
plt.subplot(236)
plt.imshow(LC815ratio_rec,cmap='Blues')
plt.title('2015年NDSI阈值分割结果图',fontsize=25)
plt.show()

④分区统计

glacier15_band=LC815ratio_rec.flatten()
glacier_bins = get_bins(glacier15_band) 

#Compute histogram 
hist, dem_bins2, glacier_bins2, bn = scipy.stats.binned_statistic_2d( dem_band, glacier15_band, 
    glacier15_band, 'count', [dem_bins, glacier_bins]) 
#Add bin info to histogram 
hist = np.insert(hist, 0, glacier_bins[:-1], 0) 
row_labels = np.insert(dem_bins[:-1], 0, 0) 
hist = np.insert(hist, 0, row_labels, 1) 

out_fn2015=r'mid/dem_glacier2015.csv'
#Save output 
np.savetxt(out_fn2015, hist, fmt='%1.0f', delimiter=',') 

dataFrame1 = pd.read_csv(r'mid/dem_glacier2015.csv') 
dataFrame1['area'] = dataFrame1['1']*0.0009
dataFrame1['area'].plot(color='red',linewidth=2,label='2015 year')
x=[1,3,5,7,9,11,13,15]
labels=['0-2000','2500-3000','3500-4000','4500-5000','5500-6000','6500-7000','7500-8000','8500-9000']
plt.xticks(x,labels,fontsize=20)
plt.title('2000和2015年不同高程带冰川分布面积',fontsize=25)
plt.legend(fontsize=23)
print(dataFrame1['area'] )