from osgeo import gdal
from osgeo import gdal_array
import math
import os
import rasterio
from matplotlib import pyplot as plt
定义影像所在路径,以及相关参数
#路径
tpath='/gdata'
# Use the same example image:
date = '2017-06-16'
url = '/gdata/'
redband = 'LC08_L1TP_117031_20150305_20170412_01_T1_B{}.TIF'.format(4)
nirband = 'LC08_L1TP_117031_20150305_20170412_01_T1_B{}.TIF'.format(5)
查看红外波段信息:
with rasterio.open(url+redband) as src:
profile = src.profile
oviews = src.overviews(1) # list of overviews from biggest to smallest
oview = oviews[1] # Use second-highest resolution overview
print('Decimation factor= {}'.format(oview))
red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))
plt.imshow(red)
plt.colorbar()
plt.title('{}\nRed {}'.format(redband, red.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')
查看近红外波段信息:
with rasterio.open(url+nirband) as src:
oviews = src.overviews(1) # list of overviews from biggest to smallest
oview = oviews[1] # Use second-highest resolution overview
nir = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))
plt.imshow(nir)
plt.colorbar()
plt.title('{}\nNIR {}'.format(nirband, nir.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')
定义计算 NDVI的函数,并进行计算:
def calc_ndvi(nir,red):
'''Calculate NDVI from integer arrays'''
nir = nir.astype('f4')
red = red.astype('f4')
ndvi = (nir - red) / (nir + red)
return ndvi
ndvi = calc_ndvi(nir,red)
查看结果:
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar()
plt.title('NDVI {}'.format(date))
plt.xlabel('Column #')
plt.ylabel('Row #')
Comment list ( 0 )