卫星数据匹配和分析

python计算脚本

这个程序是用于批量匹配FY4A和MODIS卫星数据的Python脚本。它利用时间和空间上的近似，找到MODIS数据中与FY4A数据中最接近的像素，并计算一些统计指标。以下是详细注释和解释。

原代码

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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
# -*- coding: utf-8 -*-
"""
 This code is to batch collocate the fy4a and MODIS data (distime<=5) 

 by looking for the closest pixel.

 这段代码的主要目标是对FY4A和MODIS卫星数据进行批量匹配，通过查找最接近的像素。代码主要包括以下几个部分：

函数定义：

cal_distance(lat_mod, lat_fy, lon_mod, lon_fy)：计算两个像素之间的距离。
cal_uniform(x, y)：对给定位置附近的MODIS观测值进行采样，计算统计指标。
get_days(year, month, day)：根据给定的年、月、日计算一年中的第几天（Day of Year）。
get_filename(path, filetype)：获取指定路径下指定文件类型的所有文件名。
文件路径和文件名的获取：

获取FY4A的反射率、地理信息、云类别、云顶高度等数据文件路径和文件名。
获取MODIS的反射率、地理信息、物理性质（如云顶高度、云光学厚度、云有效半径等）等数据文件路径和文件名。
数据读取和处理：

读取FY4A的地理信息、反射率等数据。
读取MODIS的地理信息、反射率等数据。
对数据进行预处理，包括坐标系的调整、缺失值的处理等。
数据匹配：

遍历FY4A的地理信息文件，获取日期和时间信息。
遍历MODIS的地理信息文件，获取日期和时间信息。
根据时间和日期信息匹配FY4A和MODIS的数据。
找到匹配的数据后，进行进一步的处理，计算一些统计指标。
一些处理步骤：

处理FY4A的时间信息，筛选出在一定时间范围内的观测。
读取FY4A的地理信息，处理经纬度数据。
读取FY4A的云类别和云顶高度数据，进行一些处理。
读取MODIS的地理信息和物理性质数据，进行一些处理。
读取MODIS的反射率数据，进行一些处理。
数据输出：

在匹配成功的情况下，将相关信息输出到文本文件中。
整体来说，这段代码是一个数据处理和匹配的脚本，用于处理FY4A和MODIS的卫星数据，计算一些统计指标，并输出到文本文件中。
"""

import numpy as np
import datetime
import h5py
import os 
import netCDF4 as nc
from pyhdf import SD
from pyhdf.SD import SD

#计算两个像素之间的距离
def cal_distance(lat_mod,lat_fy,lon_mod,lon_fy):
    dis = (lat_mod - lat_fy)**2 + (lon_mod - lon_fy)**2
    return dis

