"""Functions for working radar instances."""importcopyimportnumpyasnpfromnetCDF4importdate2numfrom..configimportget_fillvaluefrom.importdatetime_utils
[docs]defis_vpt(radar,offset=0.5):""" Determine if a Radar appears to be a vertical pointing scan. This function only verifies that the object is a vertical pointing scan, use the :py:func:`to_vpt` function to convert the radar to a vpt scan if this function returns True. Parameters ---------- radar : Radar Radar object to determine if. offset : float, optional Maximum offset of the elevation from 90 degrees to still consider to be vertically pointing. Returns ------- flag : bool True if the radar appear to be verticle pointing, False if not. """# check that the elevation is within offset of 90 degrees.elev=radar.elevation['data']returnnp.all((elev<90.0+offset)&(elev>90.0-offset))
[docs]defto_vpt(radar,single_scan=True):""" Convert an existing Radar object to represent a vertical pointing scan. This function does not verify that the Radar object contains a vertical pointing scan. To perform such a check use :py:func:`is_vpt`. Parameters ---------- radar : Radar Mislabeled vertical pointing scan Radar object to convert to be properly labeled. This object is converted in place, no copy of the existing data is made. single_scan : bool, optional True to convert the volume to a single scan, any azimuth angle data is lost. False will convert the scan to contain the same number of scans as rays, azimuth angles are retained. """ifsingle_scan:nsweeps=1radar.azimuth['data'][:]=0.0seri=np.array([radar.nrays-1],dtype='int32')radar.sweep_end_ray_index['data']=serielse:nsweeps=radar.nrays# radar.azimuth not adjustedradar.sweep_end_ray_index['data']=np.arange(nsweeps,dtype='int32')radar.scan_type='vpt'radar.nsweeps=nsweepsradar.target_scan_rate=None# no scanningradar.elevation['data'][:]=90.0radar.sweep_number['data']=np.arange(nsweeps,dtype='int32')radar.sweep_mode['data']=np.array(['vertical_pointing']*nsweeps)radar.fixed_angle['data']=np.ones(nsweeps,dtype='float32')*90.0radar.sweep_start_ray_index['data']=np.arange(nsweeps,dtype='int32')ifradar.instrument_parametersisnotNone:forkeyin['prt_mode','follow_mode','polarization_mode']:ifkeyinradar.instrument_parameters:ip_dic=radar.instrument_parameters[key]ip_dic['data']=np.array([ip_dic['data'][0]]*nsweeps)# Attributes that do not need any changes# radar.altitude# radar.altitude_agl# radar.latitude# radar.longitude# radar.range# radar.ngates# radar.nrays# radar.metadata# radar.radar_calibration# radar.time# radar.fields# radar.antenna_transition# radar.scan_ratereturn
[docs]defjoin_radar(radar1,radar2):""" Combine two radar instances into one. Parameters ---------- radar1 : Radar Radar object. radar2 : Radar Radar object. """# must have same gate spacingnew_radar=copy.deepcopy(radar1)new_radar.azimuth['data']=np.append(radar1.azimuth['data'],radar2.azimuth['data'])new_radar.elevation['data']=np.append(radar1.elevation['data'],radar2.elevation['data'])new_radar.fixed_angle['data']=np.append(radar1.fixed_angle['data'],radar2.fixed_angle['data'])new_radar.sweep_number['data']=np.append(radar1.sweep_number['data'],radar2.sweep_number['data'])new_radar.sweep_start_ray_index['data']=np.append(radar1.sweep_start_ray_index['data'],radar2.sweep_start_ray_index['data']+radar1.nrays)new_radar.sweep_end_ray_index['data']=np.append(radar1.sweep_end_ray_index['data'],radar2.sweep_end_ray_index['data']+radar1.nrays)new_radar.nsweeps+=radar2.nsweepsnew_radar.sweep_mode['data']=np.append(radar1.sweep_mode['data'],radar2.sweep_mode['data'])iflen(radar1.range['data'])>=len(radar2.range['data']):new_radar.range['data']=radar1.range['data']else:new_radar.range['data']=radar2.range['data']new_radar.ngates=len(new_radar.range['data'])# to combine times we need to reference them to a standard# for this we'll use epoch timer1num=datetime_utils.datetimes_from_radar(radar1,epoch=True)r2num=datetime_utils.datetimes_from_radar(radar2,epoch=True)new_radar.time['data']=date2num(np.append(r1num,r2num),datetime_utils.EPOCH_UNITS)new_radar.time['units']=datetime_utils.EPOCH_UNITSnew_radar.nrays=len(new_radar.time['data'])forvarinnew_radar.fields.keys():sh1=radar1.fields[var]['data'].shapesh2=radar2.fields[var]['data'].shapenew_field_shape=(sh1[0]+sh2[0],max(sh1[1],sh2[1]))new_field=np.ma.zeros(new_field_shape)new_field[:]=np.ma.maskednew_field.set_fill_value(get_fillvalue())new_field[0:sh1[0],0:sh1[1]]=radar1.fields[var]['data']new_field[sh1[0]:,0:sh2[1]]=radar2.fields[var]['data']new_radar.fields[var]['data']=new_field# radar locations# TODO moving platforms - any more?if(len(radar1.latitude['data'])==1&len(radar2.latitude['data'])==1&len(radar1.longitude['data'])==1&len(radar2.longitude['data'])==1&len(radar1.altitude['data'])==1&len(radar2.altitude['data'])==1):lat1=float(radar1.latitude['data'])lon1=float(radar1.longitude['data'])alt1=float(radar1.altitude['data'])lat2=float(radar2.latitude['data'])lon2=float(radar2.longitude['data'])alt2=float(radar2.altitude['data'])if(lat1!=lat2)or(lon1!=lon2)or(alt1!=alt2):ones1=np.ones(len(radar1.time['data']),dtype='float32')ones2=np.ones(len(radar2.time['data']),dtype='float32')new_radar.latitude['data']=np.append(ones1*lat1,ones2*lat2)new_radar.longitude['data']=np.append(ones1*lon1,ones2*lon2)new_radar.latitude['data']=np.append(ones1*alt1,ones2*alt2)else:new_radar.latitude['data']=radar1.latitude['data']new_radar.longitude['data']=radar1.longitude['data']new_radar.altitude['data']=radar1.altitude['data']else:new_radar.latitude['data']=np.append(radar1.latitude['data'],radar2.latitude['data'])new_radar.longitude['data']=np.append(radar1.longitude['data'],radar2.longitude['data'])new_radar.altitude['data']=np.append(radar1.altitude['data'],radar2.altitude['data'])returnnew_radar
[docs]defimage_mute_radar(radar,field,mute_field,mute_threshold,field_threshold=None):""" This function will split a field based on thresholds from another field. Specifically, it was designed to separate areas of reflectivity where the correlation coefficient is less than a certain threshold to discern melting precipitation. Parameters ---------- radar : Radar Radar instance which provides the fields for muting. field : str Name of field to image mute. mute_field : str Name of field to image mute by. mute_threshold : float Threshold value to mute by. field_threshold : float Additional threshold to mask. Returns ------- radar : Radar Radar object with 2 new fields from input field, one muted and one not muted. """# add checks for field availabilityiffieldnotinradar.fields.keys():raiseKeyError('Failed - ',field,' field to mute not found in Radar object.')ifmute_fieldnotinradar.fields.keys():raiseKeyError('Failed - ',mute_field,' field to mute by not found in Radar object.')# get data from fieldsdata_to_mute=radar.fields[field]['data']data_mute_by=radar.fields[mute_field]['data']# create filters# field_filter is used if user wants to use additional criteria in the# original fieldiffield_thresholdisnotNone:field_filter=data_to_mute>=field_thresholdelse:field_filter=None# mute_filter will be the primary filter for determining muted regionsmute_filter=data_mute_by<=mute_threshold# mute_mask is the combined filteriffield_filterisNone:mute_mask=mute_filterelse:mute_mask=mute_filter&field_filter# break up the field into muted regions and non muted regionsnon_muted_field=np.ma.masked_where(mute_mask,data_to_mute)non_muted_field=np.ma.masked_invalid(non_muted_field)muted_field=np.ma.masked_where(~mute_mask,data_to_mute)muted_field=np.ma.masked_invalid(muted_field)# add fields to a dictionary and save to radar objectnon_muted_dict=radar.fields[field].copy()non_muted_dict['data']=non_muted_fieldnon_muted_dict['long_name']='Non-muted '+fieldradar.add_field('nonmuted_'+field,non_muted_dict)muted_dict=radar.fields[field].copy()muted_dict['data']=muted_fieldmuted_dict['long_name']='Muted '+fieldradar.add_field('muted_'+field,muted_dict)returnradar