RadarData¶

This page contains the documentation for the RadarData class, which is the basic object in ImpDAR. If you are interacting with the API in a significant way, this is where you will find documentation from most of the things you care about, particularly how the data is stored and how to do basic processing steps on it. All of the files to define the class are in impdar/lib/Radardata, with the basic initialization and class properties found in __init__.py and addional functionality spread across _RadarDataSaving, _RadarDataFiltering, and _RadarDataProcessing.

RadarData Base¶

class impdar.lib.RadarData.RadarData(fn_mat)¶

A class that holds the relevant information for a radar profile.

We keep track of processing steps with the flags attribute. This base version’s __init__ takes a filename of a .mat file in the old StODeep format to load.

attrs_guaranteed = ['chan', 'data', 'decday', 'dt', 'pressure', 'snum', 'tnum', 'trace_int', 'trace_num', 'travel_time', 'trig', 'trig_level']¶: Attributes that every RadarData object should have. These should not be None.

attrs_optional = ['nmo_depth', 'lat', 'long', 'elev', 'dist', 'x_coord', 'y_coord', 'fn', 't_srs']¶: Optional attributes that may be None without affecting processing. These may not have existed in old StoDeep files, and they often cannot be set at the initial data load. If they exist, they all have units of meters.

chan¶: The Channel number of the data.

check_attrs()¶

Check if required attributes exist.

This is largely for development only; loaders should generally call this method last, so that they can confirm that they have defined the necessary attributes.

Raises: ImpdarError – If any required attribute is None, or any optional attribute is fully absent

data¶: np.ndarray(snum x tnum) of the actual return power.

decday¶: np.ndarray(tnum,) of the acquisition time of each trace note that this is referenced to Jan 1, 0 CE (matlabe datenum) for convenience, use the datetime attribute to access a python version of the day

dist¶: np.ndarray(tnum,) of the distances along the profile. units will depend on whether geographic coordinate transforms, as well as GPS data, are available.

dt¶: float, The spacing between samples in travel time, in seconds.

elev¶: np.ndarray(tnum,) Optional. Elevation along the profile.

fn¶: str, the input filename. May be left as None.

lat¶: np.ndarray(tnum,) latitude along the profile. Generally not in projected coordinates

long¶: np.ndarray(tnum,) longitude along the profile. Generally not in projected coords.

nmo_depth¶: np.ndarray(tnum,) Optional. Depth of each trace below the surface

pressure¶: np.ndarray(tnum,) The pressure at acquisition. ImpDAR does not use this at present.

snum¶: int number of samples per trace

tnum¶: int, the number of traces in the file

trace_int¶: float, the time between traces.

trace_num¶: np.ndarray(tnum,) The 1-indexed number of the trace

travel_time¶: np.ndarray(snum,) The two way travel time to each sample, in us

trig¶: np.ndarray(tnum,) the index in each trace where the triggered.

trig_level¶: float, The value on which the radar was triggering.

x_coord¶: np.ndarray(tnum,) Optional. Projected x-coordinate along the profile.

y_coord¶: np.ndarray(tnum,) Optional. Projected y-coordinate along the profile.

Saving RadarData¶

These are all instance methods for saving information from a RadarData object. They are defined in impdar/lib/RadarData/_RadarDataSaving.py.

RadarData.save(fn)¶

Save the radar data.

Parameters: fn (str) – Filename. Should have a .mat extension

RadarData.save_as_segy(fn)¶

Save as a (non standard-compliant) SEGY file that can be used with, e.g., SeisUNIX.

Parameters: fn (str) – The filename to save as.
Raises: ImportError – If segyio cannot be imported.

RadarData.output_shp(fn, t_srs=None, target_out=None)¶

Output a shapefile of the traces.

If there are any picks, we want to output these. If not, we will only output the tracenumber. This function requires osr/gdal for shapefile creation. I suggest exporting a csv if you don’t want to deal with gdal.

Parameters

fn (str) – The filename of the output
t_srs (int, optional) – EPSG number of the target spatial reference system. Default 4326 (wgs84)
target_out (str, optional) – Used to overwrite the default output format of picks. By default, try to write depth and if there is no nmo_depth use TWTT. You might want to use this to get the output in TWTT or sample number (options are depth, elev, twtt, snum)

Raises

ImportError – If osgeo cannot be imported

RadarData.output_csv(fn, target_out=None, delimiter=',')¶