#对给定位置附近的MODIS观测值进行采样，计算统计指标
def cal_uniform(x,y):
    M1_ref1  = float(mod_ref02[x-2,y-2])
    M1_ref2  = float(mod_ref02[x-2,y-1])
    M1_ref3  = float(mod_ref02[x-2,y])
    M1_ref4  = float(mod_ref02[x-2,y+1])
    M1_ref5  = float(mod_ref02[x-2,y+2])
    M1_ref6  = float(mod_ref02[x-1,y-2])
    M1_ref7  = float(mod_ref02[x-1,y-1])
    M1_ref8  = float(mod_ref02[x-1,y])
    M1_ref9  = float(mod_ref02[x-1,y+1])
    M1_ref10 = float(mod_ref02[x-1,y+2])
    M1_ref11 = float(mod_ref02[x,y-2])
    M1_ref12 = float(mod_ref02[x,y-1])
    M1_ref13 = float(mod_ref02[x,y])
    M1_ref14 = float(mod_ref02[x,y+1])
    M1_ref15 = float(mod_ref02[x,y+2])
    M1_ref16 = float(mod_ref02[x+1,y-2])
    M1_ref17 = float(mod_ref02[x+1,y-1])
    M1_ref18 = float(mod_ref02[x+1,y])
    M1_ref19 = float(mod_ref02[x+1,y+1])
    M1_ref20 = float(mod_ref02[x+1,y+2])
    M1_ref21 = float(mod_ref02[x+2,y-2])
    M1_ref22 = float(mod_ref02[x+2,y-1])
    M1_ref23 = float(mod_ref02[x+2,y])
    M1_ref24 = float(mod_ref02[x+2,y+1])
    M1_ref25 = float(mod_ref02[x+2,y+2])

    M2_ref1  = float(mod_ref07[x-2,y-2])
    M2_ref2  = float(mod_ref07[x-2,y-1])
    M2_ref3  = float(mod_ref07[x-2,y])
    M2_ref4  = float(mod_ref07[x-2,y+1])
    M2_ref5  = float(mod_ref07[x-2,y+2])
    M2_ref6  = float(mod_ref07[x-1,y-2])
    M2_ref7  = float(mod_ref07[x-1,y-1])
    M2_ref8  = float(mod_ref07[x-1,y])
    M2_ref9  = float(mod_ref07[x-1,y+1])
    M2_ref10 = float(mod_ref07[x-1,y+2])
    M2_ref11 = float(mod_ref07[x,y-2])
    M2_ref12 = float(mod_ref07[x,y-1])
    M2_ref13 = float(mod_ref07[x,y])
    M2_ref14 = float(mod_ref07[x,y+1])
    M2_ref15 = float(mod_ref07[x,y+2])
    M2_ref16 = float(mod_ref07[x+1,y-2])
    M2_ref17 = float(mod_ref07[x+1,y-1])
    M2_ref18 = float(mod_ref07[x+1,y])
    M2_ref19 = float(mod_ref07[x+1,y+1])
    M2_ref20 = float(mod_ref07[x+1,y+2])
    M2_ref21 = float(mod_ref07[x+2,y-2])
    M2_ref22 = float(mod_ref07[x+2,y-1])
    M2_ref23 = float(mod_ref07[x+2,y])
    M2_ref24 = float(mod_ref07[x+2,y+1])
    M2_ref25 = float(mod_ref07[x+2,y+2])

    M_COT1   = float(COT_mod_out[x-2,y-2])
    M_COT2   = float(COT_mod_out[x-2,y-1])
    M_COT3   = float(COT_mod_out[x-2,y])
    M_COT4   = float(COT_mod_out[x-2,y+1])
    M_COT5   = float(COT_mod_out[x-2,y+2])
    M_COT6   = float(COT_mod_out[x-1,y-2])
    M_COT7   = float(COT_mod_out[x-1,y-1])
    M_COT8   = float(COT_mod_out[x-1,y])
    M_COT9   = float(COT_mod_out[x-1,y+1])
    M_COT10  = float(COT_mod_out[x-1,y+2])
    M_COT11  = float(COT_mod_out[x,y-2])
    M_COT12  = float(COT_mod_out[x,y-1])
    M_COT13  = float(COT_mod_out[x,y])
    M_COT14  = float(COT_mod_out[x,y+1])
    M_COT15  = float(COT_mod_out[x,y+2])
    M_COT16  = float(COT_mod_out[x+1,y-2])
    M_COT17  = float(COT_mod_out[x+1,y-1])
    M_COT18  = float(COT_mod_out[x+1,y])
    M_COT19  = float(COT_mod_out[x+1,y+1])
    M_COT20  = float(COT_mod_out[x+1,y+2])
    M_COT21  = float(COT_mod_out[x+2,y-2])
    M_COT22  = float(COT_mod_out[x+2,y-1])
    M_COT23  = float(COT_mod_out[x+2,y])
    M_COT24  = float(COT_mod_out[x+2,y+1])
    M_COT25  = float(COT_mod_out[x+2,y+2])

    M_CER1   = float(CER_mod_out[x-2,y-2])
    M_CER2   = float(CER_mod_out[x-2,y-1])
    M_CER3   = float(CER_mod_out[x-2,y])
    M_CER4   = float(CER_mod_out[x-2,y+1])
    M_CER5   = float(CER_mod_out[x-2,y+2])
    M_CER6   = float(CER_mod_out[x-1,y-2])
    M_CER7   = float(CER_mod_out[x-1,y-1])
    M_CER8   = float(CER_mod_out[x-1,y])
    M_CER9   = float(CER_mod_out[x-1,y+1])
    M_CER10  = float(CER_mod_out[x-1,y+2])
    M_CER11  = float(CER_mod_out[x,y-2])
    M_CER12  = float(CER_mod_out[x,y-1])
    M_CER13  = float(CER_mod_out[x,y])
    M_CER14  = float(CER_mod_out[x,y+1])
    M_CER15  = float(CER_mod_out[x,y+2])
    M_CER16  = float(CER_mod_out[x+1,y-2])
    M_CER17  = float(CER_mod_out[x+1,y-1])
    M_CER18  = float(CER_mod_out[x+1,y])
    M_CER19  = float(CER_mod_out[x+1,y+1])
    M_CER20  = float(CER_mod_out[x+1,y+2])
    M_CER21  = float(CER_mod_out[x+2,y-2])
    M_CER22  = float(CER_mod_out[x+2,y-1])
    M_CER23  = float(CER_mod_out[x+2,y])
    M_CER24  = float(CER_mod_out[x+2,y+1])
    M_CER25  = float(CER_mod_out[x+2,y+2])                                    
         
    M1_ref   = [M1_ref1,M1_ref2,M1_ref3,M1_ref4,M1_ref5,M1_ref6,M1_ref7,M1_ref8,M1_ref9,M1_ref10,M1_ref11,M1_ref12,M1_ref13,M1_ref14,M1_ref15,M1_ref16,M1_ref17,M1_ref18,M1_ref19,M1_ref20,M1_ref21,M1_ref22,M1_ref23,M1_ref24,M1_ref25]
    M2_ref   = [M2_ref1,M2_ref2,M2_ref3,M2_ref4,M2_ref5,M2_ref6,M2_ref7,M2_ref8,M2_ref9,M2_ref10,M2_ref11,M2_ref12,M2_ref13,M2_ref14,M2_ref15,M2_ref16,M2_ref17,M2_ref18,M2_ref19,M2_ref20,M2_ref21,M2_ref22,M2_ref23,M2_ref24,M2_ref25]
    Mref_ss1 = np.std(M1_ref)
    Mref_ss2 = np.std(M2_ref)
    Mref_ave = np.sum(M1_ref)/len(M1_ref)
    Hp_ref   = Mref_ss1/Mref_ave                 
                                   
    M_COT    = [M_COT1,M_COT2,M_COT3,M_COT4,M_COT5,M_COT6,M_COT7,M_COT8,M_COT9,M_COT10,M_COT11,M_COT12,M_COT13,M_COT14,M_COT15,M_COT16,M_COT17,M_COT18,M_COT19,M_COT20,M_COT21,M_COT22,M_COT23,M_COT24,M_COT25]
    M_CER    = [M_CER1,M_CER2,M_CER3,M_CER4,M_CER5,M_CER6,M_CER7,M_CER8,M_CER9,M_CER10,M_CER11,M_CER12,M_CER13,M_CER14,M_CER15,M_CER16,M_CER17,M_CER18,M_CER19,M_CER20,M_CER21,M_CER22,M_CER23,M_CER24,M_CER25]
    M_COTss  = np.std(M_COT)
    M_CERss  = np.std(M_CER)
    M_COTave = np.sum(M_COT)/len(M_COT)
    M_CERave = np.sum(M_CER)/len(M_COT)
    MCOT_tau = M_COTss/M_COTave
    MCER_tau = M_CERss/M_CERave

    variable = np.array([Mref_ss1,Mref_ss2,MCOT_tau,MCER_tau,Hp_ref])
    return variable

#根据给定的年、月、日计算一年中的第几天（Day of Year）
def get_days(year,month,day):
    date = datetime.date(year,month,day)
    fy_DOY = int(date.strftime('%j'))
    return fy_DOY

#获取指定路径下指定文件类型的所有文件名
def get_filename(path,filetype):  
    filename = []  
    for root,dirs,files in os.walk(path):  
        for i in files:  
            if filetype in i:  
                filename.append(i.replace(filetype,''))
    return filename


#*************************batch read filename imformation***************************
#time = '202303/'
#ll = 1
#kk = str(ll).zfill(2)
Dir_fy_ref        = 'FY_FDI/' #'/cloud/wanzi/DATAFY4A/FY_FDI/' + time + kk + '/'
Dir_fy_geo        = 'FY_GEO/' #'/cloud/wanzi/DATAFY4A/FY_GEO/' + time + kk + '/'
Dir_fy_clp        = 'FY_CLP/' #'/cloud/wanzi/DATAFY4A/FY_CLP/' + time + kk + '/'
Dir_fy_cth        = 'FY_CTH/' #'/cloud/wanzi/DATAFY4A/FY_CTH/' + time + kk + '/'
Dir_mod_ref       = 'MYD021/' #'/student/home/wanzi/DATAMODIS/MYD021/'+ time + kk + '/'
Dir_mod_geo       = 'MYD03/'  #'/student/home/wanzi/DATAMODIS/MYD03/' + time + kk + '/'
Dir_mod_properties= 'MYD06/'  #'/student/home/wanzi/DATAMODIS/MYD06/' + time + kk + '/'
Dir_mod_sunglint  = 'MYD35/'  #'/student/home/wanzi/DATAMODIS/MYD35/' + time + kk + '/'

