Igor Code

The objective of this page is to serve a repository of useful Igor code that we keep needing for different applications. Feel free to add new code or to comment on existing code, but please don't delete old code w/o discussing with Jose.

There is a a "General Macros" igor procedure file that has a large suite of various Igor functions. It can be downloaded here. There is an Igor Help file (with clickable links) that helps to describe the various functions in this ipf that can be downloaded here. Some background: this ipf was created by Donna as a good-bye 'gift' to my NOAA friends. Two new versions of Igor have since been released, and several tools are now dated and/or obsolete. Still, I think it is a good reference for general programming techniques, and as a resource for a variety of tasks. - Donna

JGGenTools

At a meeting on July 17, 2014 with Jose, Doug, Pedro, Donna, Harald, a new scheme for organizing code was instituted and debuted to the group on July 18, 2014. There is a google document with links and function names here.

Multiplexed Experiment Tool (MPT) Code

This was also formerly known as the PAM tool, OFR tool and briefly, MET. This info can now be found at here.

Both ipfs are designed to work within a squirrel or pika experiment. For users who wish to use the code alone. If you wish to use the code not within a squirrel or pika experiments, some extra functions are needed to get the OFR code to compile. You can download a separate ipf with these functions here.

Doing a Robust Linear Regression

Linear regression gets pretty complex when one starts thinking about the details. There are 3 important decisions, which depend on the data which you are trying to fit:

• (A) Error in one or both coordinates:
  • (A1) Only one of the coordinates (Y) has error and the other does not (e.g. when fitting an aerosol measurement vs. time).
  • (A2) Or do both coordinates have errors of the same relative order-of-magnitude? (e.g. as when comparing two instruments that measure the same quantity)
• (B) Assumptions about the relative magnitudes of the errors in different points
  • (B1) All the points are assumed to have the same error
  • (B2) The error varies point-by-point.
• (C) Standard vs. Robust regression
  • (C1) The fit line is found by minimizing Chi^2, i.e. the sum of the squares of the deviations between each point and the regression line.
  • (C2) The fit line is found by minimizing the sum of the absolute values of the deviations between each point and the regression line.

This gives rise to 2 x 2 x 2 = 8 options to choose from when fitting your data. Which to use depends on the nature of the data and the use of the fit, and should be chosen intelligently. In practice if there is a good relationship and the points are scattered around a line with few outliers, all the methods will produce similar results. However when there are outliers or significant scatter, there can be significant differences and some of the methods are plainly incorrect for some datasets.

To do regressions in Igor in each case:

• A1-B1-C1: this is the default assumption of the linear fit
• A2-B1-C1: an approximation to this fit can be done with the /ODR=2 option in Igor.
• A1-B1-C2: this can be done in Igor, but you need to define a new custom fit function (which can be just a line) and then set V_FitOptions = 2 in the command line. (Going back to V_FitOptions = 0 gets you back to the regular A1-B1-C1 fits)
• A2-B1-C2: this can be accomplished approximately by doing A1-B1-C2 for Y vs X, then for X vs Y, and averaging the fitting parameters appropriately (and the fitting parameter uncertainties in quadrature). There are some papers that show that this is approximately correct.
• For all the B2 options, choose the fitting weights in the Igor Curve Fitting tab

Regridding a Wave into the Time Axis of Another Wave

To regrid a wave into another framework, use the following functions. For example you might have one curve, defined by a certain x- and y-wave (xwave1, ywave1). Let's say there are 10 points in these waves, and the maximum is 20 (and extreme case of length discrepancy). Another curve is defined by two other x- and y-waves, also with a maximum of 20, but with 20 points each (xwave2, ywave2). It is important that they have approximately the same bound values, or it will have a hard time interpolating between them.

If I want to look at the first curve in terms of the second x-wave, you have to regrid the first y-wave into the x-wave of the second. This can be done in two ways, and that is dependent on which inintial set of waves has more data points. If the wave you want to regrid into has more points than the original, then use the first of these functions (Interp_1). If the opposite is true (wave you want to regrid to has less points) use the second function (Interp_2). They can be run as is, after changing the wave names appropriately. The [p] refers to each point, and should stay with this syntax (replacing it with a number refers only to that specified cell).

Function Interp_1 ()

Wave xwave1, ywave1, xwave2 //ywave2 is not used

duplicate/o xwave2 ywave_regrid

ywave_regrid = ywave1 [binarysearchinterp(xwave1, xwave2[p])]

appendtotable ywave_regrid

End Function

Function Interp_2 ()

