Featured image of post Python批量计算Landsat NDVI并输出tiff
python遥感处理

Python批量计算Landsat NDVI并输出tiff

利用rasterio包单独完成Landsat数据的NDVI、FPAR、Cal等植物生物物理参数批量计算.

Rasterio包介绍与安装

Rasterio包能够读取和写入GeoTIFF和其他格式，并返回 Numpy N 维数组或 GeoJSON 地理json，是一个基于GDAL开发的强大的栅格数据处理python包。

Anaconda安装该包仅使用一行命令即可：

1

conda install rasterio

如果上面命令安装报错的话，可以试一下以下任意一行conda官方频道安装：

1
2
3
4
5
6

conda install -c conda-forge rasterio
conda install -c conda-forge/label/cf202003 rasterio
conda install -c conda-forge/label/dev rasterio
conda install -c conda-forge/label/broken rasterio
conda install -c conda-forge/label/rasterio_dev rasterio
conda install -c conda-forge/label/cf201901 rasterio

其他使用可以查看官方DOCS

计算NDVI代码

首先导入所需要的包

1
2
3
4
5
6

import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np
import os
from rasterio import transform

将Landsat数据预处理完成后，放在同一个文件夹下，由于本人是利用Python进行的批量辐射定标和6s大气校正，因此得到的输出只有geotiff单幅影像，将其全部放在同一文件夹后，得到的待计算数据如下所视。

文件图

解析每个文件夹下，所需要的波段文件，其中计算植物物理参数的波段总共需要B4，B5，B3三个波段。

 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
27

filepath = "D:/Data/Landsat_8_retry_cal_6s"
productIOs = []

pathlist = os.listdir(filepath) #遍历整个文件夹下文件夹的名字  并返回一个列表
#print(pathlist)
for i in pathlist:
    product_path = filepath + '/' + str(i)  
    product_list = os.listdir(product_path)
    for product in product_list:
        productIOs.append(product.split(".")[0])
    for productIO in productIOs:
        if productIO[-2:]=='B4':
            band4path = product_path + '/' + str(productIO) + '.TIF'
            
    for productIO in productIOs:
        if productIO[-2:]=='B5':
            band5path = product_path + '/' + str(productIO) + '.TIF'
    
    for productIO in productIOs:
        if productIO[-2:]=='B3':
            band3path = product_path + '/' + str(productIO) + '.TIF'   
    band4 = rasterio.open(band4path)
    band5 = rasterio.open(band5path)
    band3 = rasterio.open(band3path)
    red = band4.read(1).astype('float64')
    nir = band5.read(1).astype('float64')
    green = band3.read(1).astype('float64')

以上将三个波段读入numpy array数组以后，直接利用python 数学表达式计算各个植被指数即可，代码如下：

1
2
3
4
5
6

#植被指数计算
ndvi =( nir - red )/(nir+red)
fapar = 107.5*ndvi - 8.0
lai = 5*np.log(ndvi)+ 4.22
ci = nir/green - 1
ccc = ci*0.89+7.47

最后是批量输出。将5种计算结果分别输出到指定文件夹，并保存为geotiff格式。实际上几种植被指数的输出写法都一致，只是重复了几遍，有时间会改成函数形式，直接调用即可。

 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
27