filetype_fyL1 = '.HDF'
filetype_fyL2 = '.NC'
filetype_mod  = '.hdf'
filename_Fref = get_filename(Dir_fy_ref,filetype_fyL1)
filename_Fgeo = get_filename(Dir_fy_geo,filetype_fyL1)
filename_Fclp = get_filename(Dir_fy_clp,filetype_fyL2)
filename_Fcth = get_filename(Dir_fy_cth,filetype_fyL2)
filename_Mref = get_filename(Dir_mod_ref,filetype_mod)
filename_Mgeo = get_filename(Dir_mod_geo,filetype_mod)
filename_Mproperties = get_filename(Dir_mod_properties,filetype_mod)
filename_Msunglint   = get_filename(Dir_mod_sunglint,filetype_mod)
filename_Fref.sort()         #File Sorting
filename_Fgeo.sort()
filename_Fclp.sort()
filename_Fcth.sort()
filename_Mref.sort()
filename_Mgeo.sort()
filename_Mproperties.sort()
filename_Msunglint.sort()

for ii in range(len(filename_Fgeo)):
    fy_date  = filename_Fgeo[ii][44:44+12]
    fy_year  = fy_date[0:4]
    fy_month = fy_date[4:6]
    fy_day   = fy_date[6:8]
    fy_time  = fy_date[8:12]
    fy_DOY   = get_days(int(fy_year),int(fy_month),int(fy_day)) 
    
    collocate_name = './collocate_' + fy_year + str(fy_DOY).zfill(3) + fy_time + '.txt'
    file = open(collocate_name,'w')
    
    for jj in range(len(filename_Mgeo)):
        mod_year      = filename_Mgeo[jj][7:11]
        mod_dayofyear = filename_Mgeo[jj][11:14]
        mod_time      = filename_Mgeo[jj][15:19]
        mod_DOY       = int(mod_dayofyear)


#***********************Choose the MODIS Obs Time in AGRI Obs Time Files***********************
        dis_filetime = int(mod_time[0:2])*60+int(mod_time[2:4])-int(fy_time[0:2])*60-int(fy_time[2:4])      
       
        if mod_DOY == fy_DOY and dis_filetime >= 0 and dis_filetime < 15:
            print('MODIS Obs Time in AGRI Obs Time: fy_time',fy_year,fy_month,fy_day,fy_time,'& mod_time',mod_year,mod_dayofyear,mod_time)
            

#*************************Read FY4A Geolocation imformation************************************
            fy_Geolocation = 'Geolocation.hdf'
            FileF1 = SD(fy_Geolocation)
            Lat_fy = FileF1.select('pixel_latitude')[:]
            Lon_fy = FileF1.select('pixel_longitude')[:]
            Lon_fy = np.array(Lon_fy)
            Lon_fy[Lon_fy == -999] = np.nan
            Lat_fy = np.array(Lat_fy)
            Lat_fy[Lat_fy == -999] = np.nan

            loc = np.argwhere(Lon_fy < -160)
            for i in range(len(loc)):
                Lon_fy[loc[i,0],loc[i,1]] = Lon_fy[loc[i,0],loc[i,1]] + 360
            Lon_fy = np.array(Lon_fy)
            Lon_fy[(Lon_fy < 40) | (Lon_fy > 160)] = np.nan
            Lat_fy = np.array(Lat_fy)
            Lat_fy[(Lat_fy < -60) | (Lat_fy > 60)] = np.nan


#*************************Read FY4A L1 Data and Reflectance************************************            
            fy_ref = Dir_fy_ref + filename_Fref[ii] + filetype_fyL1
            FileF2 = h5py.File(fy_ref,'r')
            FileF2.keys()
            reflectance3  = np.copy(FileF2[u'CALChannel03'])[:]
            reflectance6  = np.copy(FileF2[u'CALChannel06'])[:]
            reflectance11 = np.copy(FileF2[u'CALChannel11'])[:]
            reflectance12 = np.copy(FileF2[u'CALChannel12'])[:]
            reflectance13 = np.copy(FileF2[u'CALChannel13'])[:]
            ref0871       = np.copy(FileF2[u'NOMChannel03'])[:]
            ref2251       = np.copy(FileF2[u'NOMChannel06'])[:]
            bt851         = np.copy(FileF2[u'NOMChannel11'])[:]
            bt111         = np.copy(FileF2[u'NOMChannel12'])[:]
            bt121         = np.copy(FileF2[u'NOMChannel13'])[:]
            
            (rows,cols) = ref2251.shape
            ref087 = np.empty(ref2251.shape)
            ref225 = np.empty(ref2251.shape)
            bt85   = np.empty(ref2251.shape)
            bt11   = np.empty(ref2251.shape)
            bt12   = np.empty(ref2251.shape)
            for i in range(rows):
                for j in range(cols):
                    if (ref0871[i,j] < 4096) :
                        ref087[i,j] = reflectance3[ref0871[i,j]]
                    if (ref2251[i,j] < 4096) :
                        ref225[i,j] = reflectance6[ref2251[i,j]]
                    if (bt851[i,j] < 4096):
                         bt85[i,j] = reflectance11[bt851[i,j]]
                    if (bt111[i,j] < 4096):
                         bt11[i,j] = reflectance12[bt111[i,j]]
                    if (bt121[i,j] < 4096):
                         bt12[i,j] = reflectance13[bt121[i,j]]
            
            fy_geo = Dir_fy_geo + filename_Fgeo[ii] + filetype_fyL1
            FileF3 = h5py.File(fy_geo,'r')
            FileF3.keys()
            fy_vaa = np.copy(FileF3[u'NOMSatelliteAzimuth'])[:]
            fy_vza = np.copy(FileF3[u'NOMSatelliteZenith'])[:]
            fy_saa = np.copy(FileF3[u'NOMSunAzimuth'])[:]
            fy_sza = np.copy(FileF3[u'NOMSunZenith'])[:]            
            fy_phi = abs(fy_saa + 180 - fy_vaa)
            
            loc = np.argwhere(fy_phi>180)
            for i in range(len(loc)):
                fy_phi[loc[i,0],loc[i,1]] = abs(360 - fy_phi[loc[i,0],loc[i,1]])
            loc = np.argwhere(fy_sza>89)
            for i in range(len(loc)):
                fy_sza[loc[i,0],loc[i,1]] = abs(fy_sza[loc[i,0],loc[i,1]] - 180)
       

