Coverage for wrainfo/process_chains.py: 69%
109 statements
« prev ^ index » next coverage.py v6.5.0, created at 2022-10-18 12:40 +0000
1"""Process chains module."""
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.
23import logging
24import os
25import glob
26import xarray as xr
27import pathlib
28import datetime as dt
29import numpy as np
30import wrainfo
33# clutter map will be calculated monthly every first day of the month
34# -------------------------------------------------------------------
36def clutter_chain(start_time,
37 path,
38 elevation_angle="dataset1",
39 days=90,
40 threshold=0,
41 status=True,
42 pattern="_000.scnx.gz",
43 data_type="cmap",
44 logger=logging.getLogger()):
45 """Create a ground clutter map and save NetCDF-file.
47 Parameters
48 ----------
49 start_time : datetime.datetime
50 start time which files will be read
51 path : str
52 path to configuration file
53 elevation_angle : str
54 choose dataset of one elevation angle
55 (for example: dataset 1 is elevation angle 0.5°)
56 days : integer
57 Number of days for which a cmap is created
58 threshold : float
59 this is a threshold for the DBZH it should not be negative
60 status : bool
61 percentage of the finished clutter map
62 pattern : str
63 extension of the scnx files
64 (scnx file: elevation angle 0.5° = "_000.scnx")
65 data_type : str
66 the clutter map (cmap) were created
67 logger : file
68 writing status messages to a file
70 Returns
71 -------
72 : files
73 NetCDF - files in output directory.
74 """
75 # read output directory for clutter map from configuration file
76 clutter_directory = wrainfo.reader.read_config_file(path=path,
77 selection="monthly_clutter_directory")
78 radar_location_identifier = wrainfo.reader.read_config_file(path=path,
79 selection="radar_location_identifier")
81 # create the stop time
82 # -----------------------
84 stop_time = start_time + dt.timedelta(days=days)
86 # create a list of all input files for the defined time period
87 # ---------------------------------------------------------------
89 flist = wrainfo.reader.create_filelist(path=path,
90 starttime=start_time,
91 endtime=stop_time,
92 pattern=pattern)
94 # create clutter map
95 # ---------------------
96 logger.info("Creating cluttermap.")
97 cmap = wrainfo.clutter.create_clutter_map_sequential(files=flist,
98 threshold=threshold,
99 grp=elevation_angle,
100 status=status,
101 logger=logger.getChild("wr_furuno.clutter"
102 ".create_clutter_map_sequential"))
104 # elevation angle
105 # ------------------
106 el = cmap.elevation.values[0]
108 # create outfilename
109 # ---------------------
110 outfilename = f"{radar_location_identifier}_{start_time:%Y%m%d}_{stop_time:%Y%m%d}_elev_{el:0.1f}_{data_type}.nc"
111 year = str(stop_time.year)
112 path1 = clutter_directory + "/" + year
114 if not os.path.exists(path1):
115 os.makedirs(path1, exist_ok=True)
116 # overwrite/remove if exist
117 f = pathlib.Path(outfilename)
118 f.unlink(missing_ok=True)
120 outfilename = os.path.join(path1, outfilename)
122 # output file as netcdf
123 # ------------------------
124 logger.info("Creating cluttermap netcdf output.")
125 print(f"-- output to {outfilename}")
126 cmap.to_netcdf(outfilename)
127 logger.info("Finished cluttermap chain.")
129 return True
132# static clutter map will be processed each month from 3 last cluttermaps
133# -----------------------------------------------------------------------
135def static_cmap(path,
136 pattern="_elev_0.5_cmap.nc"):
137 """Create a static ground clutter map and save NetCDF-file.
139 Parameters
140 ----------
141 path : str
142 path to configuration file
143 pattern : str
144 extension of the netcdf file
145 (cmaps for different elevation angles)
147 Returns
148 -------
149 : files
150 NetCDF - files in output directory.
151 """
152 # read path of clutter directory from configuration file
153 path_cmap = wrainfo.reader.read_config_file(path=path,
154 selection="monthly_clutter_directory")
155 path_cmap = path_cmap + "*"
157 # read radar location identifier from configuration file
158 radar_location_identifier = wrainfo.reader.read_config_file(path=path,
159 selection="radar_location_identifier")
161 # read output directory from configuration file
162 outdir = wrainfo.reader.read_config_file(path=path,
163 selection="static_clutter_directory")
165 # load the last 3 cmaps
166 cmaps = sorted(glob.glob(os.path.join(path_cmap, "*")))
168 # filter for elevation angle
169 cmaps_filter = []
171 for cmap in cmaps:
172 if cmap.endswith(pattern):
173 cmaps_filter.append(cmap)
174 cmaps_filter = cmaps_filter[-3:]
176 # read the last 3 cmaps
177 cmap1 = xr.open_dataset(cmaps_filter[0])
178 cmap2 = xr.open_dataset(cmaps_filter[1])
179 cmap3 = xr.open_dataset(cmaps_filter[2])
181 # add cmaps
182 cmap_final = cmap1 + cmap2 + cmap3
184 # outfilename
185 base_name = os.path.basename(cmaps_filter[1])
186 el = base_name[28:36]
188 outfilename = f"{radar_location_identifier}_static_cluttermap_{el}.nc"
190 # output directory
191 if not os.path.exists(outdir):
192 os.makedirs(outdir, exist_ok=True)
193 # overwrite/remove if exist
194 f = pathlib.Path(outfilename)
195 f.unlink(missing_ok=True)
197 outfilename = os.path.join(outdir, outfilename)
199 print(f"-- output to {outfilename}")
200 cmap_final.to_netcdf(outfilename)
202 return True
205# weather radar routine
206# ---------------------
208def wr_routine_furuno(starttime,
209 endtime,
210 delta,
211 path,
212 data_type="Level2a",
213 elevation_angle="dataset1",
214 pattern="_000.scnx.gz",
215 logger=logging.getLogger()):
216 """Process of Furuno raw data.
218 Parameters
219 ----------
220 starttime : datetime.datetime
221 endtime : datetime.datetime
222 delta : datetime.delta
223 path : str
224 path to configuration file
225 data_type : str
226 Level of processing
227 pattern : extension of the scnx file
228 (scnx file: elevation angle 0.5° = "_000.scnx")
229 elevation_angle : str
230 group of hdf5 dataset
231 logger : file
232 writing status messages to a file
234 Returns
235 -------
236 : files
237 NetCDF - files in output directory
238 """
239 date = starttime
240 stop_time = endtime
242 # read settings from configuration file
243 re_idx = wrainfo.reader.read_config_file(path=path,
244 selection="re_index_parameters")
246 re_index_parameters = np.arange(re_idx[0], re_idx[1], re_idx[2])
248 dims = wrainfo.reader.read_config_file(path=path,
249 selection="dimensions")
251 while date < stop_time:
253 # read static cmap
254 # ----------------
256 logger.info("Read static clutter map")
257 clutter_dircetory = wrainfo.reader.read_config_file(path=path,
258 selection="static_clutter_directory")
259 cmap = xr.open_dataset(clutter_dircetory)
260 cmap = cmap.DBZH_sum_bool
262 # load files
263 # ----------
264 logger.info("Create list of files to be processed")
265 flist_raw = wrainfo.reader.create_filelist(path=path,
266 starttime=date,
267 endtime=date + delta,
268 pattern=pattern)
270 # start processing
271 # ----------------
272 for file in flist_raw:
274 # extension of the file
275 # ---------------------
276 extension = os.path.splitext(file)[1]
278 if extension == ".gz":
280 # read Furuno scnx data
281 # ---------------------
282 logger.info("Read Furuno scnx file")
283 try:
284 ds = xr.open_dataset(file,
285 engine="furuno",
286 group=1,
287 backend_kwargs=dict(reindex_angle=1.0))
289 except Exception:
290 # except errors if a file could not read
291 logger.warning("Could not load scnx file:' " + file + "'.")
292 continue
294 # rename VRADH to VRAD
295 ds['VRAD'] = ds['VRADH']
296 ds = ds.drop(['VRADH'])
298 if extension == ".h5":
300 # read Furuno hdf5 data
301 # ---------------------
302 logger.info("Read Furuno hdf5 file")
303 try:
304 ds = xr.open_dataset(file,
305 engine="odim",
306 group=elevation_angle,
307 backend_kwargs=dict(reindex_angle=1.0))
309 except Exception:
310 # except errors if a file could not read
311 logger.warning("Could not load hdf5 file:' " + file + "'.")
312 continue
314 # reindexing azimuth resolution
315 # -----------------------------
316 ds = ds.reindex(azimuth=re_index_parameters,
317 method="nearest")
319 # fix errors by elevation angle
320 # -----------------------------
321 elevation_angle_ds = ds.elevation.median(dim="azimuth")
322 ds["elevation"][:] = elevation_angle_ds
324 # SNR filter by RHOHV
325 # --------------------
327 ds = ds.where(ds.RHOHV > 0.9)
329 # remove static clutter with cmap
330 # --------------------------------
332 dbzh = ds.DBZH.where((cmap == False)) # noqa E712
334 dbzh_attrs = ds.DBZH.attrs
335 ds["DBZH"] = xr.DataArray(dbzh,
336 dims=dims,
337 attrs=dbzh_attrs)
339 # clutter detection
340 # --------------------
341 logger.info("Clutter classification.")
342 ds = wrainfo.clutter.fuzzy_echo_classification(ds,
343 cmap=cmap,
344 dims=dims)
346 # set clutter pixel to NA
347 ds = wrainfo.clutter.dbzh_no_clutter(ds)
349 # phase processing & attenuation correction (ZPHI method)
350 # -------------------------------------------------------
351 logger.info("Phase processing and attenuation correction.")
352 ds = wrainfo.attenuation_corr.attenuation_correction(ds,
353 moment="DBZH_no_clutter",
354 dims=dims)
356 # rain rate retreival
357 # --------------------
358 logger.info("Create precipitation products.")
360 # QPE Z(R)
361 logger.info("Create precipitation product with Z(R) relationship.")
362 ds = wrainfo.precipitation.qpe_zr(ds,
363 "DBZH_CORR",
364 path=path,
365 a=200,
366 b=1.6)
368 # smoothing
369 ds = ds.pad(azimuth=(3, 3),
370 mode="wrap").rolling(azimuth=7,
371 min_periods=3,
372 center=True).mean().isel(azimuth=slice(3, -3))
374 # georeferencing and gridding
375 # ---------------------------
376 logger.info("Georeferencing and gridding")
378 swp_georef = wrainfo.geometry.furuno_georeferencing(ds,
379 path=path)
381 wrainfo.geometry.furuno_sweep_to_netcdf(swp_georef,
382 path=path,
383 data_type=data_type)
384 del ds
385 del swp_georef
387 date += delta
389 return True