Wave xwave1, ywave2, xwave2

duplicate/o xwave1 ywave_regrid

ywave_regrid = ywave1 [fAverageXY(xwave2, xwave1, xwave2[p-1], xwave2[p])]

End Function

Compiled by Qi Zhang and Alex Huffman, 2004

March 2012 update

The routine Interp_3 below is a modification of interp_1, which can deal with gaps on the data series. Adapted by Anna Ripoll and Jose Jimenez in March 2012.

Function Interp_3 (xwave1, ywave1, xwave2, maxdur_secs) // Regridding a Wave with gaps into the Time Axis of Another Wave without gaps

Variable maxdur_secs // We need to create a variable which define the maximum duration of the gap in seconds

Wave xwave1, ywave1, xwave2 // xwave1 is the date/time of the initial wave, ywave1 is the data of the initial wave, xwave2 is the date/time of the final wave and ywave2 is not used

variable i1,i2 // i1 is the point from xwave1 and i2 is the point from xwave2

duplicate/o xwave2, ywave_regrid2 // ywave_regrid2 is the new data wave

ywave_regrid2 = NaN // it just cleans the ywave_regrid2 in case there is anyone already

ywave_regrid2 = ywave1 [binarysearchinterp(xwave1, xwave2[p])] // it regrid the ywave1 into the xwave2

// Now we want to stop the regridding process when there is a gap, so we need to filter out

for(i1=0; i1<NumPnts(xwave1)-1; i1 += 1) // For all points of the xwave1

If ((xwave1[i1+1]-xwave1[i1])>=maxdur_secs) // If the xwave1 points are further from each other than the maximum duration gap we created

for(i2=0; i2<NumPnts(xwave2)-1; i2 += 1) // For all the xwave2 points which are inside the gap

If ((xwave2[i2]>(xwave1[i1]+maxdur_secs/2)) && (xwave2[i2]<(xwave1[i1+1]-maxdur_secs/2))) // If they are far from the bound (you can choose how far, in this case maxdur_secs/2)

ywave_regrid2[i2]=NaN // then we decide don't have that data

Endif

EndFor

EndIf

endfor

appendtotable ywave_regrid2

End Function

Looking for an Object in Experiments in a Directory

Function findThisObjectInExpInPath(ObjectType, ObjectName, PathName, ThisSubFolder)

//imu 3 Sept 2009

// This function will help you look for lost objects in your experiments.

// ObjectType must be "string", "variable", or "wave".

// ObjectName must be the exact name of the object, and can be semicolon-separated list.

// Before running this function you must create an Igor Path to the folder you want to search in.

// Command line: NewPath ThisPathName "C:Program Files:WaveMetrics:Igor Pro Folder:"

// or Misc -> New Path

// Send the pathName as a string to the function, e.g, "ThisPathName"

// ThisSubFolder refers to a subfolder to check in the experiments you're looking in.

// To look in root: use "".

// To look in all folders, use "allFolders". This is a cheat for selecting the /R flag for LoadData.

// Note that if you look in allFolders, the complete directory structure of the opened file will be

// created in the current experiment, even if the object is not found.

//

// You can further customize this function for overwriting old data and storing loaded data in a new dataFolder

// by including the /O and /T flags, respectively, with LoadData

string ObjectName, ObjectType, PathName, ThisSubFolder

variable idex, LoadFlag

// check that ObjectType is valid and set LoadFlag for LoadData

// Note that LoadData could look for multiple objectTypes at once, but I'm not using it that way

if(stringmatch(ObjectType, "string"))

LoadFlag = 4

elseif(stringmatch(ObjectType, "variable"))

LoadFlag = 2

elseif(stringmatch(ObjectType, "wave"))

LoadFlag = 1

else

abort "ObjectType must be \"string\", \"variable\", or \"wave\"."

endif

// get ; -separated list of file names in PathName

string ExperimentsInFolderList = indexedFile($PathName, -1, ".pxp")

For(idex = 0; idex < ItemsInList(ExperimentsInFolderList); idex += 1)

// Searching in a Subfolder or in root:

if(strlen(ThisSubFolder) == 0) // root

LoadData/J=objectName /L=(loadFlag)/P=$pathName StringFromList(idex,ExperimentsInFolderList)

elseif(stringmatch(ThisSubFolder, "allFolders")) // recursive look in all subfolders

LoadData/J=objectName/L=(loadFlag)/P=$pathName/R StringFromList(idex,ExperimentsInFolderList)

else // subfolder provided