#*************************Read FY4A Phase*********************************************************             
            fy_clp = Dir_fy_clp + filename_Fclp[ii] + filetype_fyL2
            FileF4 = nc.Dataset(fy_clp,'r')
            fy_clp = FileF4.variables['CLP'][:]
            fy_clp = np.array(fy_clp)
            fy_clp[fy_clp>=126] = np.nan
       

#*************************Read FY4A CTH*********************************************************             
            fy_cth = Dir_fy_cth + filename_Fcth[ii] + filetype_fyL2
            FileF5 = nc.Dataset(fy_cth,'r')
            fy_cth = FileF5.variables['CTH'][:]
            fy_cth = np.array(fy_cth)
            fy_cth[fy_cth>=20000] = np.nan
            fy_cth[fy_cth<=1]     = np.nan
            
            
#*************************Read MYD03 Data and Calculate MODIS geo information*********************
            mod_geo  = Dir_mod_geo + filename_Mgeo[jj] + filetype_mod
            FileM1   = SD(mod_geo)
            Lon_mod  = FileM1.select('Longitude')[:]
            Lat_mod  = FileM1.select('Latitude')[:]
            mod_sza0 = FileM1.select('SolarZenith')[:]
            mod_saa0 = FileM1.select('SolarAzimuth')[:]
            mod_vza0 = FileM1.select('SensorZenith')[:]
            mod_vaa0 = FileM1.select('SensorAzimuth')[:]
            
            MODIS_SZAinfo  = FileM1.select('SolarZenith')
            MODIS_SAAinfo  = FileM1.select('SolarAzimuth')
            MODIS_VZAinfo  = FileM1.select('SensorZenith')
            MODIS_VAAinfo  = FileM1.select('SensorAzimuth')
            attributes_sza = MODIS_SZAinfo.attributes()
            attributes_saa = MODIS_SAAinfo.attributes()
            attributes_vza = MODIS_VZAinfo.attributes()
            attributes_vaa = MODIS_VAAinfo.attributes()
            SZA_scales     = attributes_sza['scale_factor']
            SAA_scales     = attributes_saa['scale_factor']
            VZA_scales     = attributes_vza['scale_factor']
            VAA_scales     = attributes_vaa['scale_factor']
            mod_sza        = mod_sza0 * SZA_scales
            mod_saa        = mod_saa0 * SAA_scales
            mod_vza        = mod_vza0 * VZA_scales
            mod_vaa        = mod_vaa0 * VAA_scales
           
            
#*************************Read MYD06 Data and Calculate MODIS Official COT and CER****************               
            mod_properties = Dir_mod_properties + filename_Mproperties[jj] + filetype_mod
            FileM2    = SD(mod_properties)
            CER_mod   = FileM2.select('Cloud_Effective_Radius')[:].astype(float)
            CER_mod16 = FileM2.select('Cloud_Effective_Radius_16')[:].astype(float)
            CER_mod37 = FileM2.select('Cloud_Effective_Radius_37')[:].astype(float)
            COT_mod   = FileM2.select('Cloud_Optical_Thickness')[:].astype(float)
            CTH_mod   = FileM2.select('cloud_top_height_1km')[:].astype(float)
            phase_mod = FileM2.select('Cloud_Phase_Optical_Properties')[:]
            CER_mod[CER_mod < 0]     = np.nan
            CER_mod16[CER_mod16 < 0] = np.nan
            CER_mod37[CER_mod37 < 0] = np.nan
            COT_mod[COT_mod < 0]     = np.nan
            CTH_mod[CTH_mod < 0]     = np.nan

            MODIS_CERinfo  = FileM2.select('Cloud_Effective_Radius')
            attributes_CER = MODIS_CERinfo.attributes()
            CER_scales     = attributes_CER['scale_factor']
            CER_offsets    = attributes_CER['add_offset']
            CER_mod_out    = CER_scales * CER_mod + CER_offsets

            MODIS_COTinfo  = FileM2.select('Cloud_Optical_Thickness')
            attributes_COT = MODIS_COTinfo.attributes()
            COT_scales     = attributes_COT['scale_factor']
            COT_offsets    = attributes_COT['add_offset']
            COT_mod_out    = COT_scales * COT_mod + COT_offsets

            MODIS_CTHinfo  = FileM2.select('cloud_top_height_1km')
            attributes_CTH = MODIS_CTHinfo.attributes()
            CTH_scales     = attributes_CTH['scale_factor']
            CTH_offsets    = attributes_CTH['add_offset']
            CTH_mod_out    = CTH_scales * CTH_mod + CTH_offsets          
 
    
#*************************Read MYD021 Data and Calculate MODIS reflectance*************************                 
            mod_ref   = Dir_mod_ref + filename_Mref[jj] + filetype_mod
            FileM3    = SD(mod_ref)
            mod_ref2  = FileM3.select('EV_250_Aggr1km_RefSB')[:]
            mod_ref7  = FileM3.select('EV_500_Aggr1km_RefSB')[:]
            ref_band2 = mod_ref2[1,:,:].astype(float)
            ref_band7 = mod_ref7[4,:,:].astype(float)
            ref_band2 = np.array(ref_band2)
            ref_band7 = np.array(ref_band7)
            ref_band2[ref_band2 >= 65500] = np.nan
            ref_band7[ref_band7 >= 65500] = np.nan
                           
            MODIS_band2 = FileM3.select('EV_250_Aggr1km_RefSB')
            attributes_band2     = MODIS_band2.attributes()
            reflectance2_scales  = attributes_band2['reflectance_scales']
            reflectance2_offsets = attributes_band2['reflectance_offsets']
            MODIS_band2_scale    = reflectance2_scales[1]
            MODIS_band2_offset   = reflectance2_offsets[1]
            mod_ref02 = MODIS_band2_scale * (ref_band2 - MODIS_band2_offset)
            
            MODIS_band7 = FileM3.select('EV_500_Aggr1km_RefSB')
            attributes_band7     = MODIS_band7.attributes()
            reflectance7_scales  = attributes_band7['reflectance_scales']
            reflectance7_offsets = attributes_band7['reflectance_offsets']
            MODIS_band7_scale    = reflectance7_scales[4]
            MODIS_band7_offset   = reflectance7_offsets[4]
            mod_ref07 = MODIS_band7_scale * (ref_band7 - MODIS_band7_offset)


