当前位置:   article > 正文

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

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的掩膜掉
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

在这里插入图片描述

②次之对高程数据重分类

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
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

在这里插入图片描述

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

#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()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26

在这里插入图片描述

④分区统计

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'] )
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24

在这里插入图片描述
在这里插入图片描述

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/weixin_40725706/article/detail/81088
推荐阅读
相关标签
  

闽ICP备14008679号