Featured image of post Python计算Landsat8每个像素的太阳和传感器方位角和天顶角

Python计算Landsat8每个像素的太阳和传感器方位角和天顶角

python生成Landsat姿态角数据,介绍L8-angles python包.

经历了从Hexo迁移到hugo的庞大工程以后,发现一直以来都在整框架配置,没有把文章落到实处,因此现在把一直囤积的内容先写了。本文介绍基于python-l8-angles生成Landsat系列数据太阳高度角等姿态的逐像元数据,最终基于该包生成遥感器姿态的栅格图像。

引言

欧空局(OSA)的Sentinel 哨兵系列影像处理软件SNAP Desktop有个处理模块叫做Biophysical Processor,用于反演哨兵数据和Landsat-8 OLI数据的各种植被生理参数如叶面积指数(LAI)、叶绿素(Cab)、光和有效吸收辐射比(FAPAR)、植被覆盖度(FCOVER)、冠层含水量(CW)。该模块已用为VEGETATION、MERIS、SPOT 和 LANDSAT 传感器生成生物物理产品。它主要包括生成植被特征和相关的参数,SENTINEL2 或LANDSAT8 冠层顶部 (TOC) 反射率的综合数据库。 然后训练神经网络以根据 TOC 反射率以及定义观察配置的设置相应角度来估计冠层特征。对于 Sentinel-2,实现了2种不同的神经网络架构,即 NNET 20m 和 NNET 10m。除了辅助波段 cos(viewing_zenith)、cos(sun_zenith)、cos(relative_azimuth_angle) 之外,两者还使用Sentinel-2本身自带的12个波段的不同组合作为特征输入。更多详细内容可以去官方文档 参阅。

哨兵数据本身文档里就自带了传感器姿态波段数据而Landsat-8却没有,因此在用该模块处理Landsat数据时会产生 “Message: Product ‘LC08_ L1TP_ 124044 20200506 20200509 01 T1 resampled’ does not contain a raster with name ‘view zenith _mean’”。 snap_error1

本文根据这一需求查找了USGS官方文档,终于发现了对Landsat系列数据生成姿态栅格数据的处理包,下面对该包进行介绍。

L8-ANGLES-PYTHON包

这个包实现了一个简单的 Python 接口到 USGS Landsat 8 工具,用于计算每个像素的太阳和传感器方位角和天顶角,来自角度系数文件“ANG.TXT”。如果是用anaconda(反正我是用这),安装该包很简单,只需要一行命令:

1
2
conda install py-l8angles
conda config --add channels DHI-GRAS #或者

Github page上,官方研究人员提供了一个简单的食用例子:

1
2
3
4
import l8angles

data = l8angles.calculate_angles('./test_ANG.txt', angle_type='SOLAR',
                                 subsample=2, bands=[3,6,7])

熟悉python的程序员(大佬),到这里就能够根据自己需求往下写了,然而小白的(和我一样的菜鸡),跟我接着往下读一读吧。

基于L8-angles包的python 生成姿态栅格数据

不多说,还是基于python批处理思想,先导包:

1
2
3
4
5
6
import l8angles
import os
import rasterio
from rasterio import transform
from rasterio import plot
import numpy as np

因为是学python早期学的了,因此用了两个for循环,拿到了landsat文件路径,为后续批处理做准备:

1
2
3
4
5
6
7
8
9
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:]=='B3':
            band3path = product_path + '/' + str(productIO) + '.TIF'

接下来导入文件,提取栅格影像的参考信息,计算姿态角栅格数据,并置入生成的姿态栅格数据:

1
2
3
4
5
band3 = rasterio.open(band3path)
geo = band3.read(1).astype('float64')
data = l8angles.calculate_angles('F:/6s_test/L8angle/LC81240422014350LGN01/LC08_L1TP_124042_20141216_20170416_01_T1_ANG.txt', angle_type='SOLAR',subsample=1, bands=[3])
sa=data['sun_az']
m = np.array(sa)

最后输出姿态栅格数据文件:

1
2
3
4
5
a_output_path = "F:/6s_test/output"
output_a = a_output_path + '/' +  'sun_as.tif'
outputfile = rasterio.open(output_a,'w',driver='Gtiff',crs = band3.crs,transform=band3.transform,dtype='float64',width=band3.width,height = band3.height,count=1)
outputfile.write(m)   
outputfile.close()       

以上是所有实现的步骤生成的影像见下图,其实就是封面图。Sun_az

完整代码

最后附上完整代码,大佬们可根据自己的需求,定制批处理方法:

 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
import l8angles
import os
import rasterio
from rasterio import transform
from rasterio import plot
import numpy as np
productIOs = []
filepath = "F:/6s_test/L8angle"
pathlist = os.listdir(filepath) #遍历整个文件夹下文件夹的名字  并返回一个列表
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:]=='B3':
            band3path = product_path + '/' + str(productIO) + '.TIF'
    band3 = rasterio.open(band3path)
    geo = band3.read(1).astype('float64')
    data = l8angles.calculate_angles('F:/6s_test/L8angle/LC81240422014350LGN01/LC08_L1TP_124042_20141216_20170416_01_T1_ANG.txt', angle_type='SOLAR',subsample=1, bands=[3])
    sa=data['sun_az']
    m = np.array(sa)




    a_output_path = "F:/6s_test/output"
    output_a = a_output_path + '/' +  'sun_as.tif'
    outputfile = rasterio.open(output_a,'w',driver='Gtiff',crs = band3.crs,transform=band3.transform,dtype='float64',width=band3.width,height = band3.height,count=1)
    outputfile.write(m)   
    outputfile.close()       

附录

参考博客文献内容如下:

  1. Landsat 8 Angles Creation Tools README - USGS
  2. 一个简单的Python接口到 USGS Landsat8 工具
  3. 基于 Python 的遥感图像NDVI批处理
  4. S2ToolBox Level 2 products:LAI, FAPAR, FCOVER
Author by Jerrychoices
Built with Hugo
主题 StackJimmy 设计

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