#*************************Read MYD35 Data and Sunglint Pixels Restore to Clear***********************
            mod_lint   = Dir_mod_sunglint + filename_Msunglint[jj] + filetype_mod
            FileM4     = SD(mod_lint)
            Cloud_Mask = FileM4.select('Cloud_Mask')[:]
            NUM_sunglint_pixels = 0
            sunglint1  = Cloud_Mask[0,:,:]
#            for mm in range(0,Lat_mod.shape[0]):
#                for nn in range(0,Lat_mod.shape[1]):
#                    if sunglint1[mm,nn]>0:
#                        sunglint2 = bin(sunglint1[mm,nn])[2:].zfill(8)
#                        if sunglint2[3] == '0':
#                           NUM_sunglint_pixels = NUM_sunglint_pixels + 1
#                           mod_ref02[mm,nn] = np.nan
#                           mod_ref07[mm,nn] = np.nan  
#                    else:
#                        sunglint2 = bin(sunglint1[mm,nn]+(1<<8))[2:].zfill(8)
#                        if sunglint2[3] == '0':
#                            NUM_sunglint_pixels = NUM_sunglint_pixels + 1
#                            mod_ref02[mm,nn] = np.nan
#                            mod_ref07[mm,nn] = np.nan 


#***********************Read the Line by Line Obs_Time of FY-4A/AGRI***************************
            fy_obstime = np.copy(FileF2[u'NOMObsTime'])[:]
            fy_obstime = np.array(fy_obstime)[:,1]
            fy_linetime= [str(i) for i in fy_obstime]
            fy_linetime= [i[8:12] for i in fy_linetime]

 
#***********************Choose the Line by Line Time Difference Within 5 Minutes*************** 
            dis_linetime = np.zeros(fy_obstime.shape)
            for i in range(len(fy_linetime)):
                if len(fy_linetime[i]) == 0:
                    dis_linetime[i] = np.nan
                else:
                    dis_linetime[i] = abs(int(fy_linetime[i][0:2])*60+int(fy_linetime[i][2:4])-int(mod_time[0:2])*60-int(mod_time[2:4]))
            fy_time_index = dis_linetime <= 5
                  
            for i in range(len(fy_time_index)):
                if fy_time_index[i] == 0:
                    Lat_fy[i][:] = np.nan
                    Lon_fy[i][:] = np.nan
            
#***********************Begin Collocating******************************************************
            min_Lat  = np.nanmin(np.nanmin(Lat_mod))    #MODIS Area Corresponds to the Column of FY-4A
            max_Lat  = np.nanmax(np.nanmax(Lat_mod))
            min_Lon  = np.nanmin(np.nanmin(Lon_mod))
            max_Lon  = np.nanmax(np.nanmax(Lon_mod))
#            [row,col]= np.where((Lat_fy >= min_Lat) & (Lat_fy <= max_Lat) & (Lon_fy >= min_Lon) & (Lon_fy <= max_Lon))
#            min_row  = np.nanmin(row)
#            max_row  = np.nanmax(row)
#            min_col  = np.nanmin(col)
#            max_col  = np.nanmax(col)
#            print(min_row,max_row,min_col,max_col)
                
            FY_loc  = np.argwhere((Lat_fy >= min_Lat) & (Lat_fy <= max_Lat) & (Lon_fy >= min_Lon) & (Lon_fy <= max_Lon))
            for i in range(len(FY_loc)):
                k = FY_loc[i,0];
                l = FY_loc[i,1];
                
                if np.isnan(Lat_fy[k,l]) or np.isnan(Lon_fy[k,l]):
                    continue
                else:
                    distance = cal_distance(Lat_mod,Lat_fy[k,l],Lon_mod,Lon_fy[k,l])
                    mindis   = np.nanmin(np.nanmin(distance))
                    #print(mindis)

                    if mindis <= 1e-3:
                        [m,n]    = np.where(distance == mindis)
                        #print(m[0],n[0])
                        MCER   = CER_mod_out[m[0],n[0]]
                        MCOT   = COT_mod_out[m[0],n[0]]
                        MCTH   = CTH_mod_out[m[0],n[0]]
                        Mref02 = mod_ref02[m[0],n[0]]
                        Mref07 = mod_ref07[m[0],n[0]]
                        if MCER > 0 and MCOT > 0 and MCTH > 0 and Mref02 > 0 and Mref07 > 0:
                            Fsza   = float(fy_sza[k,l])
                            Fsaa   = float(fy_saa[k,l])
                            Fvza   = float(fy_vza[k,l])
                            Fvaa   = float(fy_vaa[k,l])
                            Fphi   = float(fy_phi[k,l])
                            FCTH   = fy_cth[k,l]
                            Fphase = fy_clp[k,l]
                            Fref1  = float(ref087[k,l])
                            Fref2  = float(ref225[k,l])
                            Fbt1   = float(bt85[k,l])
                            Fbt2   = float(bt11[k,l])
                            Fbt3   = float(bt12[k,l])
                            Msza   = float(mod_sza[m[0],n[0]])
                            Msaa   = float(mod_saa[m[0],n[0]])
                            Mvza   = float(mod_vza[m[0],n[0]])
                            Mvaa   = float(mod_vaa[m[0],n[0]])
                            Mphase = phase_mod[m[0],n[0]]

#*******************************Mark Sunglint Pixels(introduction: 1:sunglint_pixel; 0:not sunglint_pixel)******************************
                            if sunglint1[m[0],n[0]] > 0:
                                sunglint2 = bin(sunglint1[m[0],n[0]])[2:].zfill(8)
                                if sunglint2[3] == '0':
                                    sunglint_pixel = 1
                                else:
                                    sunglint_pixel = 0
                            elif sunglint1[m[0],n[0]] < 0:
                                sunglint2 = bin(sunglint1[m[0],n[0]]+(1<<8))[2:].zfill(8)
                                if sunglint2[3] == '0':
                                    sunglint_pixel = 1
                                else:
                                    sunglint_pixel = 0

#*******************************Mark Uniform Pixels(introduction: 1:uniform_pixel; 0:inhomogenous_pixel)******************************
                            if m[0] < 2028 and m[0] > 1 and n[0] < 1352 and n[0] > 1:
                                variable= cal_uniform(m[0],n[0])
                                Hp = variable[4]
                                if variable[0] <= 0.03 and variable[1] <= 0.03 and variable[2] < 0.2 and variable[3] < 0.2:
                                    uniform_pixel = 1
                                else:
                                    uniform_pixel = 0
                            else:
                                uniform_pixel = 0
                                Hp = np.nan

                            file.write('%5d %5d %5d %5d %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %5d %12.4f %12.4f %12.4f %12.4f %5d %12.4f %12.4f %12.2f %12.2f %12.2f %12.4f %12.4f %12.0f %12.0f %5d %5d' % (k,l,m[0],n[0],Lat_fy[k,l],Lon_fy[k,l],Fsza,Fsaa,Fvza,Fvaa,Fphi,Fphase,Msza,Msaa,Mvza,Mvaa,Mphase,Fref1,Fref2,Fbt1,Fbt2,Fbt3,MCOT,MCER,MCTH,FCTH,uniform_pixel,sunglint_pixel))
                            file.write('\n')

    if os.path.getsize(collocate_name) == 0:
        os.remove(collocate_name)

