Coverage for wrainfo/clutter.py: 88%

1"""Clutter detection module.""" 

2 

3# WRaINfo, Is a software to process FURUNO weather radar data. 

4# 

5# Copyright (c) 2022, FernLab (GFZ Potsdam, fernlab@gfz-potsdam.de) 

6# 

7# This software was developed within the context of the RaINfo ("Potential use of 

8# high resolution weather data in agriculture") project of FernLab funded by 

9# the Impulse and Networking Fund of the Helmholtz Association. 

10# 

11# Licensed under the Apache License, Version 2.0 (the "License"); 

12# you may not use this file except in compliance with the License. 

13# 

14# You may obtain a copy of the License at 

15# http://www.apache.org/licenses/LICENSE-2.0 

16# 

17# Unless required by applicable law or agreed to in writing, software 

18# distributed under the License is distributed on an "AS IS" BASIS, 

19# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

20# See the License for the specific language governing permissions and 

21# limitations under the License. 

22 

23# imports 

24# ------- 

25 

26import logging 

27import xarray as xr 

28import wradlib as wrl 

29import numpy as np 

30from wrainfo.reader import read_single_file 

31 

32 

33# fuzzy echo classification based on wRadlib 

34# ------------------------------------------- 

35 

36def fuzzy_echo_classification(ds, cmap, dims=["azimuth", "range"]): 

37 """Fuzzy echo classification and clutter identification based on polarimetric variables. 

38 

39 Parameters 

40 ---------- 

41 ds : xarray.Dataset 

42 cmap : xarray.DataArray 

43 processed clutter map 

44 dims : list 

45 dimension of the dataset 

46 

47 Returns 

48 ------- 

49 : xarray.Dataset 

50 xarray.Dataset with a clutter map. 

51 """ 

52 trpz_default = {"zdr": [0.7, 1.0, 9999, 9999], 

53 "rho": [0.1, 0.15, 9999, 9999], 

54 "phi": [15, 20, 10000, 10000], 

55 "dop": [-0.2, -0.1, 0.1, 0.2], 

56 "map": [1, 1, 9999, 9999], 

57 "rho2": [-9999, -9999, 0.95, 0.98], 

58 "dr": [-20, -12, 9999, 9999], 

59 "cpa": [0.6, 0.9, 9999, 9999], } 

60 

61 weights = {"zdr": 0.4, 

62 "rho": 0.4, 

63 "rho2": 0.4, 

64 "phi": 0.1, 

65 "dop": 0.1, 

66 "map": 0.5} 

67 

68 dat = dict(rho=ds.RHOHV.values, 

69 phi=ds.PHIDP.values, 

70 ref=ds.DBZH.values, 

71 rho2=ds.RHOHV.values, 

72 dop=ds.VRAD.values, 

73 zdr=ds.ZDR.values, 

74 map=cmap.values) 

75 

76 clmap, nanmask = wrl.clutter.classify_echo_fuzzy(dat, 

77 weights=weights, 

78 trpz=trpz_default, 

79 thresh=0.5) 

80 

81 ds = ds.assign(dict(FUZZ=(dims, clmap))) 

82 

83 return ds 

84 

85 

86# remove clutter 

87# -------------- 

88 

89def dbzh_no_clutter(ds, dims=["azimuth", "range"]): 

90 """Set DBZH values, where clutter were identified, to NaN-values. 

91 

92 Parameters 

93 ---------- 

94 ds : xarray.Dataset 

95 dims : list 

96 dimensions of the dataset 

97 

98 Returns 

99 ------- 

100 : xarray.Dataset 

101 xarray.Dataset with removed clutter. 

102 """ 

103 dbzh_no_clutter = ds.DBZH.where(ds.FUZZ == False) # noqa E712 

104 

105 dbzh_no_clutter_attrs = {"standard_name": "radar_equivalent_reflectivity_factor_h", 

106 "long_name": "Equivalent reflectivity factor H", 

107 "unit": "dBZ", } 

108 

109 ds["DBZH_no_clutter"] = xr.DataArray(dbzh_no_clutter, 

110 dims=dims, 

111 attrs=dbzh_no_clutter_attrs) 

112 

113 return ds 

114 

115 

116def create_clutter_map_sequential(files, threshold=0, grp="dataset1", status=False, logger=logging.getLogger()): 

117 """Create a clutter map sequential in order to not overload the working memory. 

118 

119 Parameters 

120 ---------- 

121 files : list 

122 threshold : integer 

123 threshold of the reflectivity 

124 status : bool 

125 percent of progress 

126 

127 Returns 

128 ------- 

129 : xarray.Dataset 

130 xarray.Dataset with clutter values. 

131 """ 

132 file_count = 0 

133 elapsed = 0 

134 logger.info("Processing " + str(len(files)) + " files sequentially.") 

135 file = files[0] 

136 sum_dbzh_values = None 

137 

138 for file in files: 

139 try: 

140 # read furuno files 

141 data = read_single_file(file=file, grp=grp) 

142 

143 except Exception: 

144 # except errors if a file could not read 

145 logger.warning("Could not load file:' " + file + "'.") 

146 continue 

147 

148 # create a data array of reflectivity values 

149 if sum_dbzh_values is None: 

150 sum_dbzh_values = data.DBZH[0, :, :].values 

151 

152 else: 

153 values = data.DBZH[0, :, :].values 

154 

155 # extract only correct values without NaN values 

156 good = ~np.isnan(values) 

157 sum_dbzh_values[good] += data.DBZH[0, :, :].values[good] 

158 

159 file_count += 1 

160 

161 new = round(20 * file_count / len(files), 0) 

162 

163 if status: 

164 if new > elapsed: 

165 elapsed = new 

166 print(str(elapsed * 5) + "% done") 

167 

168 logger.info("Creating boolean clutter map.") 

169 sum_dbzh_bool = sum_dbzh_values > threshold 

170 

171 # set attributes for the new data array 

172 sum_dbzh_values_attrs = {"standard_name": "radar_equivalent_reflectivity_factor_h", 

173 "long_name": "Equivalent reflectivity factor H", 

174 "unit": "dBZ", } 

175 

176 # create dataset 

177 logger.info("Creating arrays.") 

178 data["DBZH_sum_values"] = xr.DataArray(sum_dbzh_values, 

179 dims=["azimuth", "range"], 

180 attrs=sum_dbzh_values_attrs) 

181 

182 data["DBZH_sum_bool"] = xr.DataArray(sum_dbzh_bool, 

183 dims=["azimuth", "range"], 

184 attrs=sum_dbzh_values_attrs) 

185 

186 # delete all other moments in the dataset except for clutter map 

187 variables = ['DBZH_sum_bool', 'DBZH_sum_values'] 

188 cmap = data[variables] 

189 

190 return cmap