Output a csv of the traces.

If there are any picks, we want to output these. If not, we will only output the tracenumber.

Parameters

fn (str) – The filename of the output
target_out (str, optional) – Used to overwrite the default output format of picks. By default, try to write depth and if there is no nmo_depth use TWTT. You might want to use this to get the output in TWTT or sample number (options are depth, elev, twtt, snum)
delimiter (str, optional) – Delimiter for csv. Default ‘,’.

Processing RadarData¶

These are all instance methods for processing data on a RadarData object. They are defined in impdar/lib/RadarData/_RadarDataProcessing.py.

RadarData.reverse()¶

Reverse radar data

Essentially flip the profile left-right. This is desirable in a number of instances, but is particularly useful for concatenating profiles acquired in opposite directions. The St Olaf version of this function had a bunch of options. They seemed irrelevant to me. We just try to flip everything that might need flipping.

RadarData.nmo(ant_sep, uice=169000000.0, uair=300000000.0, const_firn_offset=None, rho_profile=None, permittivity_model=<function firn_permittivity>, const_sample=True)¶

Normal move-out correction.

Converts travel time to distance accounting for antenna separation. This potentially affects data and snum. It also defines nmo_depth, the moveout-corrected depth

Updated to accomodate vertical velocity profile. B Hills 8/2019

Parameters

ant_sep (float) – Antenna separation in meters.
uice (float or np.ndarray (2 x m), optional) – Speed of light in ice, in m/s. (different from StoDeep!!). Default 1.69e8.
uair (float, optional) – Speed of light in air. Default 3.0e8
const_firn_offset (float, optional) – Offset all depths by this much to account for firn-air. THIS IS THE TWO-WAY THICKNESS (not the firn-air). Useful if data start below thesurface so you can skip complicated, depth-dependent corrections. Default None.
rho_profile (str,optional) – Filename for a csv file with density profile (depths in first column and densities in second) Units should be meters for depth, kgs per meter cubed for density. Note that using a variable uice will break the linear scaling between depth and time, so we are forced to choose whether the y-axis is linear in speed or time. I chose time, since this eases calculations like migration. For plotting vs. depth, however, the functions just use the bounds, so the depth variations are averaged out. You can use the helper function constant_sample_depth_spacing() in order to correct this, but you should call that after migration.
permittivity_model (fun, optional) – density to permittivity model from the literature
const_sample (bool, optional) – interpolate to constant sample spacing after the nmo correction

RadarData.crop(lim, top_or_bottom='top', dimension='snum', uice=169000000.0, rezero=True, zero_trig=True)¶

Crop the radar data in the vertical. We can take off the top or bottom.

This will affect data, travel_time, and snum.

Parameters

lim (float (int if dimension=='snum')) – The value at which to crop.
top_or_bottom (str, optional) – Crop off the top (lim is the first remaining return) or the bottom (lim is the last remaining return).
dimension (str, optional) – Evaluate in terms of sample (snum), travel time (twtt), or depth (depth). If depth, uses nmo_depth if present and use uice with no transmit/receive separation. If pretrig, uses the recorded trigger sample to crop.
rezero (bool, optional) – Set the zero on the y axis to the cropped value (if cropping off the top). Default True. This is desirable if you are zeroing to the surface.
zero_trig (bool, optional) – Reset the trigger to zero. Effectively asserts that the crop was to the surface. Default True.
uice (float, optional) – Speed of light in ice. Used if nmo_depth is None and dimension==’depth’

RadarData.restack(traces)¶

Restack all relevant data to the given number of traces.

This function just takes the average of the given number of traces. This reduces file size and can get rid of noise. There are fancier ways to do this— if you have GPS, you probably want to restack to constant trace spacing instead.

Parameters: traces (int) – The (odd) number of traces to stack

RadarData.rangegain(slope)¶

Apply a range gain.

Parameters: slope (float) – The slope of the linear range gain to be applied. Maybe try 1.0e-2?

RadarData.agc(window=50, scaling_factor=50)¶

Try to do some automatic gain control

This is from StoDeep–I’m not sure it is useful but it was easy to roll over so I’m going to keep it. I think you should have most of this gone with a bandpass, but whatever.

Parameters

window (int, optional) – The size of window we use in number of samples (default 50)
scaling_factor (int, optional) – The scaling factor. This gets divided by the max amplitude when we rescale the input. Default 50.