print('Collocate pixels end!')

目录结构

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
.                             # 根目录
│  .DS_Store                  # macOS系统用于存储文件夹自定义属性的隐藏文件
│  collocateCTH.py            # 执行卫星数据匹配的Python脚本文件
│  collocate_20230600545.txt  # 匹配结果输出的文本文件
│  Geolocation.hdf            # FY4A卫星的地理位置信息文件
│
├─FY_CLP                      # FY4A云类别数据目录
│      FY4A-_AGRI--...NC      # FY4A云类别数据文件
│
├─FY_CTH                      # FY4A云顶高度数据目录
│      FY4A-_AGRI--...NC      # FY4A云顶高度数据文件
│
├─FY_FDI                      # FY4A反射率数据目录
│      FY4A-_AGRI--...HDF     # FY4A反射率数据文件
│
├─FY_GEO                      # FY4A地理信息数据目录
│      FY4A-_AGRI--...HDF     # FY4A地理信息数据文件
│
├─MYD021                      # MODIS反射率数据目录
│      MYD021KM...hdf         # MODIS反射率数据文件
│
├─MYD03                       # MODIS地理信息数据目录
│      MYD03...hdf            # MODIS地理信息数据文件
│
├─MYD06                       # MODIS云物理性质数据目录
│      MYD06_L2...hdf         # MODIS云物理性质数据文件
│
└─MYD35                       # MODIS太阳耀斑数据目录
        MYD35_L2...hdf         # MODIS太阳耀斑数据文件

程序解释和注释

在程序的开始部分，首先导入了必要的库，如 numpy、datetime、h5py、os、netCDF4 和 pyhdf。这些库用于数据读取、处理和数学计算。

1
2
3
4
5
6
7
import numpy as np
import datetime
import h5py
import os 
import netCDF4 as nc
from pyhdf import SD
from pyhdf.SD import SD

之后定义了一些函数，用于计算距离、采样统计、获取年内天数和文件名获取。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 计算两个像素之间的距离
def cal_distance(lat_mod, lat_fy, lon_mod, lon_fy):
    # ...

# 对给定位置附近的MODIS观测值进行采样，计算统计指标
def cal_uniform(x, y):
    # ...

# 根据给定的年、月、日计算一年中的第几天（Day of Year）
def get_days(year, month, day):
    # ...

# 获取指定路径下指定文件类型的所有文件名
def get_filename(path, filetype):
    # ...

接着，程序设置了数据文件的路径和类型，并使用 get_filename 函数获取了所有的数据文件名。这些文件名将用于后续的数据匹配和处理。

1
2
3
4
5
6
7
8
# 设置数据文件的路径和类型
Dir_fy_ref = 'FY_FDI/'
Dir_fy_geo = 'FY_GEO/'
# ...

# 获取所有的数据文件名
filename_Fref = get_filename(Dir_fy_ref, filetype_fyL1)
# ...

程序的主要部分是一个双层循环，它迭代遍历FY4A和MODIS的地理信息文件。在每次迭代中，程序检查MODIS的观测时间是否在FY4A观测时间的15分钟内。如果是，它读取FY4A和MODIS的相关数据，包括地理位置、反射率、云顶高度等，并尝试找到最近的MODIS像素与FY4A数据进行匹配。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
for ii in range(len(filename_Fgeo)):
    # FY4A数据处理
    # ...

    for jj in range(len(filename_Mgeo)):
        # MODIS数据处理
        # ...

        # 选择MODIS观测时间在FY4A观测时间内的数据
        if mod_DOY == fy_DOY and dis_filetime >= 0 and dis_filetime < 15:
            # 读取FY4A和MODIS的数据
            # ...

            # 开始卫星数据匹配
            # ...

在匹配过程中，程序检查了太阳耀斑和云的均匀性，并将相关信息写入到一个文本文件中。如果匹配的数据没有产生任何输出，相应的输出文件会被删除。

最后，程序打印出匹配结束的信息。

1
print('Collocate pixels end!')

整个程序的目标是建立FY4A和MODIS数据之间的对应关系，这对于卫星数据的分析和应用非常重要。程序通过时间和空间的近似来实现匹配，其输出可以用于进一步的科学研究或气象分析。

Matlab绘图比较

该MATLAB脚本用于比较AGRI（FY-4A卫星）和MODIS卫星获取的云顶高度（CTH）和云相态（CLP）的数据。它载入了匹配的数据点，并绘制了相关的散点图和概率密度图来显示两个卫星数据集之间的比较结果。下面是对脚本中各部分的解释和注释：

原代码

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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
%% 此程序用于绘制AGRI和MODIS官方CTH空间分布与点对点的对比
clear
clc

%% read data
collocate = load('/Users/wanzi/Desktop/tsetCTH/collocate_20230600545.txt');

MODISfile = '/Users/wanzi/Desktop/tsetCTH/MYD06_L2.A2023060.0550.061.2023132225720.hdf';
MOD_CTH   = double(hdfread(MODISfile,'/cloud_top_height_1km'));
MOD_CTH   = double(MOD_CTH);
MOD_CTH(MOD_CTH == 0 | MOD_CTH == -999)=nan;

AGRIfile_CTH = '/Users/wanzi/Desktop/tsetCTH/FY4A-_AGRI--_N_DISK_1047E_L2-_CTH-_MULT_NOM_20230301054500_20230301055959_4000M_V0001.NC';
AGRI_CTH     = ncread(AGRIfile_CTH,'/CTH');
AGRI_CTH     = AGRI_CTH';
AGRI_CTH(AGRI_CTH > 20000 | AGRI_CTH == -999) = nan;

