/****************************************************************************/ /* */ /* IAGA_remove_spikes.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program removes various kinds of spikes in IMAGE data files and */ /* replaces these erronuous data values by missing values. */ /* */ /* Usage: */ /* IAGA_remove_spikes [-s] [-e] [-p] [-d] [-a] [-b] [-r] [-c] [-i] [-v] */ /* [-o] [-h] [-t] [<] IAGAFile > NewIAGAFile */ /* [-s StartTime] First record to be searched for spikes. */ /* If missing then start of file is assumed. */ /* [-e EndTime] Last record to be searched for spikes. */ /* If missing then end of file is assumed. */ /* [-p] This is used when reading data from pipe */ /* [-c CompList] List of components to be searched. If missing */ /* then all components are searched. */ /* [-o StationList] List of three letter station ID's to be */ /* searched. If missing then all stations in the */ /* data file are searched. */ /* [-d Trigger] Trigger value (nT). A field jump is a possible */ /* spike if |B(t)-B(t-dt)| > Trigger */ /* Default value for Trigger is 5.0 nT. */ /* However, also general field variation is taken */ /* into account when considering whether a jump is */ /* a spike or not. So Trigger gives only the */ /* limit for a possible spike. Field jumps less */ /* than Trigger are never considered as a spike. */ /* Option -a controls how the field variations are */ /* taken into account. */ /* [-a AmplFactor] AmplFactor controls how field variations are */ /* taken into account when searching spikes. */ /* The method is to look at differences between */ /* successive field values and comparing them to */ /* average field differences of 10 previous points.*/ /* If the difference is larger than AmplFactor* */ /* (average difference) then we have a possible */ /* spike. Default value for AmplFactor is 10. */ /* [-b AmplFactor2] AmplFactor2 is used in searching the end of */ /* spike structure. It defines the factor by how */ /* much the difference of two good data points */ /* after the spike may differ from the average */ /* difference (see -a option) before the spike. */ /* The default value for AmplFactor2 is 3. It is */ /* thus assumed that the activity of field is */ /* about the same before and after the spike. */ /* [-r FieldSlope] FieldSlope describes how quickly the field is */ /* allowed to change during the spike. Unit of */ /* FieldSlope is nT/min and the default value is 9.*/ /* [-v] Verbose. Display information about spikes which */ /* are found and removed. */ /* [-h] Don't write the header lines. */ /* [-i] GapLength Interpolate gaps which are produced by spike */ /* removals and whose length is less than or equal */ /* to GapLength. */ /* [-t] SpikeLength Maximum time length of the spike structure in */ /* seconds. If spike structure seems to be longer */ /* than SpikeLength then nothing is removed. */ /* IAGAFile Name of IAGA-format daily file. */ /* NewIAGAFile Name of corrected IAGA-format file. */ /* */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Geophysical Research Division */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.14 14.8.2006 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.14 14.08.2006 Fixed Spike.h to handle correctly situations where */ /* missing data was encountered. */ /* 1.13 09.09.1999 Fixed a Y2K bug in NewTime.h file which resulted in */ /* incorrect year in date strings if year >= 2000. */ /* 1.12 19.05.1999 Added -h option. */ /* 1.11 06.11.1998 Added -a, -b, -r and -t options. */ /* 1.0 10.03.1997 First official release */ /****************************************************************************/ #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" #include "Spike.h" char *version = "1.14"; char *date = "14.8.2006"; #define HeaderLineCount 20 #define UsageLineCount 48 char *HeaderText[HeaderLineCount] = { "**************************************************************************", " This program locates and removes spikes and similar structures in IAGA ", " format data files and replaces these erronuous data values with missing ", " values. If desired, it will fill the resulting data gaps by interpolation", " The algorithm for searching spikes consists of two parts and both are ", " controlled by two parameters: ", " 1. Locate spike start. Spike start is located by looking at differences ", " of successive field values. When the difference is larger than given ", " trigger value (-d option) and larger than average field difference ", " of 10 previous points by given factor (-a option) then we have found ", " a spike structure. ", " 2. Locate spike end. The spike end is found by imposing two conditions. ", " First, the field difference of two field values immediately after the ", " spike must be within a given factor (-b option) of the average field ", " difference before the spike (see above). Second, the first field value", " after the spike must be reasonably close to the field value just ", " before the spike. This is controlled by 'slope' parameter (-r option) ", " which describes how much the field is allowed to change during the ", " spike. The unit of slope parameter is nT/min. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-s] [-e] [-d] [-a] [-b] [-r] [-c] [-o] [-i] [-v] [-h] [-p] [-t]", " [<] IAGAFile > NewIAGAFile ", " [-s YYMMDDHHMMSS] First record to be searched for a spike. ", " If missing then start of IAGAFile is assumed. ", " [-e YYMMDDHHMMSS] Last record to be searched for a spike. ", " If missing then end of IAGAFile is assumed. ", " [-d TriggerValue] Trigger value (nT). A field value is a possible ", " spike if |B(t)-B(t-dt)| > TriggerValue. Default ", " value of TriggerValue is 5.0 nT. However, general", " field variations (-a option) are also taken into ", " account when deciding whether a field jump is a ", " spike or not. TriggerValue gives only the lowest ", " limit for a field jump to be considered as a ", " spike. Jumps smaller than TriggerValue are never ", " considered as spikes. TriggerValue is a float ", " number so it may have decimals. ", " [-a AmplFactor] AmplFactor is the factor by how much a field jump", " must exceed the average difference of 10 previous", " field values before it can be regarded as a ", " spike. The default value is 10. ", " [-b AmplFactor2] AmplFactor2 is used in searching the end of spike", " structure. It defines the factor by how much the ", " difference of two good data points after the ", " spike may differ from the average difference (see", " -a option) before the spike. The default value ", " for AmplFactor2 is 3. It is thus assumed that the", " activity of field is about the same before and ", " after the spike. ", " [-r FieldSlope] FieldSlope describes how quickly the field is ", " allowed to change during the spike. Unit of ", " FieldSlope is nT/min and the default value is 9. ", " [-c Components] Spikes are searched only in these components. ", " If missing then all components are searched. ", " [-o Stations] List of three letter station ID's. Only these ", " stations are searched. If missing then all ", " stations are searched. ", " [-i GapLength] Interpolate all data gaps caused by spike removal", " which are shorter or equal to GapLength. ", " [-v] Verbose. Display information about spikes which ", " are found and removed. ", " [-h] Don't write the header lines. ", " [-p] Data is read from pipe. This prevents the info ", " text from appearing. ", " [-t] SpikeLength Maximum length of spike structure in seconds. ", " If the spike seems to be longer than SpikeLength ", " then nothing is removed. ", " IAGAFile Name of IAGA-format file. ", " NewIAGAFile Name of fixed IAGA-format file. ", }; /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; long status = 0; long FileCount = 0; char FileName[100] = ""; /* Name of the original data file */ Network_struct IMAGE; /* Structure containing the data */ /* for all stations. */ StationPtr Station; /* Pointer to the current station */ char StationList[100] = ""; /* List of stations in command line */ char CompStr[5] = "XYZ"; /* List of components */ long Comp; /* Current magnetic component */ long VerboseFlag = 0; /* Display information about gaps */ long HeaderFlag = 1; /* Write the header lines */ Time_sec StartTime = MIN_TIME; /* Time of first data block (secs) */ Time_sec EndTime = MAX_TIME; /* Time of last data block */ Time_sec SpikeStart,SpikeEnd; /* Start and end times of spikes */ Time_sec T; /* Current time */ long SpikeLength = 86400; /* Maximum length of spike = 1 day */ long GapLength = 0; /* Maximum gap length which is */ /* still interpolated */ long TriggerValue = 50; /* unit = 0.1 nT */ long AmplBefore = 10; long AmplAfter = 3; float SlopeParam = 15.0; /* = 9 nT/min */ /*==========================*/ /* Analyse the command line */ /*==========================*/ if (argc == 1) { /* No arguments, show the usage text */ PrintUsage(argv[0],0,HeaderLineCount,UsageLineCount); return OK; } for (params = 1; params < argc; params++) { if (*argv[params] != '-') { strcpy(FileName,argv[params]); FileCount++; } else { switch (*(argv[params]+1)) { case 's' : StartTime = StrToSecs(argv[++params],0); break; case 'e' : EndTime = StrToSecs(argv[++params],0); break; case 'c' : strcpy(CompStr,argv[++params]); break; case 'o' : strcpy(StationList,argv[++params]); break; case 'd' : TriggerValue = 10*atof(argv[++params]); break; case 'a' : AmplBefore = atol(argv[++params]); break; case 'b' : AmplAfter = atol(argv[++params]); break; case 'r' : SlopeParam = (10*atof(argv[++params]))/60; break; case 'v' : VerboseFlag = 1; break; case 't' : SpikeLength = atol(argv[++params]); break; case 'i' : GapLength = atol(argv[++params]); break; case 'p' : break; case 'h' : HeaderFlag = 0; break; default : fprintf(stderr,"\n### %s: \"%s\" is not an option.\n\n", argv[0], argv[params]); PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; break; } } } /*===========================================*/ /* Check that there is at most one data file */ /*===========================================*/ if (FileCount>1) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } /*====================================================*/ /* Read the data from an IAGA format file into memory */ /*====================================================*/ status = ReadIAGA(FileName,&IMAGE,NULL,NULL,NULL); if (status != 0) { if (status == FileError) fprintf(stderr,"Error in reading file %s\n",FileName); if (status == OutOfMemory) fprintf(stderr,"Out of memory while reading file %s\n",FileName); if (status == FileFormatError) fprintf(stderr,"%s is not a IAGA format file\n",FileName); return FAIL; } /*===================*/ /* Write the header */ /*===================*/ if ((VerboseFlag == 1) && (HeaderFlag == 1)) { fprintf(stderr,"-----------------------------------------------------------------------\n"); fprintf(stderr," Obs Comp Start time - End time Length Jump Interpol\n"); fprintf(stderr,"-----------------------------------------------------------------------\n"); } /*======================================*/ /* Go through all stations in NETWORK */ /*======================================*/ Station = IMAGE.StationList; if (StartTime == MIN_TIME) StartTime = Station->StartTime; if (EndTime == MAX_TIME) EndTime = Station->EndTime; while (Station != NULL) { if (StationInList(Station->StationID,StationList)) { /*======================================*/ /* Go through all specified components */ /*======================================*/ for (Comp=0;Comp SpikeLength) { SpikeEnd = SpikeStart + SpikeLength; if (SpikeEnd > EndTime) SpikeEnd = Station->EndTime; } else { if (VerboseFlag == 1) WriteSpikeInfo(Station,CompStr[Comp],SpikeStart, SpikeEnd,10,GapLength); if ((SpikeEnd != Station->EndTime) && (SpikeEnd < EndTime)) RemoveData(Station,CompStr[Comp],SpikeStart,SpikeEnd); if ((SpikeEnd-SpikeStart < GapLength) && (SpikeEnd < EndTime)) InterpolateData(Station,SpikeStart,SpikeEnd,CompStr[Comp]); } } T = SpikeEnd+Station->TimeStep; } /* End of while loop */ } /* End of component loop */ } /* End of if */ Station = Station->Next; } /* End of Station Loop */ /*======================================*/ /* Write corrected data into stdout */ /*======================================*/ WriteIAGA(NULL,&IMAGE,NULL,NULL,NULL); FreeNetwork(&IMAGE); return OK; }