RadarData.constant_space(spacing, min_movement=0.01, show_nomove=False)¶

Restack the radar data to a constant spacing.

This method uses the GPS information (i.e. the distance, x, y, lat, and lon), to do a 1-d interpolation to get new values in the radargram. It also updates related variables like lat, long, elevation, and coordinates. To avoid retaining sections of the radargram when the antenna was in fact stationary, some minimum movement between traces is enforced. This value is in meters, and should change to be commensurate with the collection strategy (e.g. skiing a radar is slower than towing it with a snowmobile).

This function comprises the second half of what was done by StoDeep’s interpdeep. If you have GPS data from an external, high-precision GPS, you would first want to call impdar.lib.gpslib.kinematic_gps_control so that the GPS-related variables are all improved, then you would want to call this method. impdar.lib.gpslib provides some wrappings for doing both steps and for loading in the external GPS data.

Parameters

spacing (float) – Target trace spacing, in meters
min_movement (float, optional) – Minimum trace spacing. If there is not this much separation, toss the next shot. Set high to keep everything. Default 1.0e-2.
show_nomove (bool, optional) – If True, make a plot shading the areas where we think there is no movement. This can be really helpful for diagnosing what is wrong if you have lingering stationary traces. Untested.

RadarData.elev_correct(v_avg=169000000.0)¶

Move the surface down in the data array to account for surface elevation.

NMO depth attribute must have been created before elev_correct is called. This method should generally be called after you have interpolated precision GPS onto the data, otherwise the noise a handheld GPS will make the results look pretty bad.

Be aware that this is not a precise conversion if your nmo correction had antenna separation, or if you used a depth-variable velocity structure. This is because we have a single vector describing depths in the data array, so only one value for each depth step.

Parameters: v_avg (float, optional) – Average velocity. This is what will define the depth slices in the new data array. Default 1.69e8.
Raises: ValueError: – If there is no nmo_depth since this is used for calculating depths

Filtering Radar Data¶

These are all instance methods for filtering data to remove noise. They are defined in impdar/lib/RadarData/_RadarDataFiltering.py.

RadarData.migrate(mtype='stolt', vtaper=10, htaper=10, tmig=0, vel_fn=None, vel=168000000.0, nxpad=10, nearfield=False, verbose=0)¶

Migrate the data.

This is a wrapper around all the migration routines in migration_routines.py.

Parameters: mtype (str, optional) – The chosen migration routine. Options are: kirch, stolt, phsh. Default: stolt

RadarData.vertical_band_pass(low, high, order=5, filttype='butter', cheb_rp=5, fir_window='hamming', *args, **kwargs)¶

Vertically bandpass the data

This function uses a forward-backward filter on the data. Returns power that is not near the wavelength of the transmitter is assumed to be noise, so the limits for the filtering should generally surround the radar frequency. Some experimentation may be needed to see what gives the clearest results for any given data

There are a number of options for the filter type. Depending on the type of filter chosen, some other p

Parameters

low (float) – Lowest frequency passed, in MHz
high (float) – Highest frequency passed, in MHz
order (int) – Filter order (default 5)
filttype (str, optional) – The filter type to use. Options are butter(worth), cheb(yshev type I), bessel, or FIR (finite impulse response). Default is butter.
cheb_rp (float, optional) – Maximum ripple, in decibels, of Chebyshev filter. Only used if filttype==’cheb’. Default is 5.
fir_window (str, optional) – The window type passed to scipy.signal.firwin. Only used if filttype==’fir’. Default is hamming’

RadarData.adaptivehfilt(window_size, *args, **kwargs)¶

Adaptively filter to reduce noise in upper layers

This subtracts the average of traces around an individual trace in order to filter it. You can call this method directly, or it can be called by sending the ‘adaptive’ option to RadarData.hfilt()

Parameters

window_size (int) – number of traces to include in the moving average to be removed
Documentation (Original StoDeep) –

HFILTDEEP-This StoDeep subroutine processes bandpass filtered
or NMO data to reduce the horizontal noise in the upper layers. The user need not specify any frequencies. This program simply takes the average of all of the traces and subtracts it from the bandpassed data. It will remove most horizontally-oriented features in a radar profile: ringing, horizontal reflectors. It will also create artifacts at travel-times corresponding to any bright horizontal reflectors included in the average trace.

