/****************************************************************************/ /* */ /* IAGA_check_data.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program removes various kinds of spikes and jumps in IMAGE data */ /* files and replaces these erronuous data values by missing values. */ /* */ /* Usage: */ /* IAGA_check_data [-s] [-e] [-p] [-c] [-o] [-i] [-v] [-h] */ /* [-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. */ /* [-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. */ /* 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.0 10.11.2003 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.0 10.11.2003 First official release */ /****************************************************************************/ #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" #include "Spike.h" char *version = "1.0"; char *date = "10.11.2003"; #define HeaderLineCount 5 #define UsageLineCount 19 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", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-s] [-e] [-p] [-c] [-o] [-i] [-v] [-h] [<] 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. ", " [-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. ", " 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[200] = ""; /* 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[200] = ""; /* 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 = 300; /* Maximum length of spike = 1 day */ long GapLength = 0; /* Maximum gap length which is */ /* still interpolated */ long TriggerValue = 300; /* unit = 0.1 nT */ long AmplBefore = 10; long AmplAfter = 3; float SlopeParam = 1.0; /* nT/min */ long SpikeFound = 0; /* Flag for found spikes */ /*==========================*/ /* 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 an 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]); } SpikeFound = 1; } 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 (SpikeFound); }