Coverage for wrainfo/geometry.py: 99%

1"""Geometry 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 wradlib as wrl 

27import numpy as np 

28import os 

29import xarray as xr 

30from wrainfo.reader import read_config_file 

31 

32 

33# functions to georeferencing and gridding the Furuno data 

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

35 

36def get_target_grid(ds, nb_pixels): 

37 """Get target grid.""" 

38 xgrid = np.linspace(ds.x.min(), ds.x.max(), nb_pixels, dtype=np.float32) 

39 ygrid = np.linspace(ds.y.min(), ds.y.max(), nb_pixels, dtype=np.float32) 

40 grid_xy_raw = np.meshgrid(xgrid, ygrid) 

41 grid_xy_grid = np.dstack((grid_xy_raw[0], grid_xy_raw[1])) 

42 return xgrid, ygrid, grid_xy_grid 

43 

44 

45def get_target_coordinates(grid): 

46 """Get target coordinates.""" 

47 grid_xy = np.stack((grid[..., 0].ravel(), grid[..., 1].ravel()), axis=-1) 

48 return grid_xy 

49 

50 

51def get_source_coordinates(ds): 

52 """Get source coordinates.""" 

53 xy = np.stack((ds.x.values.ravel(), ds.y.values.ravel()), axis=-1) 

54 return xy 

55 

56 

57def create_dataarray(ds, data, moment): 

58 """Create data array.""" 

59 # get moment attrributes 

60 mom_attrs = ds[moment].attrs 

61 # remove _Undetect since it's not in the original data 

62 mom_attrs.pop("_Undetect", None) 

63 # create moment DataArray 

64 mom = xr.DataArray(data.astype(np.float32), 

65 dims=["y", "x"], 

66 attrs=mom_attrs 

67 ) 

68 # fix encoding 

69 mom.encoding = dict(zlib=True, _FillValue=0., least_significant_digit=2) 

70 return mom 

71 

72 

73def create_dataset(ds, moments, xgrid, ygrid): 

74 """Create xarray dataset.""" 

75 # create x,y DataArrays 

76 x = xr.DataArray(xgrid, dims=["x"], attrs=dict(standard_name="projection_x_coordinate", 

77 long_name="x coordinate of projection", 

78 units="m")) 

79 y = xr.DataArray(ygrid, dims=["y"], attrs=dict(standard_name="projection_y_coordinate", 

80 long_name="y coordinate of projection", 

81 units="m")) 

82 data_vars = dict() 

83 data_vars.update(moments) 

84 data_vars.update({"time": ds.reset_coords().time}) 

85 

86 # create Dataset 

87 ds_out = xr.Dataset( 

88 data_vars=data_vars, 

89 coords={"x": x, "y": y}, 

90 ) 

91 # apply CF Conventions 

92 ds_out.attrs['Conventions'] = 'CF-1.5' 

93 return ds_out 

94 

95 

96# georeferencing and gridding of Furuno data 

97# ------------------------------------------ 

98 

99def furuno_georeferencing(ds, path): 

100 """Georeference FURUNO Dataset, Grid/Project moments to raster dataset with EPSG Code. 

101 

102 Parameters 

103 ---------- 

104 ds : xarray.Dataset 

105 moments : list 

106 list of selected data variables from the dataset 

107 path : str 

108 path to configuration file 

109 

110 Returns 

111 ------- 

112 : xarray.Dataset 

113 Georeferenced xarray.Dataset. 

114 """ 

115 # read configuration file for epsg_code and nb_pixels 

116 epsg_code = read_config_file(path=path, selection="epsg_code") 

117 nb_pixels = read_config_file(path=path, selection="nb_pixels") 

118 moments = read_config_file(path=path, selection="moments_in_processed_files") 

119 

120 # create projection 

121 proj = wrl.georef.epsg_to_osr(epsg_code) 

122 

123 # georeference single sweep 

124 ds = ds.pipe(wrl.georef.georeference_dataset, proj=proj) 

125 

126 # get source coordinates 

127 src = get_source_coordinates(ds) 

128 

129 # create target grid 

130 xgrid, ygrid, trg_grid = get_target_grid(ds, nb_pixels) 

131 

132 # get target coordinates 

133 trg = get_target_coordinates(trg_grid) 

134 

135 # setup interpolator 

136 ip_near = wrl.ipol.Nearest(src, trg) 

137 

138 # calculate moments, create DataArrays 

139 moment_dict = dict() 

140 for mom in moments: 

141 res = ip_near(ds[mom].values.ravel()).reshape((len(ygrid), len(xgrid))) 

142 moment_dict[mom] = create_dataarray(ds, res, mom) 

143 

144 # create output dataset 

145 ds_out = create_dataset(ds, moment_dict, xgrid, ygrid) 

146 

147 # write crs (projection) to netcdf_file 

148 ds_out = ds_out.rio.write_crs(epsg_code, grid_mapping_name="transverse_mercator") 

149 ds_out.attrs["elevation_angle"] = ds.elevation.values[0] 

150 

151 # add time dimension 

152 time1 = [1] 

153 ds_out = ds_out.expand_dims(time=time1) 

154 

155 return ds_out 

156 

157 

158# output as netcdf 

159# ---------------- 

160 

161def furuno_sweep_to_netcdf(ds, 

162 path, 

163 data_type=None): 

164 """Create and save georeferenced NetCDF files. 

165 

166 Parameter 

167 --------- 

168 ds : xarray.Dataset 

169 path: str 

170 path to configuration file 

171 data_type : str 

172 level of processing (Level1, Level2, or Level3) 

173 

174 Return 

175 ------ 

176 : NetCDF files 

177 NetCDF files in output directory. 

178 """ 

179 # read configuration file 

180 radar_location_identifier = read_config_file(path=path, selection="radar_location_identifier") 

181 outdir = read_config_file(path=path, selection="output_path_processed_files") 

182 moments = read_config_file(path=path, selection="moments_in_processed_files") 

183 

184 # choose moments, which will save in the netcdf 

185 ds = ds[moments] 

186 el = ds.elevation_angle 

187 

188 # use time, elevation, and epsg code to create filename 

189 t0 = ds.time.values.astype("M8[s]").astype("O")[0] 

190 

191 if data_type is not None: 

192 outfilename = f"{radar_location_identifier}_{t0:%Y%m%d}_{t0:%H%M%S}UTC_elev_{el}_{data_type}.nc" 

193 else: 

194 outfilename = f"{radar_location_identifier}_{t0:%Y%m%d}_{t0:%H%M%S}UTC_elev_{el}.nc" 

195 

196 # extract year, month, day from filename 

197 year = outfilename[10:14] 

198 month = outfilename[14:16] 

199 day = outfilename[16:18] 

200 el = outfilename[29:37] 

201 

202 # create output directory 

203 path1 = outdir + "/" + year + "/" + month + "/" + day + "/" + el + "/" 

204 

205 # absolute output path 

206 outfilename = os.path.join(path1, outfilename) 

207 

208 # output path will create if not exists 

209 if not os.path.exists(path1): 

210 os.makedirs(path1, exist_ok=True) 

211 

212 # check whether output file exists 

213 if not os.path.isfile(outfilename): 

214 # output file 

215 print(f"-- output to {outfilename}") 

216 ds.to_netcdf(outfilename) 

217 

218 return True