You will want to experiment with the creation of the average trace. Ideally choose an area where all reflectors are sloped and relatively dim so that they average out while the horizontal noise is amplified. Note that generally there is no perfect horizontal filter that will work at all depths. You will have to experiment to get the best results for your area of interest.

WARNING: Do not use hfiltdeep on elevation-corrected data!!!

Created: Logan Smith - 6/12/02 Modified by: L. Smith, 5/27/03. K. Dwyer, 6/3/03. B. Welch, 5/2/06. B. Youngblood, 6/13/08. J. Olson, 7/10/08

RadarData.horizontalfilt(ntr1, ntr2, *args, **kwargs)¶

Remove the average trace.

Parameters

ntr1 (int) – Leftmost trace for averaging
ntr2 (int) – Rightmost trace for averaging

Original StoDeep Documentation:

HFILTDEEP - This StoDeep subroutine processes bandpass filtered or NMO data to reduce the horizontal noise in the upper layers. The user need not specify any frequencies. This program simply takes the average of all of the traces and subtracts it from the bandpassed data. It will remove most horizontally-oriented features in a radar profile: ringing, horizontal reflectors. It will also create artifacts at travel-times corresponding to any bright horizontal reflectors included in the average trace.

You will want to experiment with the creation of the average trace. Ideally choose an area where all reflectors are sloped and relatively dim so that they average out while the horizontal noise is amplified. Note that generally there is no perfect horizontal filter that will work at all depths. You will have to experiment to get the best results for your area of interest.

RadarData.highpass(wavelength)¶

High pass in the horizontal for a given wavelength.

This only works if the data have constant trace spacing; we check the processing flags to enforce this.

Parameters: wavelength (int) – The wavelength to pass, in meters.

Original StoDeep Documentation:

HIGHPASSDEEP - This is NOT a highpass frequency filter–rather it is a horizontal filter to be used after interpolation because our data now has constant spacing and a constant time. Note that this horizontal filter requires constant trace-spacing in order to be effective.

You will want to experiment with the creation of the average trace. Ideally choose an area where all reflectors are sloped and relatively dim so that they average out while the horizontal noise is amplified. Note that generally there is no perfect horizontal filter that will work at all depths. You will have to experiment to get the best results for your area of interest.

WARNING: Do not use highpassdeep on elevation-corrected data!!!

Created by L. Smith and modified by A. Hagen, 6/15/04. B. Welch, 5/3/06. J. Werner, 6/30/08. J. Olson, 7/10/08

RadarData.winavg_hfilt(avg_win, taper='full', filtdepth=100)¶

Uses a moving window to find the average trace, then subtracts this from the data.

Parameters

avg_win (int) – The size of moving window. Must be odd and less than tnum. We will correct these rather than raise an exception.
taper (str) – How the filter varies with depth. Options are full or pexp. For full, the power tapers exponentially. For pexp, the filter stops after the sample number given by filtdepth.

Original StoDeep Documentation:

WINAVG_HFILTDEEP - This StoDeep subroutine performs a horizontal filter on the data to reduce the ringing in the upper layers. It uses a moving window to create an average trace for each individual trace, applies an exponential taper to it and then subtracts it from the trace. This is based off of the original hfiltdeep and uses the moving window designed in cosine_win_hfiltdeep. The extent of the taper can help to minimize artifacts that are often produced near regions of bright bed reflectors.

You will want to experiment with the creation of the average trace. Ideally choose an area where all reflectors are sloped and relatively dim so that they average out while the horizontal noise is amplified. Note that generally there is no perfect horizontal filter that will work at all depths. You will have to experiment to get the best results for your area of interest.

WARNING: Do not use winavg_hfiltdeep on elevation-corrected data!!!

Created: Kieran Dwyer 6/18/03 Modified by B. Welch, 5/2/06. J. Olson, 7/10/08.

RadarData.hfilt(ftype='hfilt', bounds=None, window_size=None)¶

Horizontally filter the data.

This is a wrapper around other filter types. Horizontal filters are implemented (and documented) in the impdar.lib.horizontal_filters module.

Parameters

ftype (str, optional) – The filter type. Options are hfilt and adaptive. Default hfilt
bounds (tuple, optional) – Bounds for the hfilt. Default is None, but required if ftype is hfilt.
window_size (int, optional) – number of traces in the moving average. Default is None, but required if ftype is adaptive.