"""pyrad.proc.process_dem======================Functions to manage DEM data.. autosummary:: :toctree: generated/ process_dem process_visibility process_gecsx"""fromcopyimportdeepcopyfromwarningsimportwarnimportnumpyasnpimportpyartfrom..io.io_auximportget_datatype_fields,get_fieldname_pyartfrom..io.read_data_demimportread_dem,dem2radar_data# from memory_profiler import profile
[docs]defprocess_dem(procstatus,dscfg,radar_list=None):""" Gets DEM data and put it in radar coordinates Parameters ---------- procstatus : int Processing status: 0 initializing, 1 processing volume, 2 post-processing dscfg : dictionary of dictionaries data set configuration. Accepted Configuration Keywords:: datatype : string. Dataset keyword arbitrary data type supported by pyrad keep_in_memory : int. Dataset keyword if set keeps the COSMO data dict, the COSMO coordinates dict and the COSMO field in radar coordinates in memory. Default False regular_grid : int. Dataset keyword if set it is assume that the radar has a grid constant in time and there is no need to compute a new COSMO field if the COSMO data has not changed. Default False dem_field : str. Dataset keyword name of the DEM field to process demfile : str. Dataset keyword Name of the file containing the DEM data radar_list : list of Radar objects Optional. list of radar objects Returns ------- new_dataset : dict dictionary containing the output field with name corresponding to dem_field ind_rad : int radar index """ifprocstatus!=1:returnNone,None# debugging# start_time = time.time()fordatatypedescrindscfg["datatype"]:radarnr,_,_,_,_=get_datatype_fields(datatypedescr)breakind_rad=int(radarnr[5:8])-1ifradar_list[ind_rad]isNone:warn("No valid radar")returnNone,Noneradar=radar_list[ind_rad]keep_in_memory=dscfg.get("keep_in_memory",0)regular_grid=dscfg.get("regular_grid",0)field_name=get_fieldname_pyart(dscfg["dem_field"])fname=dscfg["dempath"][ind_rad]+dscfg["demfile"]ifkeep_in_memory:ifdscfg["initialized"]==0:dem_data=read_dem(fname,field_name)dscfg["global_data"]={"dem_data":dem_data,"dem_field":None}ifregular_grid:dscfg["global_data"]["dem_field"]=dem2radar_data(radar,dem_data,field_name)dscfg["initialized"]=1dem_data=dscfg["global_data"]["dem_data"]else:dem_data=read_dem(fname,field_name)ifdem_dataisNone:warn("DEM data not found")returnNone,Noneifregular_grid:print("DEM field already in memory")dem_field=dscfg["global_data"]["dem_fields"]else:dem_field=dem2radar_data(radar,dem_data,field_name=field_name)ifdem_fieldisNone:warn("Unable to obtain DEM fields")returnNone,None# prepare for exitnew_dataset={"radar_out":deepcopy(radar)}new_dataset["radar_out"].fields=dict()new_dataset["radar_out"].add_field(field_name,dem_field)returnnew_dataset,ind_rad
[docs]defprocess_visibility(procstatus,dscfg,radar_list=None):""" Gets the visibility in percentage from the minimum visible elevation. Anything with elevation lower than the minimum visible elevation plus and offset is set to 0 while above is set to 100. Parameters ---------- procstatus : int Processing status: 0 initializing, 1 processing volume, 2 post-processing dscfg : dictionary of dictionaries data set configuration. Accepted Configuration Keywords:: datatype : string. Dataset keyword arbitrary data type supported by pyrad offset : float. Dataset keyword The offset above the minimum visibility that must be filtered radar_list : list of Radar objects Optional. list of radar objects Returns ------- new_dataset : dict dictionary containing the output field "visibility" ind_rad : int radar index """ifprocstatus!=1:returnNone,Nonefordatatypedescrindscfg["datatype"]:radarnr,_,datatype,_,_=get_datatype_fields(datatypedescr)ifdatatype=="minvisel":minvisel_field=get_fieldname_pyart(datatype)breakind_rad=int(radarnr[5:8])-1ifradar_list[ind_rad]isNone:warn("No valid radar")returnNone,Noneradar=radar_list[ind_rad]offset=dscfg.get("offset",0.0)minvisel_data=radar.fields[minvisel_field]["data"]+offsetele_data=np.broadcast_to(radar.elevation["data"].reshape(radar.nrays,1),(radar.nrays,radar.ngates))vis_dict=pyart.config.get_metadata("visibility")vis_dict["data"]=100.0*np.ma.greater_equal(ele_data,minvisel_data,dtype=float)# if a gate has visibility 0 all the subsequent gates in the ray# are set to 0forrayinrange(radar.nrays):ind=np.where(vis_dict["data"][ray,:]==0.0)[0]ifind.size>0:vis_dict["data"][ray,ind[0]:]=0.0# prepare for exitnew_dataset={"radar_out":deepcopy(radar)}new_dataset["radar_out"].fields=dict()new_dataset["radar_out"].add_field("visibility",vis_dict)returnnew_dataset,ind_rad
[docs]defprocess_gecsx(procstatus,dscfg,radar_list=None):""" Computes ground clutter RCS, radar visibility and many others using the GECSX algorithmn translated from IDL into python Parameters ---------- procstatus : int Processing status: 0 initializing, 1 processing volume, 2 post-processing dscfg : dictionary of dictionaries data set configuration. Accepted Configuration Keywords:: datatype : list of string. Dataset keyword arbitrary data type supported by pyrad range_discretization : float. Dataset keyword Range discretization used when computing the Cartesian visibility field the larger the better but the slower the processing will be az_discretization : float. Dataset keyword Azimuth discretization used when computing the Cartesian visibility field, the larger the better but the slower the processing will be ke : float. Dataset keyword Equivalent earth-radius factor used in the computation of the radar beam refraction atm_att : float. Dataset keyword One-way atmospheric refraction in db / km mosotti_kw : float. Dataset keyword Clausius-Mosotti factor K, depends on material (water) and wavelength for water = sqrt(0.93) raster_oversampling : int. Dataset keyword The raster resolution of the DEM should be smaller than the range resolution of the radar (defined by the pulse length). If this is not the case, this keyword can be set to increase the raster resolution. The values for the elevation, sigma naught, visibility are repeated. The other values are recalculated. Values for raster_oversampling: 0 or undefined: No oversampling is done 1: Oversampling is done. The factor N is automatically calculated such that 2*dx/N < pulse length 2 or larger: Oversampling is done with this value as N sigma0_method : string. Dataset keyword Which estimation method to use, either 'Gabella' or 'Delrieu' clip : int. Dataset keyword If set to true, the provided DEM will be clipped to the extent of the polar radar domain. Increases computation speed a lot but Cartesian output fields will be available only over radar domain radar_list : list of Radar objects Optional. list of radar objects Returns ------- new_dataset : list of dict list of dictionaries containing the polar data output and the Cartesian data output in this order The first dictionary (polar) contains the following fields: "rcs_clutter", "dBm_clutter", "dBZ_clutter" and "visibility_polar" The second dictionary (cart) contains the following fields: "bent_terrain_altitude", "terrain_slope", "terrain_aspect", "elevation_angle", "min_vis_elevation", "min_vis_altitude", "incident_angle", "sigma_0", "effective_area" ind_rad : int radar index """ifprocstatus!=1:returnNone,Nonefordatatypedescrindscfg["datatype"]:radarnr,_,_,_,_=get_datatype_fields(datatypedescr)ind_rad=int(radarnr[5:8])-1fname=dscfg["dempath"][ind_rad]+dscfg["demfile"]demproj=Noneif"demproj"indscfg.keys():demproj=dscfg["demproj"]try:demproj=int(demproj)exceptValueError:# demproj is not an EPSG intpassdem_data=read_dem(fname,projparams=demproj)# If no radar data is provided we create empty radar object from user# specificationiflen(radar_list)==0:ranges=np.arange(dscfg["range_resolution"]/2,dscfg["rmax"],dscfg["range_resolution"])azimuths=np.arange(dscfg["azmin"],dscfg["azmax"],dscfg["anglestep"])elevations=dscfg["antenna_elevations"]radar=pyart.testing.make_empty_ppi_radar(len(ranges),len(azimuths),len(elevations))radar.latitude["data"]=np.array(dscfg["RadarPosition"]["latitude"])radar.longitude["data"]=np.array(dscfg["RadarPosition"]["longitude"])radar.altitude["data"]=np.array(dscfg["RadarPosition"]["altitude"])radar.azimuth["data"]=np.array(list(azimuths)*len(elevations))radar.range["data"]=rangesradar.fixed_angle["data"]=np.array(elevations)radar.elevation["data"]=np.array([len(azimuths)*[e]foreinelevations]).ravel()# change radar nameradar.metadata["instrument_name"]=dscfg["RadarName"]else:radar=radar_list[0]if"antenna_elevations"indscfg:# Intersection between radar elevation angles and config choice# using a certain numerical tolerance as radar angles are# sometimes coded as 0.699996 for 0.7 degrees in the radar files# for exampleel1=radar.fixed_angle["data"].astype(float)el2=dscfg["antenna_elevations"].astype(float)idx_to_process=[iforiinrange(len(el1))ifnp.any(np.isclose(el1[i],el2))]print("Radar elevations angles redefined in config file")print("Elevation angles {:s} will be processed".format(str([el1[i]foriinidx_to_process])))radar=radar.extract_sweeps(idx_to_process)# Create dict with radar specificationsradar_specs={}radar_specs["frequency"]=dscfg["frequency"][ind_rad]radar_specs["loss"]=dscfg["lrxh"][ind_rad]+dscfg["mflossh"][ind_rad]radar_specs["power"]=dscfg["txpwrh"][ind_rad]radar_specs["tau"]=dscfg["pulse_width"][ind_rad]radar_specs["beamwidth"]=dscfg["radar_beam_width_h"][ind_rad]radar_specs["gain"]=dscfg["AntennaGainH"][ind_rad]az_conv=dscfg.get("AzimTol",0)[ind_rad]ke=dscfg.get("refcorr",4/3.0)[ind_rad]atm_att=dscfg.get("attg",0.012)[ind_rad]mosotti_kw=dscfg.get("mosotti_factor",0.9644)[0]sigma0_method=dscfg.get("sigma0_method","Gabella")raster_oversampling=dscfg.get("raster_oversampling",1)verbose=dscfg.get("verbose",1)clip=dscfg.get("clip",1)daz=dscfg.get("az_discretization",0.2)dr=dscfg.get("range_discretization",100)gecsx_grid,gecsx_radar=pyart.retrieve.gecsx(radar,radar_specs,dem_data,fill_value=None,az_conv=az_conv,dr=dr,daz=daz,ke=ke,atm_att=atm_att,mosotti_kw=mosotti_kw,raster_oversampling=raster_oversampling,sigma0_method=sigma0_method,clip=clip,return_pyart_objects=True,verbose=verbose,)new_dataset=[{"radar_out":gecsx_grid},{"radar_out":gecsx_radar}]returnnew_dataset,ind_rad