Preface
graduation project has done a streamlined batch processing of Landsat data radiation calibration, atmospheric correction, cropping, NDVI and other index calculations. Now the task is to realize the correction of GF-7 data, some of which need to be modified, just take this opportunity to put Sort this whole thing out. The name of the 6S atmospheric correction python package is Py6s, which is developed and maintained by Dr. RT Wilson* of Southampton. The specific introduction and use can be viewed by yourself, and there are standardized examples and practical application cases on the official website. At present, the processing parameters of various sensors of Landsat data are integrated in the package , which can be called directly, but the correction parameters of domestic satellites such as gain (gain) and bias value (bias), It is only published in China Resources Satellite Application Center, and the website also publishes the spectral response function(* SRF*), you need to sort it out yourself to perform 6s correction processing. As for the principle of more atmospheric correction, please check the paper or other blogs if necessary. Here are a few recommended, such as The Heart of Remote Sensing, ice grapefruit and so on.
text
Only by standing on the shoulders of giants can we move forward better. The predecessors in China have already done the 6s atmospheric correction engineering module of python. The author is a Ph.D. of the State Key Laboratory of Remote Sensing of Wuhan University. View in github repository. The files in this warehouse can be directly corrected by running the corresponding python file, but the current script can only do single-frame processing. However, on the basis of this script, Jiangbei stores multiple images to be processed in a txt file, and uses the txt file to perform BAT Script batch processing. Although the 6s batch atmospheric correction is also implemented in disguise, I personally feel that it is not very pythonic, so I still want to go through the pure python code of batch processing with my own ideas.
Secondary functions
Give a few secondary functions first, the functions that serve the main function.
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
|
def MeanDEM ( pointUL , pointDR , GMTEDdem ):
'''
Calculates the average elevation of the area where the image is located.
'''
try:
DEMIDataSet = gdal.Open ( GMTEDdem )
except Exception as e:
pass
DEMBand = DEMIDataSet.GetRasterBand (1)
geotransform = DEMIDataSet.GetGeoTransform ()
# DEM resolution
pixelWidth = geotransform [1]
pixelHight = geotransform [5]
# DEM starting point: upper left corner, X: longitude, Y: latitude
originX = geotransform [0]
originY = geotransform [3]
# The position of the upper left corner of the study area in the DEM matrix
yoffset1 = int(( originY - pointUL [' lat ']) / pixelWidth )
xoffset1 = int(( pointUL [' lon '] - originX ) / (- pixelHight ))
# The position of the lower right corner of the study area in the DEM matrix
yoffset2 = int(( originY - pointDR [' lat ']) / pixelWidth )
xoffset2 = int(( pointDR [' lon '] - originX ) / (- pixelHight ))
# The number of rows and columns of the study area matrix
xx = xoffset2 - xoffset1
yy = yoffset2 - yoffset1
# Read the data in the study area and calculate the elevation
DEMRasterData = DEMBand.ReadAsArray (xoffset1, yoffset1, xx, yy )
MeanAltitude = np.mean ( DEMRasterData )
return MeanAltitude
def untar ( fname , dirs ):
''' Decompress tar file function'''
t = tarfile.open ( fname )
t.extractall (path= dirs )
|
These two functions are basically unchanged, only the untar is changed to pass the parameter system.
main function
The following describes the changes of several main functions.
Radiometric Calibration
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
def RadiometricCalibration ( BandId,SatelliteID,SensorID,Year,config ):
''' Get the radiometric correction gain and offset value;
BANDID : nth band
SatelliteID : satellite platform, stored in json
SensorID : sensor type, stored in json
of the scaling factor, stored in json
config : scaling parameter storage json file path'''
if SensorID [0:3] == " bwd ":
Gain_ =config["Parameter"][ SatelliteID ][ SensorID ][Year]["gain"]
Bias_ =config["Parameter"][ SatelliteID ][ SensorID ][Year]["offset"]
else:
Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"][BandId-1]
Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"][BandId-1]
return Gain_, Bias_
|
It is roughly the same as the original script, but the discriminant function has been changed in order to be used for different sensors with high score data. In the future, the discrimination and json parameters of the forward-looking sensor (FWD) will be added. In fact, it is enough to add elif to form multiple discriminations.
Atmospheric Correction
Next comes the atmospheric correction:
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
def AtmosphericCorrection(BandId, metedata, config, SatelliteID, SensorID, GMTEDdem):
'''
Atmospheric Correction Function
'''
dom = xml.dom.minidom.parse ( metedata )
# 6S model
s = SixS ()
# Sensor type customization
s.geometry = Geometry.User ()
s.geometry.solar_z = float(dom.getElementsByTagName('SunAltitude')[0].firstChild.data)
s.geometry.solar_a = float( dom.getElementsByTagName (' SunAzimuth ')[0] .firstChild.data )
# s.geometry.view_z = float(dom.getElementsByTagName('SatelliteZenith')[0].firstChild.data)
# s.geometry.view_a = float(dom.getElementsByTagName('SatelliteAzimuth')[0].firstChild.data)
s.geometry.view_z = 0
s.geometry.view_a = 0
# date
DateTimeparm = dom.getElementsByTagName (' CenterTime ')[0] .firstChild.data
DateTime = DateTimeparm.split (' ')
Date = DateTime [0].split('-')
s.geometry.month = int(Date[1])
s.geometry.day = int(Date[2])
# print( s.geometry )
# Center latitude and longitude\
LeftTop = dom.getElementsByTagName (' LeftTopPoint ')[0]
TopLeftLat = float(LeftTop.getElementsByTagName('Latitude')[0].firstChild.data)
TopLeftLon = float(LeftTop.getElementsByTagName('Longtitude')[0].firstChild.data)
RightTop = dom.getElementsByTagName (' RightTopPoint ')[0]
TopRightLat = float(RightTop.getElementsByTagName('Latitude')[0].firstChild.data)
TopRightLon = float(RightTop.getElementsByTagName('Longtitude')[0].firstChild.data)
BottomRight = dom.getElementsByTagName (' RightBottomPoint ')[0]
BottomRightLat = float(BottomRight.getElementsByTagName('Latitude')[0].firstChild.data)
BottomRightLon = float(BottomRight.getElementsByTagName('Longtitude')[0].firstChild.data)
Bottomleft = dom.getElementsByTagName (' LeftBottomPoint ')[0]
BottomLeftLat = float(Bottomleft.getElementsByTagName('Latitude')[0].firstChild.data)
BottomLeftLon = float(Bottomleft.getElementsByTagName('Longtitude')[0].firstChild.data)
ImageCenterLat = ( TopLeftLat + TopRightLat + BottomRightLat + BottomLeftLat ) / 4
# Atmospheric mode type
if ImageCenterLat > -15 and ImageCenterLat < 15:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
if ImageCenterLat > 15 and ImageCenterLat < 45:
if s.geometry.month > 4 and s.geometry.month < 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)
if ImageCenterLat > 45 and ImageCenterLat < 60:
if s.geometry.month > 4 and s.geometry.month < 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)
# aerosol type continent
s.aero_profile = AtmosProfile.PredefinedType ( AeroProfile.Continental )
# Underlying surface type
s.ground_reflectance = GroundReflectance.HomogeneousLambertian (0.36)
# 550nm aerosol optical thickness, corresponding to a visibility of 40km
s.aot550 = 0.14497
# Find the DEM height by studying the range of the removed area.
pointUL = dict ()
pointDR = dict ()
pointUL [" lat "] = max( TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat )
pointUL [" lon "] = min( TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon )
pointDR [" lat "] = min( TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat )
pointDR [" lon "] = max( TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon )
meanDEM = ( MeanDEM ( pointUL , pointDR , GMTEDdem )) * 0.001
# Study area altitude, satellite sensor orbit altitude
s.altitudes = Altitudes()
s.altitudes.set_target_custom_altitude ( meanDEM )
s.altitudes.set_sensor_satellite_level ()
# Correction band (according to band name) panchromatic 440-905
if BandId == 1:
SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["1"]
s.wavelength = Wavelength(0.450,0.520,SRFband)
elif BandId == 2:
SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["2"]
s.wavelength = Wavelength(0.520,0.590,SRFband)
elif BandId == 3:
SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["3"]
s.wavelength = Wavelength(0.630,0.690,SRFband)
elif BandId == 4:
SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["4"]
s.wavelength = Wavelength(0.770,0.890,SRFband)
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance (-0.1)
# run 6s atmospheric model
s.run ()
xa = s.outputs.coef_xa
xb = s.outputs.coef_xb
xc = s.outputs.coef_xc
# x = s.outputs.values
return ( xa , xb , xc)
|
It is worth mentioning here that the solar zenith angle ( SolarZenith ) and solar azimuth ( SolarAzimuth ) parameters of xml in GF1\GF2 In the newer data of high score, the xml tags have been changed to SunAltitude and SunAzimuthThat is, the zenith angle has been changed to the height angle, so that in the original script
1
|
s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data)
|
For the writing method, 90- must be removed. The specific principle can be viewed here.
Another point is that in the data above GF6, the longitude and latitude labels of the image vertex coordinates of the header file xml are assigned to their respective position labels, comparing the original writing method and the new writing method:
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
|
#Original spelling
# center latitude and longitude
TopLeftLat = float(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data)
TopLeftLon = float(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data)
TopRightLat = float(dom.getElementsByTagName('TopRightLatitude')[0].firstChild.data)
TopRightLon = float(dom.getElementsByTagName('TopRightLongitude')[0].firstChild.data)
BottomRightLat = float(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data)
BottomRightLon = float(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data)
BottomLeftLat = float(dom.getElementsByTagName('BottomLeftLatitude')[0].firstChild.data)
BottomLeftLon = float(dom.getElementsByTagName('BottomLeftLongitude')[0].firstChild.data)
#new writing
# center latitude and longitude
LeftTop = dom.getElementsByTagName (' LeftTopPoint ')[0]
TopLeftLat = float(LeftTop.getElementsByTagName('Latitude')[0].firstChild.data)
TopLeftLon = float(LeftTop.getElementsByTagName('Longtitude')[0].firstChild.data)
RightTop = dom.getElementsByTagName (' RightTopPoint ')[0]
TopRightLat = float(RightTop.getElementsByTagName('Latitude')[0].firstChild.data)
TopRightLon = float(RightTop.getElementsByTagName('Longtitude')[0].firstChild.data)
BottomRight = dom.getElementsByTagName (' RightBottomPoint ')[0]
BottomRightLat = float(BottomRight.getElementsByTagName('Latitude')[0].firstChild.data)
BottomRightLon = float(BottomRight.getElementsByTagName('Longtitude')[0].firstChild.data)
Bottomleft = dom.getElementsByTagName (' LeftBottomPoint ')[0]
BottomLeftLat = float(Bottomleft.getElementsByTagName('Latitude')[0].firstChild.data)
BottomLeftLon = float(Bottomleft.getElementsByTagName('Longtitude')[0].firstChild.data)
|
In fact, it is equivalent to turning a corner, but it must be written to remind other developers to pay attention to viewing the xml of the image, so as to extract the correct content in time.
Block
With the above two most important functions, the block processing can be performed. This is also directly used in the original script, but there are several more discriminant functions to automatically discriminate the remote sensing platform and sensor, so as to correctly call the parameters:
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
def Block(IDataSet, SensorID, Year, SatelliteID, atcfiles, config, metedata, outFileName, GMTEDdem ):
''' block handler
With block operation efficiency, a GF-7-MUX data can be corrected within 2 minutes
'''
# global cols,rows,atcfiles
cols = IDataSet.RasterXSize
rows = IDataSet.RasterYSize
#set output band
Driver = IDataSet.GetDriver ()
geoTransform1 = IDataSet.GetGeoTransform ()
ListgeoTransform1 = list(geoTransform1)
ListgeoTransform1[5] = -ListgeoTransform1[5]
newgeoTransform1 = tuple(ListgeoTransform1)
proj1 = IDataSet.GetProjection ()
OutRCname = os.path.join (atcfiles, outFileName + " .tif ")
outDataset = Driver.Create (OutRCname,cols,rows,4,gdal.GDT_Int32)
outDataset.SetGeoTransform (newgeoTransform1)
outDataset.SetProjection (proj1)
#Read 4 bands separately, multi-band file correction
if SensorID == "MUX":
for m in range(1,5):
ReadBand = IDataSet.GetRasterBand (m)
outband = outDataset.GetRasterBand (m)
outband.SetNoDataValue (-9999)
#Get the gain and offset bias of the corresponding band
Gain,Bias = RadiometricCalibration ( m,SatelliteID,SensorID,Year,config )
#Get atmospheric correction coefficient
AtcCofa , AtcCofb , AtcCofc = AtmosphericCorrection ( m, metedata , config, SatelliteID , SensorID , GMTEDdem )
nBlockSize = 2048
i = 0
j = 0
b = cols*rows
print("The %d band correction: "%m)
try:
while i <rows:
while j <cols:
#Save chunk size
nXBK = nBlockSize
nYBK = nBlockSize
#The last area that is not enough to be divided into blocks, how many reads are there
if i+nBlockSize >rows:
nYBK = rows - i
if j+nBlockSize >cols:
nXBK =cols - j
#Read the image in chunks
Image = ReadBand.ReadAsArray ( j,i,nXBK,nYBK )
outImage = np.where (Image>0,Image*Gain + Bias,-9999)
y = np.where ( outImage !=-9999,AtcCofa * outImage - AtcCofb,-9999)
atcImage = np.where (y!=-9999,(y / (1 + y * AtcCofc ))*10000,-9999)
outband.WriteArray ( atcImage,j,i )
j= j+nXBK
time.sleep (1)
j=0
i = i+nYBK
except KeyboardInterrupt :
raise
#Single band file correction
elif SensorID == " nad " or SensorID == " bwd " or SensorID == " fwd ":
ReadBand = IDataSet.GetRasterBand (1)
outband = outDataset.GetRasterBand (1)
outband.SetNoDataValue (-9999)
Gain, Bias = RadiometricCalibration(ReadBand, SatelliteID, SensorID, Year, config)
nBlockSize = 2048
i = 0
j = 0
b = cols*rows
print("Single band correction")
while i <rows:
while j <cols:
#Save chunk size
nXBK = nBlockSize
nYBK = nBlockSize
#The last area that is not enough to be divided into blocks, how many reads are there
if i+nBlockSize >rows:
nYBK = rows - i
if j+nBlockSize >cols:
nXBK =cols - j
#Read the image in chunks
Image = ReadBand.ReadAsArray ( j,i,nXBK,nYBK )
outImage = np.where (Image>0,Image*Gain + Bias,-9999)
outband.WriteArray ( outImage,j,i )
j= j+nXBK
time.sleep (1)
j=0
i = i+nYBK
print("Radiometric calibration completed, no atmospheric correction required")
|
The writing here is still not rigorous enough. Although the effect of atmospheric correction in the panchromatic band is not so significant, since the official also gave the spectral response function and gain bias, it should also be judged to perform radiation and atmospheric correction, but it is stolen here. I’m lazy, I haven’t finished writing all the sensor judgments.
Main function
Finally, in order to be more concise and make the parameters input more flexibly, here is a main function written by myself :
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
74
75
76
77
78
79
80
81
82
83
84
|
def read_output_file(input_tarfile,config_file,OutputFilePath,outFileName,GMTEDdem):
'''
main function
input_tarfile : tar compressed file path;
config_file : path of radiation correction parameter file;
OutputFilePath : the path to save the output file;
outFileName : output file save name;
GMTEDdem : DEM file required for correction.
'''
config = json.load (open( config_file ))
filename = os.path.basename ( input_tarfile )
fileType = filename[0:3]
filename_split = filename.split ("_")
base_dir = os.path.dirname ( input_tarfile )
if fileType == 'GF7':
#GFType = filename[4:7]
#GFType = filename_split [5]
#GF701_011628_E107.2_N22.9_20211202114945_BWD_01_SC0_0001_2112039662
SensorID = filename_split [5]
Year = filename_split [4][:4]
SatelliteID = filename_split [0]
untarfilename = filename[:-7]
outpath = os.path.join ( base_dir, untarfilename )
atcfiles = os.path.join ( OutputFilePath,outFileName )
untar ( input_tarfile , outpath )
tiffFile = glob.glob ( os.path.join ( outpath ,"*. tif "))[0]
metedata = glob.glob ( os.path.join ( outpath ,"*.xml"))[0]
rpb_file = glob.glob ( os.path.join ( outpath ,"*_rpc.txt"))[0]
try:
os.mkdir ( atcfiles )
except Exception as e:
pass
#Copy the header file and geometric correction file to the atmospheric correction result file
metedata_basename = os.path.basename ( metedata )
copy_metedata = os.path.join ( atcfiles,metedata_basename )
shutil.copy ( metedata, copy_metedata )
rpb_basename = os.path.basename ( rpb_file )
copy_rpb_file = os.path.join ( atcfiles,rpb_basename )
shutil.copy ( rpb_file, copy_rpb_file )
try:
IDataSet = gdal.Open ( tiffFile )
except Exception as e:
print("Failed to open file %S" % tiffFile )
elif fileType == 'zy3':
SensorID = filename_split [1]
Year = filename_split [4][:4]
SatelliteID = filename_split [0]
untarfilename = filename[:-7]
outpath = os.path.join ( base_dir, untarfilename )
atcfiles = os.path.join ( OutputFilePath,outFileName )
untar ( input_tarfile , outpath )
tiffFile = glob.glob ( os.path.join ( outpath ,"*. tif "))[0]
metedata = glob.glob ( os.path.join ( outpath ,"*.xml"))[0]
rpb_file = glob.glob ( os.path.join ( outpath ,"*_rpc.txt"))[0]
try:
os.mkdir ( atcfiles )
except Exception as e:
pass
#Copy the header file and geometric correction file to the atmospheric correction result file
metedata_basename = os.path.basename ( metedata )
copy_metedata = os.path.join ( atcfiles,metedata_basename )
shutil.copy ( metedata, copy_metedata )
rpb_basename = os.path.basename ( rpb_file )
copy_rpb_file = os.path.join ( atcfiles,rpb_basename )
shutil.copy ( rpb_file, copy_rpb_file )
try:
IDataSet = gdal.Open ( tiffFile )
except Exception as e:
print("Failed to open file %S" % tiffFile )
Block(IDataSet,SensorID,Year,SatelliteID,atcfiles,config,metedata,outFileName,GMTEDdem)
|
The idea of writing is also referring to the if __main__ part of the script. After encapsulating it into a function, directly giving the 5 parameters of the function, you can easily perform customized batch processing or other process applications in python.
1
2
3
4
5
6
|
Dem = R".GMTED2km.tif"
outFileName = ".123"
OutputFilePath =r" .outtar "
config_file = r" .RadiometricCorrectionParameter.json "
input_tarfile = r".GFXXXXXX.tar.gz "
read_output_file ( input_tarfile , config_file , OutputFilePath , outFileName , Dem)
|
supplementary data
Atmospheric correction is a complex physical mechanism, which requires many parameters. Most of the parameters are given in the 6s package, and there are still some parameters that need to be input by yourself, such as DEM and spectral response function. The use of auxiliary data is introduced here.
DEM
The elevation data is very simple. If you need high-precision, you can directly upload high-resolution dem. The original script uses the USGS GMTED elevation data with a resolution of 2km. If you have the conditions, you can directly use SRTM or Japan ALOS2. Yes, but make sure to cover the study area.
Radiometric calibration and spectral response json file
The radiometric calibration parameters and the parameters required by the spectral response function are stored in the RadiometricCorrectionParameter.json file in json format , and its structure is as follows:
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
|
{
"Parameter": {
"GF1": {
"WFV1": {
"2013": {
"gain": [
5.851,
7.153,
8.368,
7.474
],
"offset": [
0.0039,
0.0047,
0.003,
0.0274
]
},
"SRF": {
"1": [
0.46233000000000113,
...
0.10741000000000049
],
"2": [
0.7663399999999977,
....
]
}
},
|
Note here that the spectral response function of the Py6s package:
1
|
s.wavelength = Wavelength(0.770,0.890,SRFband)
|
Among them, the correction band range is 0.770~0.890nm, and its data resolution is 2.5nm, and the official parameters are often 1nm resolution, so it needs to be interpolated to 0.5, and then the interval is 2.5, which is also what needs to be done later. After downloading the required SRF and radiometric calibration parameters from the [Satellite Center Official Website] (http://www.cresda.com/CN/Downloads/dbcs/index.shtml), arrange them in json data format. It can be easily called, and data can be added and customized according to different sensors and remote sensing platforms.
Afterword
At present, the open source atmospheric correction algorithm, only 6s is relatively reliable, and the python package makes the implementation of atmospheric correction simpler, so it is also very important to Southampton Professor Williamson and thanks to Dr. Zhao for open source. In the future, consider adding the discrimination of GF and ZY series sensors, sort out the json and then commit to github . On the other hand, the rapid atmospheric correction and deep learning atmospheric correction of domestic satellites also need to be studied. These are the current thoughts. If you have any interesting content that you would like to discuss with me, please send an email to communicate with me.