AGRIfile_CLP = '/Users/wanzi/Desktop/tsetCTH/FY4A-_AGRI--_N_DISK_1047E_L2-_CLP-_MULT_NOM_20230301054500_20230301055959_4000M_V0001.NC';
AGRI_CLP0    = ncread(AGRIfile_CLP,'/CLP');
AGRI_CLP0    = AGRI_CLP0';
AGRI_CLP0(AGRI_CLP0 >= 126) = nan;
[m,n]=size(AGRI_CLP0);
AGRI_CLP=zeros(m,n);
for i=1:m
    for j=1:n
        if AGRI_CLP0(i,j) == 0
            AGRI_CLP(i,j) = 0;     %clear
        end
        if AGRI_CLP0(i,j) == 1 || AGRI_CLP0(i,j) == 2
            AGRI_CLP(i,j) = 1;     %water
        end
        if AGRI_CLP0(i,j) == 3 || AGRI_CLP0(i,j) == 4
            AGRI_CLP(i,j) = 2;     %ice
        end
        if AGRI_CLP0(i,j) == 5
            AGRI_CLP(i,j) = 3;     %uncertain
        end
    end
end

lines1 = 1;
lines2 = length(collocate);
Lat = collocate(lines1:lines2,5);
Lon = collocate(lines1:lines2,6);
loc = collocate(lines1:lines2,1:4);          %AGRI_row, AGRI_col, MODIS_row, MODIS_col
CLP_collocateMODIS = collocate(lines1:lines2,12);
CTH_collocateMODIS = zeros(lines2-lines1+1,1);
CLP_collocateAGRI  = zeros(lines2-lines1+1,1);
CTH_collocateAGRI  = zeros(lines2-lines1+1,1);

for i=1:lines2-lines1+1
    CLP_collocateAGRI(i,1)  = AGRI_CLP(loc(i,1)+1,loc(i,2)+1);
    CTH_collocateAGRI(i,1)  = AGRI_CTH(loc(i,1)+1,loc(i,2)+1);
    CTH_collocateMODIS(i,1) = MOD_CTH(loc(i,3)+1,loc(i,4)+1);
end

CLP_collocateMODIS(CLP_collocateMODIS(:,1)==2) = 1;
CLP_collocateMODIS(CLP_collocateMODIS(:,1)==3) = 2;

MOD_info= hdfinfo(MODISfile);
MOD_COT = double(hdfread(MODISfile,'/Cloud_Optical_Thickness'));
MOD_CER = double(hdfread(MODISfile,'/Cloud_Effective_Radius'));
MOD_COT(MOD_COT == -9999) = 0;
MOD_CER(MOD_CER == -9999) = 0;
M_COT_scale = MOD_info.Vgroup.Vgroup(2).SDS(73).Attributes(5).Value;
M_COT_offset= MOD_info.Vgroup.Vgroup(2).SDS(73).Attributes(6).Value;
M_CER_scale = MOD_info.Vgroup.Vgroup(2).SDS(67).Attributes(5).Value;
M_CER_offset= MOD_info.Vgroup.Vgroup(2).SDS(67).Attributes(6).Value;
MODIS_COT   = MOD_COT.*M_COT_scale+M_COT_offset;
MODIS_CER   = MOD_CER.*M_CER_scale+M_CER_offset;

COT_collocateMODIS = zeros(lines2-lines1+1,1);
COT = collocate(lines1:lines2,20);
for i=1:lines2-lines1+1
    COT_collocateMODIS(i,1) = MODIS_COT(loc(i,3)+1,loc(i,4)+1);
end
figure
scatter(COT_collocateMODIS,COT)

%% merge
% CTH_collocateAGRI_all  = [];
% CTH_collocateMODIS_all = [];
% CLP_collocateMODIS_all =[];
% CLP_collocateAGRI_all =[];
CLP_collocateMODIS_all = [CLP_collocateMODIS_all;CLP_collocateMODIS];
CLP_collocateAGRI_all = [CLP_collocateAGRI_all;CLP_collocateAGRI];
CTH_collocateAGRI_all  = [CTH_collocateAGRI_all;CTH_collocateAGRI];
CTH_collocateMODIS_all = [CTH_collocateMODIS_all;CTH_collocateMODIS];

% Latmin = min(Lat); 
% Latmax = max(Lat); 
% Lonmin = min(Lon); 
% Lonmax = max(Lon); 

%% plot AGRI CLP
% scatter(CLP_collocateAGRI(:,1),CLP_collocateMODIS(:,1))
set(gcf,'units','pixels','position',[50,50,900,800]);

position_ax1 = [0.06,0.55,0.35,0.39];
colorbar_ax1 = [0.43,0.55,0.02,0.39];
ax1 = subplot('Position',position_ax1);
scatter(Lon,Lat,7,CLP_collocateAGRI,'filled')
colormap(ax1,[1 1 1;0 0 1;0.3 0.75 0.93]);
caxis(ax1,[-0.5 2.5]);
hcb = colorbar('position',colorbar_ax1);
set(hcb,'ytick',-1:1:3,'tickdir','out','YTickLabel',{'','Clear','Water','Ice',''},'FontName','Arial','FontSize',15),
set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out');
set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out');
xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ...
    '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ...
    '160°E','165°E','170°E','175°E','180°E'})
yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ...
    '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'})
box on
title('(a) AGRI CPH')
set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);

%% plot MODIS CLP
position_ax4 = [0.06,0.06,0.35,0.39];
colorbar_ax4 = [0.43,0.06,0.02,0.39];
ax4 = subplot('Position',position_ax4);
scatter(Lon,Lat,7,CLP_collocateMODIS,'filled')
colormap(ax4,[1 1 1;0 0 1;0.3 0.75 0.93]);
caxis(ax4,[-0.5 2.5]);
hcb = colorbar('position',colorbar_ax4);
set(hcb,'ytick',-1:1:3,'tickdir','out','YTickLabel',{'','Clear','Water','Ice',''},'FontName','Arial','FontSize',15),
set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out');
set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out');
xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ...
    '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ...
    '160°E','165°E','170°E','175°E','180°E'})
yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ...
    '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'})
box on
title('(c) MODIS CPH')
set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);

%% plot AGRI CTH
INDEX_CTH = 4;
INTERVAL  = 1;
position_ax2 = [0.57,0.55,0.35,0.39];
colorbar_ax2 = [0.94,0.55,0.02,0.39];
ax2 = subplot('Position',position_ax2);
scatter(Lon,Lat,7,CTH_collocateAGRI*0.001,'filled')
% map=mymap(mapname);
% colormap(ax2,map)
hcb=colorbar('position',colorbar_ax2);
caxis(ax2,[0 INDEX_CTH])
set(hcb,'ytick',0:INTERVAL:20,'tickdir','out','FontName','Arial','FontSize',15)
set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out');
set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out');
xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ...
    '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ...
    '160°E','165°E','170°E','175°E','180°E'})
yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ...
    '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'})