#ndvi
ndvi_output_path = "D:/Data/NDVI_files"
output_ndvi = ndvi_output_path + '/' + str(i) + '_ndvi' +'.TIF'
ndviImage = rasterio.open(output_ndvi,'w',driver='Gtiff',width=band4.width,height = band4.height,\
    count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
ndviImage.write(ndvi,1)
ndviImage.close()
#FAPAR
fapar_output_path = "D:/Data/FAPAR_files"
output_fapar = fapar_output_path + '/' + str(i) + '_fapar' +'.TIF'
faparImage = rasterio.open(output_fapar,'w',driver='Gtiff',width=band4.width,height = band4.height,\
    count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
faparImage.write(fapar,1)
faparImage.close()
#LAI
LAI_output_path = "D:/Data/LAI_files"
output_LAI = LAI_output_path + '/' + str(i) + '_LAI' +'.TIF'
LAIImage = rasterio.open(output_LAI,'w',driver='Gtiff',width=band4.width,height = band4.height,\
    count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
LAIImage.write(lai,1)
LAIImage.close()
#CCC
CCC_output_path = "D:/Data/CCC_files"
output_CCC = CCC_output_path + '/' + str(i)  +'.TIF'
CCCImage = rasterio.open(output_CCC,'w',driver='Gtiff',width=band4.width,height = band4.height,count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
CCCImage.write(ccc,1)
CCCImage.close()

完整代码

这里给出完整代码工程，很久之前写的了，有些地方完全可以封装为函数，但是当时没有这种意识，当作学习记录吧。

 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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

#LC8计算ndvi正式版

import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np
import os
from rasterio import transform




filepath = "D:/Data/Landsat_8_retry_cal_6s"
productIOs = []

pathlist = os.listdir(filepath) #遍历整个文件夹下文件夹的名字  并返回一个列表
#print(pathlist)
for i in pathlist:
    product_path = filepath + '/' + str(i)  
    product_list = os.listdir(product_path)
    for product in product_list:
        productIOs.append(product.split(".")[0])
    for productIO in productIOs:
        if productIO[-2:]=='B4':
            band4path = product_path + '/' + str(productIO) + '.TIF'
            
    for productIO in productIOs:
        if productIO[-2:]=='B5':
            band5path = product_path + '/' + str(productIO) + '.TIF'
    
    for productIO in productIOs:
        if productIO[-2:]=='B3':
            band3path = product_path + '/' + str(productIO) + '.TIF'   
    band4 = rasterio.open(band4path)
    band5 = rasterio.open(band5path)
    band3 = rasterio.open(band3path)
    red = band4.read(1).astype('float64')
    nir = band5.read(1).astype('float64')
    green = band3.read(1).astype('float64')
    #植被生理参数
    ndvi =( nir - red )/(nir+red)
    fapar = 107.5*ndvi - 8.0
    lai = 5*np.log(ndvi)+ 4.22
    ci = nir/green - 1
    ccc = ci*0.89+7.47
    #ndvi
    ndvi_output_path = "D:/Data/NDVI_files"
    output_ndvi = ndvi_output_path + '/' + str(i) + '_ndvi' +'.TIF'
    ndviImage = rasterio.open(output_ndvi,'w',driver='Gtiff',width=band4.width,height = band4.height,\
        count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
    ndviImage.write(ndvi,1)
    ndviImage.close()
    #FAPAR
    fapar_output_path = "D:/Data/FAPAR_files"
    output_fapar = fapar_output_path + '/' + str(i) + '_fapar' +'.TIF'
    faparImage = rasterio.open(output_fapar,'w',driver='Gtiff',width=band4.width,height = band4.height,\
        count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
    faparImage.write(fapar,1)
    faparImage.close()
    #LAI
    LAI_output_path = "D:/Data/LAI_files"
    output_LAI = LAI_output_path + '/' + str(i) + '_LAI' +'.TIF'
    LAIImage = rasterio.open(output_LAI,'w',driver='Gtiff',width=band4.width,height = band4.height,\
        count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
    LAIImage.write(lai,1)
    LAIImage.close()
    #CCC
    CCC_output_path = "D:/Data/CCC_files"
    output_CCC = CCC_output_path + '/' + str(i)  +'.TIF'
    CCCImage = rasterio.open(output_CCC,'w',driver='Gtiff',width=band4.width,height = band4.height,count=1,crs = band4.crs,transform=band4.transform,dtype='float64')
    CCCImage.write(ccc,1)
    CCCImage.close()

附录

参考文章：

  1. Rasterio: access to geospatial raster data
  2. python 栅格处理利器之Rasterio——郭晓隐
  3. S2ToolBox Level 2 products:LAI, FAPAR, FCOVER
© 2022 - 2023 jerrychoices 的博客
Author by Jerrychoices
Built with Hugo
主题 Stack 由 Jimmy 设计

本站访客数人次 总访问量次 本文阅读量次