get_many_data.py
11 KB
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
import urllib.request
import json
import os
import datetime
from imageio import imread
import cv2
import numpy as np
import loca_data
png_path = r"d:\PYTHON\tests\test1\png"
lenght_lan = 79497.3714882
lenght_lon = 111162.6
base_pixel_len = 152.8740566
str_limiter = '||'
d1plan_g = base_pixel_len / lenght_lan
d1plon_g = base_pixel_len / lenght_lon
x286lon = 38.72287
x286lan = 44.417242
base_lon = x286lon - (128 * d1plon_g)
min_lon = x286lon - (128 * d1plon_g)
max_lon = x286lon + (128 * d1plon_g)
base_lan = x286lan + (128 * d1plan_g)
min_lan = x286lan - (128 * d1plan_g)
max_lan = x286lan + (128 * d1plan_g)
pol_store = []
prec_meas = []
# Генерировать измерители всех типов
def GenerateMeasurements():
#l_dlan = d1plan_g * 48
#l_dlon = d1plon_g * 48
# Точки для полигонов
result = []
result.append(GetMeas(0, 'TEST_00', x286lan, x286lon))
result.append(GetMeas(1, 'TEST_01', 44.3251277777778, 38.6788694444444))
result.append(GetMeas(2, 'TEST_02', 44.3727555555556, 38.7005))
result.append(GetMeas(3, 'TEST_03', 44.4186305555556, 38.6651361111111))
result.append(GetMeas(4, 'TEST_04', 44.4545888888889, 38.6328638888889))
result.append(GetMeas(5, 'TEST_05', 44.4508222222222, 38.6833333333333))
result.append(GetMeas(6, 'TEST_06', 44.4329277777778, 38.7327722222222))
result.append(GetMeas(7, 'TEST_07', 44.4106166666667, 38.7633277777778))
result.append(GetMeas(8, 'TEST_08', 44.3846138888889, 38.7533722222222))
result.append(GetMeas(9, 'TEST_09', 44.3416583333333, 38.735175))
prec_meas.append(GetMeas(5, 'X286|s0_d0', x286lan+0, x286lon-0))
pol_store.append(GeneratePolygons([result[0], result[3], result[4], result[5], result[6]], 1, "Verh"))
pol_store.append(GeneratePolygons([result[0], result[6], result[7], result[8], result[9], result[1], result[2]], 2, "Nige_Polkovnichego"))
return
# получить измеритель (имя, номер, х-коорд картинки и у-коорд картинки) из широты и долготы
def GetMeas(_number, _name, _lan, _lon):
if _lan <= min_lan or _lan >= max_lan:
print("MeasUnit named {} have not current lan ({})".format(_name, _lan))
return None
if _lon <= min_lon or _lon >= max_lon:
print("MeasUnit named {} have not current lon ({})".format(_name, _lon))
return None
lan_delta = _lan - base_lan
lon_delta = _lon - base_lon
y_delta = lan_delta * lenght_lan
x_delta = lon_delta * lenght_lon
my = y_delta // base_pixel_len
mx = x_delta // base_pixel_len
result = dict(number = _number, name = _name, x = abs(mx), y = abs(my))
return result
# Получить осадки по значению dbZ
def GetPrecByDBZ_2(_dbz):
result = 0
if _dbz == 0 or _dbz > 127:
return result
step1 = (_dbz - 32) / 10
step2 = 10 ** step1
step3 = step2 / 200
step4 = step3 ** 0.625
result = step4
return result
# Сгенерировать полигон по набору геокоординат
def GeneratePolygons(_meas_list, _number, _name):
# Максимальные значения координат Х и У для полигона
max_x = 0
max_y = 0
min_x = 256
min_y = 256
#Временный полигон
tmp_pol = []
# Перебрать измерители последовательно и сформировать полигон как набор пар (х, у) плюс данные о границах полигона (прямоугольных)
for imeas in _meas_list:
mx =imeas["x"]
my = imeas["y"]
# Ищем прямоугольные границы полигона
if mx > max_x:
max_x = mx
if my > max_y:
max_y = my
if mx < min_x:
min_x = mx
if my < min_y:
min_y = my
tmp_pol.append((mx, my))
result = dict(number = _number, name = _name, max_x = max_x, max_y = max_y, min_x = min_x, min_y = min_y, pol = tmp_pol)
result["pol"] = tmp_pol
return result
# Получить "обрезанную" матрицу по полигону
def GetCuttedMatrix(_image, _mpolygon):
matrix = _mpolygon["pol"]
max_x = int(_mpolygon["max_x"])
max_y = int(_mpolygon["max_y"])
min_x = int(_mpolygon["min_x"])
min_y = int(_mpolygon["min_y"])
mask = np.zeros(_image.shape, dtype=np.uint8)
roi_corners = np.array([matrix], dtype=np.int32)
channel_count = _image.shape[2] # i.e. 3 or 4 depending on your image
ignore_mask_color = (255,)*channel_count
cv2.fillPoly(mask, roi_corners, ignore_mask_color)
masked_image = cv2.bitwise_and(_image, mask)
masked_image = masked_image[min_y:max_y, min_x:max_x]
return masked_image
# Получаем все файлы из директории с сохранёнными изображениями
png_files = os.listdir(png_path)
print(png_files)
print("Lan border = {} - {}".format(min_lan, max_lan))
print("Lon border = {} - {}".format(min_lon, max_lon))
# Сформировать массив "матрица картинки-время" (время из названия файла)
store = []
for png_f in png_files:
dt_value = png_f.replace('dj_286m_', '').replace('.png', '')
dt_value = int(dt_value)
dt_value = datetime.datetime.fromtimestamp(dt_value)
date_diff = datetime.timedelta(hours=0)
dt_value = dt_value + date_diff
file_data = imread(png_path + '\\' + png_f)
st_unit = dict(
datetime = dt_value,
matrix = file_data
)
store.append(st_unit)
# Генерируем измерители
GenerateMeasurements()
# Формируем данные
main_data = "Дата||АГК-286М||АГК-34||РАД-286М||RAD286(cor)||Сумма верх||Сумма верх(кор)||Объём верх||Объём верх(кор)||Сумма низ||Сумма низ(кор)||Объём низ||Объём низ(кор)\n"
# Последниё интервал расчёта
last_datetime = store[0]["datetime"]
# Дельта (служебная) между интервалами расчёта
dt = datetime.timedelta(minutes=10)
# Интервал времени
datetime_min = datetime.datetime(2020, 1, 29, 22, 50, 0)
datetime_max = datetime.datetime(2020, 2, 4, 15, 40, 0)
# Реальные 286М и 34
real286 = loca_data.GetPrec286M()
real34 = loca_data.GetLevel34()
# Пробегаемся по всем датам в хранилище реальных осадков
for rdata_key in real286.keys():
koef_multy_re_d_rad = 1
koef_addit_re_m_rad = 0
prec_data = None
lmatrix = None
ldatetime = rdata_key
print(ldatetime)
local_data = ""
local_data += ldatetime.__str__() + str_limiter
# 0. Осадки реальные
real_prec = real286[rdata_key]
local_data += real_prec.__str__() + str_limiter
# 0.1 Уровень реальный реальные
real_level = real34[rdata_key]
local_data += real_level.__str__() + str_limiter
for imgmatrix in store:
if imgmatrix["datetime"] == rdata_key:
lmatrix = imgmatrix["matrix"]
if lmatrix is None:
local_data += "0||0||0||0||0||0||0||0||0||0||"
main_data += local_data + "\n"
continue
# 1. Расчёт осадков по координатам
# Формируем для каждого измерителя оасдки в точке, согласно его прямоугольным координатам относительно картинки
mes = prec_meas[0]
y = int(mes["y"])
x = int(mes["x"])
value = lmatrix[y][x][0]
dmrl_prec = GetPrecByDBZ_2(value)
dmrl_prec = dmrl_prec / 6
local_data += dmrl_prec.__str__() + str_limiter
# 1.1 Считаем коэффициенты
if dmrl_prec > 0 and real_prec > 0:
koef_multy_re_d_rad = real_prec / dmrl_prec
koef_addit_re_m_rad = 0
elif dmrl_prec == 0 and real_prec > 0:
koef_multy_re_d_rad = 1
koef_addit_re_m_rad = real_prec - dmrl_prec
elif dmrl_prec > 0 and real_prec == 0:
koef_multy_re_d_rad = 1
koef_addit_re_m_rad = 0
#local_data += koef_multy_re_d_rad.__str__() + str_limiter
#local_data += koef_addit_re_m_rad.__str__() + str_limiter
# 1.2 Считаем осадки после корректировки
dmrl_prec_cor = (dmrl_prec * koef_multy_re_d_rad) + koef_addit_re_m_rad
local_data += dmrl_prec_cor.__str__() + str_limiter
# 2. Расчет площадей и осадков на площадях
for pol in pol_store:
masked_image = GetCuttedMatrix(lmatrix, pol)
sum_prec = 0
sum_squere = 0
sum_prec_cor = 0
vol_water = 0
vol_water_cor = 0
for yy in masked_image:
for xx in yy:
lm_value = xx[0]
m_prec = GetPrecByDBZ_2(lm_value)
m_prec = m_prec / 6
if m_prec > 0:
sum_squere += 1
sum_prec += m_prec
m_prec_cor = ((m_prec * koef_multy_re_d_rad) + koef_addit_re_m_rad)
sum_prec_cor = sum_prec_cor + m_prec_cor
vol_local = base_pixel_len * base_pixel_len * 1/1000 * m_prec
vol_water += vol_local
vol_local_cor = base_pixel_len * base_pixel_len * 1/1000 * m_prec_cor
vol_water_cor += vol_local_cor
local_data += sum_prec.__str__() + str_limiter
local_data += sum_prec_cor.__str__() + str_limiter
local_data += vol_water.__str__() + str_limiter
local_data += vol_water_cor.__str__() + str_limiter
#local_data += sum_squere.__str__() + str_limiter
# Вываливаем данные
main_data += local_data + "\n"
# Заменяем временные разделители и прочие символы для правильного отобраежния в Эксель
main_data = main_data.replace("||", "\t")
#main_data = main_data.replace(".", ",")
# Пишем результат в текстовый файл
f = open('result.txt', 'w')
f.write(main_data)
f.close()
"""
1. Для каждого измерителя ищем разницу в широте и долготе лот базовых координат
1.1 Проверяем на выход за предельные значения
2. По здначениям дельты и длинам градусов ищем разницу в метрах
3. По разнице в метрах и базовому размеру пикселя ищем нужный квадрат как количество целых делений дельты в метрах на базовую длину градуса
"""