LoadData/J=objectName /L=(loadFlag)/P=$pathName/S=ThisSubFolder StringFromList(idex,ExperimentsInFolderList)

endif

endfor

end

Add experiment info to graph

Use this simple function to add experiment and graph info to any graph. It can be useful for talks to later find the igor experiment and graph that was used to make slides.

Function AddExpinfo2Graph()
  pathinfo home
  String message = "\\Z06Location:"+S_path
  message +="\rname: \\{igorinfo(1)}"
  message +="\rdate: \\{date()}, \\{time()}"
  message +="\rgraphname: \\{winlist(\"*\",\"\",\"WIN:\")}"
  TextBox/C/N=expinfo/F=0/A=RT/E=2 message
End Function

Added by Harald Stark, Oct 4th, 2013

Other Tricks

• Making a radical in an axis label: \B\S\F'\F'Wingdings'l

RK4 and Euler's Methods

These functions are adapted from Numerical Recipes and allow the integration of a system of ordinary differential equations in time using the Euler method and the 4th Order Runge-Kutta Method. Adapted by Christoph Knote (NCAR) and edited by Jose-Luis Jimenez on Feb. 2013.

A pretty flexible decile binning function

Download it here Note: Need to delete all NaNs first, and also need to be sure data is on the same time grid.

Bring to front all windows with a given wave

Howard Rodstein - wavemetrics

Function DisplayGraphsAndTablesWithWave(w)

Wave w

String list = WinList("*", ";", "WIN:3") // List of all graphs and tables

Variable numWindows = ItemsInList(list)

Variable i

for(i=0; i<numWindows; i+=1)

String windowName = StringFromList(i, list)

CheckDisplayed /W=$windowName w

if (V_flag != 0)

DoWindow /F $windowName

endif

endfor

End

Pedro's general utilities

Last Revision Date: 2/25/2014 (v0.95). Download the ipf here

Functions included:

• function NaNMe(ZWave, Code) Removes Zeros (Code=0) and ICARTT Codes (Code=1) or both (Code=2)
• function/WAVE TieDyeMe(BinMax, BinNumber, Hue, CWaveName) Returns a 3 point wave called CWaveName specifying a the rgb code for a color. Colors are linearly interpolated with index BinNumber and max BinMax on a rainbow colorscale. So same effect as coloring a dataset by f(z), except this applies to multiple datasets along a linear dimension.
• function GenStopWave(YEnd, TimeCode, WYStart, CLength) Utility Function for Resample/Desample functions.
• function NaNCheck(MyWave, Tout) Checks for NaNs in a Wave. Returns 1 for valid values, 0 for NaNs
• function NaNCheck2D(My2DWave, IndexN, Dim, Tout) Same as NaNCheck for 2D Waves, will return NaN if one value is not valid
• function NaNCheckAll2D(My2DWave, IndexN, Dim, Tout) Returns 1 if there is at least 1 valid value, NaN otherwise
• function/WAVE SumRows(TDWave, LRow, RRow, Result) Adds all rows of a 2D Wave between indices LRow and RRow and returns the a 1D Wave with name Result
• function/WAVE AveColumns(TDWave, LRow, RRow, Result) Averages all columns of a 2D Wave between indices LRow and RRow and returns the a 1D Wave with name Result
• function Display2Dvs2D(WaveY, WaveX, WaveLabel) Plot 2 dimensional waves as pairs of 1D waves against each other in a stacked plot
• function/WAVE ReSampleMe(WXStart, WXEnd, WXVar, WWYStart, YEnd, WYVar, ZVar, MinOLap, TimeCode, QCFlag) Box Averages a faster trace (YVar) onto a slower timescale (XVar, can be NaN) and outputs the wave onto ZVar. If less than MinOLap in one X bin is comprised of valid Y values a NaN or zero (depeding on global constant NaNResample) will be returned. TimeCode follows the ICARTT convention (1=StartUTC, 1s data, -1:StopUTC,1s data, 0:Start and StopUTC required). QCFlag=1 will keep the wave were the overlap of both timedomains is calculates, useful for diagnostics.
• function/WAVE DeSampleMe(WXStart, WXEnd, WXVar, WWYStart, YEnd, WYVar, ZVar, TimeCode) Opposite of Resample, a slow wave (X) gets mapped (unsmoothed) onto a faster wave (Y).
• function/WAVE TBoxSmooth(TStart, TEnd, MyVar, SmthVar, TSmooth, MinOLap, TimeCode) Boxsmooths Data by succesively desampling and resampling it back to the original timescale.
• function NaNMean(MyWave, LowL, HighL) Compute the average of a NaN containing wave between indices LowL and HighL.
• function NaNSum(MyWave, LowL, HighL) Compute the sum of a NaN containing wave between indices LowL and HighL.
• function/WAVE CleanNaNs(MyTime, MyWave, suffix, ReMixWave) Removes NaNs from both MyTime and MyWave, saves the results in Waves Name(MyTime)+Suffix and Name(MyWave)+Suffix and return the later wave. The original positions of the NaNs are saved in the wave RemixWave.
• function/WAVE MCleanNaNs(MyWave, suffix, ReMix) Remove NaNs from a wave according to a previously computed Remix wave (ie this can be used for cleaning up multiple waves of the same length based on one "master" wave)
• function ReBuildNaNWave(NaNWave, MyTime, MyWaveN, ReMix) Inserts the NaNs back into a previously DeNaN'd wave. The original Timewave (pre-NaNClean) is only used for correctly dimensioning the output wave. Unlike CleanNaNs, the TimeWave needs to be Re'NaN'd by as separate function call.
• function StackMyImage(MyImage, MyAxis, MyMask, DateFlag): Will generate a rainbow-coloured set of traces from a 2DMatrix with the Dimensions MyAxis and MyMask. Use NaNs in MyMask to limit which parts of MyImage are plotted. A legend with the contents of MyMask is added, set DateFlag=1 if MyMask is a Timewave.
• function SmoothCorrs(MyTime, DataA, DataB, SExt, SmoothT, SmoothC, CInterval, MaxSmooth, SmoothType) Computes the correlation of 2 Datasets as a function of the length of the coaveraging interval. Outputs the length of the averaging interval into Wave SmoothT and the r2 for each step into SmoothC. CInterval and MaxSmooth determine the step and maximum length of the averaging interval. SExt is the extension added to DataA and DataB for the averaged results. Smoothtype=0 chooses TBoxSmooth as the averaging function (pure Boxcar), Smoothtype=1 chooses a running average smooth instead
• function OldTBoxSmooth(MyTime, MyData, SData TimeStep, SmoothType) Smoothes MyData=f(MyTime) with a constant TimeStep and outputs result into SData. Smoothtype=1 is a running average, 0 does a boxcar resample. Use if you really need the running average capability, otherwise TBoxsmooth should be more efficient/better mantained.
• function DiurnalCycle(SStartTime, StopTime, MyWave, TimeCode, MyDay, MyCycle, DevType) Computes a diurnal average of MyWave into MyCycle according to SStartime and Stoptime. MyDay has the hour of the day timwave, TimeCode follows the ICARTT convention (1=StartUTC, 1s data, -1:StopUTC,1s data, 0:Start and StopUTC required), Devtype defines what variance measure to output (0=none, 1 = stddev, 2 = stderr)
• function WaterP(Temp) Calculates Water Pressure in mbars as function of Temperature, useful for RH computations. Temp in degrees C
• function ParseTimeDataStr(TDate) Will return a TimeDate Value from an inputed TimeString
• function GetMyParams(VarId, VarType, FileId, SFolder, AppendMe, WOutput) Will retrieve 1D waves (Vartype=1) or nvar(Vartype=2)/svars (Vartype=4) from all Igor experiments in a chosen directory and import them into the current experiment. VarId is the name of Igor object, FileId is a temporary nvar containing the file handle of the current experiment being imported. SFolder specifies a Subfolder that exists in the experiments to be parsed. AppendMe specifies for 1D wave object (Vartype=1) if the waves should be concatenated into a 1D Wave (AppendMe=1) or if a 2D Wave with Columns=Experiments should be created. WOutput specifies the Name/Path of the target wave in which the imported data will be written to. DeleteMe will delete the temporary wave of name VarId at the end of the import process. MAKE SURE that WOutput != VarId, or you will lose your data.
• function SetUpBinning (BinFolder, BinWave, BinRes) Will create the necessary Indexwaves to make histograms of any number of variables (to be defined in text wave BinEm) as a function of BWave. If BinType=1, BinRes will specify the absolute step in BWave per Bin. If BinType=0, BinRes will specify the number of bins (ie 10 for deciles)
• function BinEm(BinFolder, BinWave, BinRes, WavesToBin) Create Histograms of all the waves listed in WavestoBin by the variable BWave. Binning for BWave has to be defined previously by SetupBinning. Use the TextWave NamesToBin to specify abbreviated names for the output waves (or otherwise use WavesToBin again). If FullOutput=1 it will also create stdev and Npnts waves for each histogram