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
|
'''
@Time ： 2023/9/7 16:08
@Auth ： Guo Jiaxiang
@Blog : https://www.guojxblog.cn
@GitHub : https://github.com/guojx0820
@Email : [email protected]
'''
import os
import time
import numpy as np
from osgeo import gdal, osr
class Retrieval_AOD_AE_Classification:
def __init__(self):
self.geo_resolution = 0.01
def _read_tiff_mcd_(self, aod_550_file, aod_470_file):
aod_550_dataset = gdal.Open(aod_550_file)
aod_470_dataset = gdal.Open(aod_470_file)
aod_550_data = aod_550_dataset.ReadAsArray() * 0.001
aod_470_data = aod_470_dataset.ReadAsArray() * 0.001
geo_transform = aod_550_dataset.GetGeoTransform()
projection = aod_550_dataset.GetProjection()
return aod_550_data, aod_470_data, geo_transform, projection
def _ae_470_550_cal_(self, aod_550_data, aod_470_data, geo_transform, projection):
lon_min = geo_transform[0]
lat_max = geo_transform[3]
geo_resolution_row = geo_transform[1]
geo_resolution_col = geo_transform[5]
ae_470_550_data = -(np.log(aod_470_data / aod_550_data)) / (np.log(470 / 550))
return ae_470_550_data, lon_min, lat_max, geo_resolution_row, geo_resolution_col
def _traditional_method_classification_(self, aod550, ae470_550):
tra_type = np.empty((aod550.shape[0], aod550.shape[1]), dtype=np.float32)
for i_raw in range(tra_type.shape[0]):
for j_col in range(tra_type.shape[1]):
aod_temp = np.round(aod550[i_raw, j_col], 3)
ae_temp = np.round(ae470_550[i_raw, j_col], 3)
if aod_temp == 0. and ae_temp == 0.:
tra_type[i_raw, j_col] = 0.
elif aod_temp > 0.5 and ae_temp > 1.0:
tra_type[i_raw, j_col] = 1.
elif aod_temp < 0.5 and ae_temp > 1.0:
tra_type[i_raw, j_col] = 2.
elif aod_temp < 0.5 and ae_temp < 1.0:
tra_type[i_raw, j_col] = 3.
elif aod_temp > 0.5 and ae_temp < 0.7:
tra_type[i_raw, j_col] = 4.
else:
tra_type[i_raw, j_col] = 5.
return tra_type
def _write_tiff_(self, data, lon_min, lat_max, geo_resolution_row, geo_resolution_col, out_name):
if data.ndim == 3:
band_count, rows, cols = data.shape
else:
band_count, (rows, cols) = 1, data.shape
driver = gdal.GetDriverByName('GTiff')
out_raster = driver.Create(out_name, cols, rows, band_count, gdal.GDT_Float32)
out_raster.SetGeoTransform((lon_min, geo_resolution_row, 0, lat_max, 0, geo_resolution_col))
out_raster_SRS = osr.SpatialReference()
out_raster_SRS.ImportFromEPSG(4326)
out_raster.SetProjection(out_raster_SRS.ExportToWkt())
if band_count == 1:
out_raster.GetRasterBand(1).WriteArray(data)
else:
for i_band_count in range(band_count):
out_raster.GetRasterBand(i_band_count + 1).WriteArray(data[i_band_count])
out_raster.FlushCache()
del out_raster
out_raster = None
if __name__ == '__main__':
start_time = time.time()
input_directory = '/mnt/d/Experiments/Aerosol_Classification/Data/MCD19A2/'
output_directory = '/mnt/d/Experiments/Aerosol_Classification/Data/Results/MCD_Cls/'
if os.path.exists(output_directory) == False:
os.makedirs(output_directory)
for root, dirs, files in os.walk(input_directory):
aod550_file_list = [input_directory + i_tif for i_tif in files if
i_tif.endswith('20220328_avr.tif') and i_tif.startswith('aod_055', 0)]
aod470_file_list = [input_directory + i_tif for i_tif in files if
i_tif.endswith('20220328_avr.tif') and i_tif.startswith('aod_047', 0)]
retr_aod_ae = Retrieval_AOD_AE_Classification()
for i_aod550, j_aod470 in zip(aod550_file_list, aod470_file_list):
single_start_time = time.time()
out_name = output_directory + str(os.path.basename(i_aod550)[12:20]) + '_MCD_Cls.tiff'
aod_550_data, aod_470_data, geo_transform, projection = retr_aod_ae._read_tiff_mcd_(i_aod550, j_aod470)
ae_470_550_data, lon_min, lat_max, geo_resolution_row, geo_resolution_col = retr_aod_ae._ae_470_550_cal_(
aod_550_data, aod_470_data, geo_transform, projection)
tra_type = retr_aod_ae._traditional_method_classification_(aod_550_data, ae_470_550_data)
retr_aod_ae._write_tiff_(tra_type, lon_min, lat_max, geo_resolution_row, geo_resolution_col, out_name)
single_end_time = time.time()
single_run_time = np.round(single_end_time - single_start_time, 3)
print('The data ' + str(os.path.basename(i_aod550)[12:20]) + ' are classified, the time cosuming is ' +
str(single_run_time) + ' s.')
end_time = time.time()
total_run_time = np.round(end_time - start_time, 3)
print('All the data are classified, the time cosuming is ' + str(total_run_time) + ' s.')
|