box on
title('(b) AGRI CTH (km)')
set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);

%% plot MODIS CTH
position_ax5 = [0.57,0.06,0.35,0.39];
colorbar_ax5 = [0.94,0.06,0.02,0.39];
ax5 = subplot('position',position_ax5);
scatter(Lon,Lat,7,CTH_collocateMODIS*0.001,'filled')
% map=mymap(mapname);
% colormap(ax5,map)
hcb=colorbar('position',colorbar_ax5);
caxis(ax5,[0 INDEX_CTH])
set(hcb,'ytick',0:INTERVAL:20,'tickdir','out','FontName','Arial','FontSize',15)
set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out');
set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out');
xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ...
    '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ...
    '160°E','165°E','170°E','175°E','180°E'})
yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ...
    '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'})
box on
title('(d) MODIS CTH (km)')
set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);

%% plot compare
clear
load('./testCTH.mat')
CTH_collocateMODIS_all(CLP_collocateMODIS_all ~= CLP_collocateAGRI_all,:)=[];
CTH_collocateAGRI_all(CLP_collocateMODIS_all ~= CLP_collocateAGRI_all,:)=[];

% scatter(CTH_collocateAGRI_all*0.001,CTH_collocateMODIS_all*0.001)

pixel_num=length(CTH_collocateMODIS_all);
n=100;
xmin = 0; xmax = 18;
ymin = 0; ymax = 18;
x = CTH_collocateMODIS_all*0.001;  y = CTH_collocateAGRI_all*0.001;
[ Nbin ] = Prob( x, y, xmin, xmax, ymin, ymax, n);
pdata = Nbin'./pixel_num;
max=max(max(pdata));
min=min(min(pdata));
pdata1 = (pdata - min)/(max-min);
pdata1(pdata1==0)=nan;
x_lable=xmin:0.18:xmax; y_lable=ymin:0.18:ymax;

pdata1(pdata1<0)=nan;
pcolor(x_lable(1,1:end-1),y_lable(1,1:end-1),pdata1)
shading flat
set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7); % set figure
set(gca,'xlim',[0 18],'xtick',0:3:50,'tickdir','out');
set(gca,'ylim',[0 18],'ytick',0:3:50,'tickdir','out');
xlabel('Aqua/MODIS  CTH (km)') %Annual
ylabel('FY-4A/AGRI  CTH (km)') %Annual
axis square
hcb=colorbar;
caxis([0 1])
set(hcb,'ytick',0:0.2:1,'tickdir','out')
set(hcb,'TickLabels',{'0.0','0.2','0.4','0.6','0.8','1.0'})
hold on
x1 = [0 100];
y1 = [0 100];
plot(x1,y1,'--','color','black','LineWidth',1.7);


目录结构

1
2
3
4
5
collocate_20230600545.txt
plot_cth_compare.m
plot_cth_compare_cp.m
Prob.m
TestCTH.zip

脚本解读

数据读取

脚本首先清理工作空间，然后读取之前匹配好的AGRI和MODIS数据点信息，包括云顶高度和云相态等。

1
2
3
4
5
clear
clc

% Load the collocated data points between AGRI and MODIS
collocate = load('/Users/wanzi/Desktop/tsetCTH/collocate_20230600545.txt');

数据处理

读取MODIS和AGRI的CTH数据，并将异常值替换为NaN，方便后续处理。AGRI的云相态数据也进行了读取和分类处理。

1
2
3
4
5
6
7
8
9
% Read and preprocess MODIS CTH data
MODISfile = '/Users/wanzi/Desktop/tsetCTH/MYD06_L2.A2023060.0550.061.2023132225720.hdf';
MOD_CTH   = double(hdfread(MODISfile,'/cloud_top_height_1km'));
MOD_CTH(MOD_CTH == 0 | MOD_CTH == -999)=nan;

% Read and preprocess AGRI CTH and CLP data
AGRIfile_CTH = '/Users/wanzi/Desktop/tsetCTH/FY4A-_AGRI--_N_DISK_1047E_L2-_CTH-_MULT_NOM_20230301054500_20230301055959_4000M_V0001.NC';
AGRI_CTH     = ncread(AGRIfile_CTH,'/CTH');
AGRI_CTH(AGRI_CTH > 20000 | AGRI_CTH == -999) = nan;

提取匹配点数据

对于匹配点数据（即AGRI和MODIS的共同观测点），提取CTH和CLP信息，为绘制散点图做准备。

1
2
3
4
5
6
7
8
9
10
11
12
13
% Initialize variables for collocated CTH and CLP data
CTH_collocateMODIS = zeros(lines2-lines1+1,1);
CLP_collocateAGRI  = zeros(lines2-lines1+1,1);
CTH_collocateAGRI  = zeros(lines2-lines1+1,1);

% Extract CTH and CLP data for collocated points
for i=1:lines2-lines1+1
    % ...
end

% Modify CLP data for consistency
CLP_collocateMODIS(CLP_collocateMODIS(:,1)==2) = 1;
CLP_collocateMODIS(CLP_collocateMODIS(:,1)==3) = 2;

绘图

脚本使用MATLAB的绘图功能来展示AGRI和MODIS的云相态和云顶高度数据。这些绘图包括散点图和概率密度图，用于可视化两个数据集之间的差异和相似性。

1
2
3
4
5
6
7
8
% Plotting AGRI and MODIS CLP in subplots
% ...

% Plotting AGRI and MODIS CTH in subplots
% ...

% Plotting the probability density of collocated CTH between AGRI and MODIS
% ...

概率密度函数

使用自定义的 Prob 函数来计算和绘制AGRI和MODIS CTH数据点的概率密度图。

1
2
3
4
% Calculate and plot the probability density function for collocated CTH data
[ Nbin ] = Prob( x, y, xmin, xmax, ymin, ymax, n);
pdata = Nbin'./pixel_num;
% ...

概率密度图归一化

归一化概率密度图的数值，移除概率为零的点，并设置图的颜色范围和刻度。

1
2
3
4
% Normalize and prepare the data for plotting
pdata1 = (pdata - min)/(max-min);
pdata1(pdata1==0)=nan;
% ...

绘制一致性线

在概率密度图上绘制一条对角线（一致性线），表示AGRI和MODIS CTH数据的一致性。

1
2
3
4
% Plot the consistency line on the probability density plot
x1 = [0 100];
y1 = [0 100];
plot(x1,y1,'--','color','black','LineWidth',1.7);

整个脚本的目的是可视化AGRI和MODIS卫星数据之间的比较，通过这些图形用户可以更直观地理解两个数据集的相似性和差异性。