Python calculates sun and sensor azimuth and zenith angles for each pixel of Landsat8

After experiencing the huge project of migrating from Hexo to hugo, I found that I have been working on the whole framework configuration without putting the article into practice, so now I write the content that I have been hoarding first. This article introduces the generation of image-by-image metadata based on python-l8-angles for the Landsat series data solar altitude angle and other poses, and finally generates raster images of remote sensor poses based on this package.

Introduction

The Sentinel image processing software SNAP Desktop of ESA has a processing module called Biophysical Processor, which is used to invert various vegetation physiological parameters such as leaf area index (LAI), chlorophyll (Cab), light and effective absorbed radiation ratio (FAPAR), vegetation cover (FCOVER), and canopy water content (CW). The module has been used to generate biophysical products for the VEGETATION, MERIS, SPOT and LANDSAT sensors. It consists mainly of generating a comprehensive database of vegetation characteristics and associated parameters, SENTINEL2 or LANDSAT8 top-of-canopy (TOC) reflectance. The neural network is then trained to estimate canopy features based on TOC reflectance and the corresponding angles of the settings defining the observation configuration. For Sentinel-2, two different neural network architectures were implemented, NNET 20m and NNET 10m. In addition to the auxiliary bands cos(viewing_zenith), cos(sun_zenith), and cos(relative_azimuth_angle), both use Sentinel -2 itself comes with different combinations of 12 bands as feature inputs. More details can be found in the official documentation 。

Sentry data itself comes with sensor attitude band data in the documentation while Landsat-8 does not, so processing Landsat data with this module generates “Message: Product ‘LC08_ L1TP_ 124044 20200506 20200509 01 T1 resampled’ does not contain a raster with name ‘view zenith _mean’”. snap_error1

we searched the official USGS documentation for this requirement and finally found the processing package for generating attitude raster data for Landsat series data, which is described below.

L8-ANGLES-PYTHON packages

This package implements a simple Python interface to the USGS Landsat 8 tool for calculating solar and sensor azimuth and zenith angles per pixel from the angle coefficient file “ANG.TXT”. If you are using anaconda (which is what I use anyway), installing the package is as simple as a single command.

1
2


conda install py-l8angles
conda config --add channels DHI-GRAS #or

On the Github page, the official researchers provide a simple edible example.

1
2
3
4


import l8angles

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

Programmers familiar with python will be able to customize the content according to their needs here

Generating pose raster data in python based on L8-angles package

We need to import the relevant packages first

1
2
3
4
5
6


import l8angles
import os
import rasterio
from rasterio import transform
from rasterio import plot
import numpy as np

Because it is learned python early learned, so used two for loop, got landsat file path, for the subsequent batch processing to prepare for.

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'

Next, import the file, extract the reference information from the raster image, calculate the attitude angle raster data, and place the generated attitude raster data in.

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)

Final output attitude raster data file.

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()       

The above are all the steps to achieve the generated image is shown in the figure below, which is actually the cover image. Sun_az

Full Code

Finally, the complete code is attached, so that the big boys can customize the batch method according to their needs:.

 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) #Iterates through the names of the folders in the entire folder and returns a list
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()       

Appendix

The reference blog literature reads as follows：