-
Notifications
You must be signed in to change notification settings - Fork 1
/
gfssi_b02_ssi_1km.py
156 lines (134 loc) · 5.92 KB
/
gfssi_b02_ssi_1km.py
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time : 2019/8/27
@Author : AnNing
"""
import os
import h5py
import numpy as np
from scipy.interpolate import griddata
from lib.lib_read_ssi import FY4ASSI
from lib.lib_database import add_result_data, exist_result_data
from lib.lib_constant import FULL_VALUE
def lat2row(lat):
row = int(((lat - 9.995) / 0.01))
return row
def modiGHI(a, b, r):
c = a * (1 + (r[0] * b / 1000 + r[1]) * 0.01)
return c
def topoCorrection(radiaArray, deltHgt):
ghi_ri = []
rr = [[2.6036, 0.0365], [2.6204, 0.0365], [2.6553, 0.0362], [2.6973, 0.0356], [2.7459, 0.0343],
[2.8012, 0.0324], [2.8616, 0.0299], [2.9236, 0.0257], [2.9870, 0.0204]]
if len(deltHgt) == len(radiaArray):
for i in range(len(deltHgt)):
if i >= lat2row(52.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[8]))
if i >= lat2row(47.5) and i < lat2row(52.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[7]))
if i >= lat2row(42.5) and i < lat2row(47.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[6]))
if i >= lat2row(37.5) and i < lat2row(42.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[5]))
if i >= lat2row(32.5) and i < lat2row(37.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[4]))
if i >= lat2row(27.5) and i < lat2row(32.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[3]))
if i >= lat2row(22.5) and i < lat2row(27.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[2]))
if i >= lat2row(17.5) and i < lat2row(22.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[1]))
if i < lat2row(17.5):
ghi_ri.append(modiGHI(np.array(radiaArray[i]), np.array(deltHgt[i]), rr[0]))
return np.array(ghi_ri)
def _write_out_file(out_file, result):
valid_count = 0
for key in result:
if result[key] is None:
continue
else:
valid_count += 1
if valid_count == 0:
print('没有足够的有效数据，不生成结果文件')
return
out_dir = os.path.dirname(out_file)
if not os.path.isdir(out_dir):
os.makedirs(out_dir)
try:
compression = 'gzip'
compression_opts = 5
shuffle = True
with h5py.File(out_file, 'w') as hdf5:
for dataset in result.keys():
data = result[dataset]
if data is not None:
data[np.isnan(data)] = FULL_VALUE
hdf5.create_dataset(dataset,
dtype=np.float32, data=result[dataset], compression=compression,
compression_opts=compression_opts,
shuffle=shuffle)
print('>>> 成功生成HDF文件{}'.format(out_file))
except Exception as why:
print(why)
print('HDF写入数据错误')
os.remove(out_file)
def fy4a_ssi_4km_to_1km(in_file, out_file, resultid=None, planid=None, datatime=None, resolution_type=None):
print('<<< itcal: {}'.format(in_file))
area_type = 'Full_DISK'
if os.path.isfile(out_file):
print('数据已经存在: {}'.format(out_file))
if not exist_result_data(resultid=resultid, datatime=datatime,
resolution_type=resolution_type,
area_type=area_type):
add_result_data(resultid=resultid, planid=planid, address=out_file, datatime=datatime,
resolution_type=resolution_type, area_type=area_type)
return
datas = FY4ASSI(in_file)
data_get = {
'SSI': datas.get_ssi,
'DirSSI': datas.get_ib,
'DifSSI': datas.get_id,
'G0': datas.get_g0,
'Gt': datas.get_gt,
'DNI': datas.get_dni,
}
result = {}
elements = data_get.keys()
lats_4km = FY4ASSI.get_latitude_4km()
lons_4km = FY4ASSI.get_longitude_4km()
lats_1km = FY4ASSI.get_latitude_1km()
lons_1km = FY4ASSI.get_longitude_1km()
ddem = FY4ASSI.get_ddem_1km()
lats_min = np.nanmin(lats_1km) - 5
lats_max = np.nanmax(lats_1km) + 5
lons_min = np.nanmin(lons_1km) - 5
lons_max = np.nanmax(lons_1km) + 5
print(lats_min, lats_max, lons_min, lons_max)
index1 = np.logical_and.reduce((lons_4km <= lons_max, lons_4km >= lons_min,
lats_4km <= lats_max, lats_4km >= lats_min,))
for element in elements:
print(element)
values = data_get.get(element)()
index2 = np.isfinite(values)
index = np.logical_and(index1, index2)
valid_count = index.sum()
print(index1.sum())
print(index2.sum())
print('有效点数量:{}'.format(valid_count))
if valid_count <= 0:
data = np.full_like(lons_1km, np.nan)
else:
lats_ = lats_4km[index].reshape(-1, 1)
lons_ = lons_4km[index].reshape(-1, 1)
points = np.concatenate((lons_, lats_), axis=1)
values = values[index]
data = griddata(points, values, (lons_1km, lats_1km), method='linear')
data = topoCorrection(data, ddem)
result[element] = data
_write_out_file(out_file, result)
if os.path.isfile(out_file) and not exist_result_data(resultid=resultid, datatime=datatime,
resolution_type=resolution_type,
area_type=area_type):
add_result_data(resultid=resultid, planid=planid, address=out_file, datatime=datatime,
resolution_type=resolution_type, area_type=area_type)