/****************************************************************************/ /* */ /* IAGA_IL_index.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program computes the IL index for given stations. */ /* */ /* Usage: */ /* IAGA_IL_index [-s] [-e | -h] [-o ] [-p] [-b] [-c] [-f] [-g] */ /* [<] IAGAFile > IL_index_file */ /* [-s YYYYMMDDHH] Time of the first record included. If missing */ /* then start of file is assumed. */ /* [-e YYYYMMDDHH] Time of the record not included anymore. If */ /* missing then end of file is assumed. */ /* [-h HH] Number of hours included. Either -e or -h can */ /* be specified, not both. */ /* [-o 'Station list'] List of stations delimited by aposthropes. */ /* If missing then all stations included. */ /* Stations are identified by three-letter code */ /* and separated by exactly one space. */ /* Format : 'SJG KAK HON HER' */ /* |-b BaseLineInfo] Determines the baselines used in the plots. */ /* BaseLineInfo may have the following two values: */ /* 1. H1-H2 The baseline is determined as the */ /* average field value from hour H1 to hour H2 */ /* (H2 is not included). The hour is counted */ /* from the start time of the plot. */ /* E.g -s 94120106 -b 6-7 implies that the */ /* baseline is the one hour average from */ /* 94120112 to 94120113. Note that H1 and H2 */ /* may be either positive or negative (e.g. */ /* -b -5--3 ) */ /* 2. QN The baseline is determined as the */ /* average of the quietest N hour interval in */ /* the dataplot. E.g -b Q2 implies that the */ /* program computes the max-min values for the */ /* intervals 0-2,1-3,2-4, ... and uses the */ /* interval with smallest variations to compute */ /* the baseline value. */ /* If -b is missing then the baseline is set to */ /* the average value of the field over the whole */ /* period of the plot. */ /* [-c] Compute and write into output the contributions */ /* [-g] Write into output the latitudes and longitudes */ /* of stations which determine each IL and IU index*/ /* [-f HeaderFormat] Format of the header. 1 = old format (default) */ /* 2 = new format with baseline time info etc. */ /* of various stations into the IL and IU indices. */ /* [-a] A special option designed for Eija Tanskanen */ /* to write time information as minutes since */ /* the start of the year. */ /* [-v] Don't write header information into output file */ /* IAGAFile Name of IAGA-format file. */ /* IL_index_file Name of IL_index_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.04 10.01.2020 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.04 10.01.2020 Date strings must be YYYYMMDD (four digits for year) */ /* Replaced StrToSecs -> StrToSecsC and */ /* SecsToStr -> SecsToStrC. */ /* 1.03 15.05.2019 Added -f and -g options */ /* 1.02 13.12.2018 Added baseline period times to header of output file */ /* 1.01 08.10.2007 Added -a and -v options */ /* 1.0 12.09.2006 First official release */ /****************************************************************************/ #include #include #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" char *version = "1.03"; char *date = "15.05.2019"; #define HeaderLineCount 4 #define UsageLineCount 40 char *HeaderText[HeaderLineCount] = { "**************************************************************************", " This program computes the IL, IU and IE indices for given time interval ", " and given stations. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-s] [-e | -h] [-o] [-p] [-b] [-c] [<] IAGA_file > IL_index_file", " [-s YYYYMMDDHH] Time of the first record included. If missing ", " then start of file is assumed. ", " [-e YYYYMMDDHH] Time of the record not included anymore. If ", " missing then end of file is assumed. ", " [-h HH] Number of hours included. Either -e or -h can ", " be specified, not both. ", " [-o 'Station list'] List of stations enclosed in quotes. ", " [-p] Use this if data is read from pipeline ", " [-b BaseLineInfo] Determines the baselines used in the plots. ", " BaseLineInfo may have the following two values: ", " 1. H1-H2 The baseline is determined as the ", " average field value from hour H1 to hour H2 ", " (H2 is not included). The hour is counted ", " from the start time of the plot. ", " E.g -s 94120106 -b 6-7 implies that the ", " baseline is the one hour average from ", " 94120112 to 94120113. ", " H1 and H2 may be either positive or negative.", " (e.g. -b -5--3). ", " 2. QN The baseline is determined as the ", " average of the quietest N hour interval in ", " the dataplot. E.g -b Q2 implies that the ", " program computes the max-min values for the ", " intervals 0-2,1-3,2-4, ... and uses the ", " interval with smallest variations to compute ", " the baseline value. ", " If -b is missing then the baseline is set to ", " the average value of the field over the whole ", " period of the plot. ", " [-c] Compute and write into output the contributions ", " of various stations into the IL and IU indices. ", " [-g] Write into output the latitudes and longitudes ", " of stations which determine each IL and IU index.", " [-f HeaderFormat] Format of the header. 1 = old format (default) ", " 2 = new format with baseline time info etc. ", " [-v] Don't write header (station info etc.) into ", " output file ", " IAGA_file Name of input IAGA format file. ", " IL_index_file Name of IL index file. " }; /*--------------------------------------------------------------------------*/ /* Get the start and end hours from the BaseStr. The format of BaseStr is */ /* HH-HH. */ /*--------------------------------------------------------------------------*/ void GetHours(char *BaseStr,long *t0, long *t1) { long i = 0; long t = 0; long negative = 0; if (BaseStr[0] == '-') { negative = 1; i++; } while ((BaseStr[i] != '-') && (i < strlen(BaseStr))) { t = 10*t + (BaseStr[i]-48); i++; } if (negative) t = -t; *t0 = t; t = 0; negative = 0; if (BaseStr[++i] == '-') { negative = 1; i++; } while (i < strlen(BaseStr)) t = 10*t + (BaseStr[i++]-48); if (negative) t = -t; *t1 = t; } /*--------------------------------------------------------------------------*/ /* Get the hour count from the BaseStr. The format of the string is QHH. */ /*--------------------------------------------------------------------------*/ long GetQHour(char *BaseStr) { long i = 0; long t = 0; while (++i < strlen(BaseStr)) t = 10*t + (BaseStr[i]-48); return t; } /*--------------------------------------------------------------------------*/ /* Find the baseline value for given station and given field component. */ /* If the baseline value is not defined in the parameter file then compute */ /* the value as the average value over the whole dataplot period. */ /*--------------------------------------------------------------------------*/ void ComputeBaselines(NetworkPtr IMAGE, Time_sec StartTime,Time_sec EndTime, char *BaseStr, Time_sec *baseStart, Time_sec *baseEnd) { long Base = -1; Time_sec t,t0,t1; long Hours,Max,Min; long h1,h2; long MaxMinSum; /* Sum of Max - Min values over all stations */ long StatCount; /* Number of stations included in MaxMinSum */ long Good; long CurrentMin; StationPtr s; /* Check the BaseStr (-b option) and compute the baseline */ /* accordingly. */ if (strlen(BaseStr) != 0) { if (BaseStr[0] == 'Q') { /* Quiet period average */ Hours = GetQHour(BaseStr); t0 = StartTime; CurrentMin = 999999; for (t = StartTime;tStationList; MaxMinSum = 0; StatCount = 0; while (s != NULL) { FindMaxMin(s,'X',t,3600*Hours,&Max,&Min); Good = FindPercentage(s,'X',t,3600*Hours); if ((Max != MissingValue) && (Good > 80)) { MaxMinSum += Max-Min; StatCount++; } s = s->Next; } if (StatCount > 0) { MaxMinSum = MaxMinSum/StatCount; if (MaxMinSum < CurrentMin) { CurrentMin = MaxMinSum; t0 = t; } } } t1 = t0 + 3600*Hours; } else { /* Average over given hours */ GetHours(BaseStr,&h1,&h2); t0 = StartTime + h1*3600; t1 = StartTime + h2*3600; } } if (strlen(BaseStr) == 0) { t0 = StartTime; t1 = EndTime; } *baseStart = t0; *baseEnd = t1; /* Finally compute baselines as averages over determined time interval */ s = IMAGE->StationList; while (s != NULL) { s->Xbase = ComputeAverage(s,'X',t0,t1-t0); s->Dbase = 0; /* Dbase is used in computation of station weights in IL index */ s->Fbase = 0; /* Fbase is used in computation of station weights in IU index */ s = s->Next; } } /*--------------------------------------------------------------------------*/ /* Write the station names in one line and baseline values in another line */ /* into stdout. */ /* The first character of the line is a percent sign '%' so that these */ /* lines are interpreted as command lines in MatLab. */ /*--------------------------------------------------------------------------*/ void WriteHeader(NetworkPtr IMAGE, Time_sec baseStart, Time_sec baseEnd, long HeaderFormat) { StationPtr s; time_t t; /* time in seconds since 1900 Jan 1st */ struct tm TM; /* time in years, months,.., secs */ long count; char TimeStartStr[20]; char TimeEndStr[20]; printf("%% The following %d stations are used:\n",IMAGE->StationCount); /* Station names */ s = IMAGE->StationList; printf("%%"); while (s != NULL) { printf(" %s",s->StationID); s = s->Next; } printf("\n"); /* Contributions to IU */ if (HeaderFormat == 2) { printf("%% Station contributions to IU and IL:\n"); } s = IMAGE->StationList; count = (s->EndTime - s->StartTime)/(s->TimeStep) + 1; printf("%%"); while (s != NULL) { printf(" %4d",s->Fbase); s = s->Next; } printf("\n"); /* Contributions to IL */ s = IMAGE->StationList; printf("%%"); while (s != NULL) { printf(" %4d",s->Dbase); s = s->Next; } printf("\n"); /* baselines */ if (HeaderFormat == 2) { printf("%% Station baseline values:\n"); } s = IMAGE->StationList; printf("%%"); while (s != NULL) { printf(" %d",s->Xbase); s = s->Next; } if (HeaderFormat == 2) { printf("\n"); /* baseline times */ printf("%% Baseline period:\n"); printf("%% "); // if (baseStart < YEAR2000) // printf("19"); // else // printf("20"); SecsToStrC(baseStart,TimeStartStr); printf("%s - ", TimeStartStr); // if (baseEnd < YEAR2000) // printf("19"); // else // printf("20"); SecsToStrC(baseEnd,TimeEndStr); printf("%s", TimeEndStr); } printf("\n\n"); /* File creation time */ time(&t); /* get current time in seconds */ TM = *localtime(&t); /* convert into years,...,secs */ printf("%% File created on "); printf("%02d-%02d-20%02d ",TM.tm_mday,TM.tm_mon+1,TM.tm_year-100); printf("%02d:%02d:%02d\n\n",TM.tm_hour,TM.tm_min,TM.tm_sec); } /*--------------------------------------------------------------------------*/ /* Compute the IL index from stations in structure IMAGE and store the data */ /* into structure IL. */ /*--------------------------------------------------------------------------*/ void Compute_IL_index(NetworkPtr IL_index, NetworkPtr IMAGE) { StationPtr w,s; StationPtr vL,vU; /* Stations with lowest (L) and uppermost (U) field values */ Time_sec t; long IU,IL,XX; long station_index; long IU_station_index; long IL_station_index; long IU_latitude; long IL_latitude; w = IL_index->StationList; t = w->StartTime; while (t < w->EndTime) { IU = -999999; IL = 999999; s = IMAGE->StationList; IU_latitude = s->Latitude; IL_latitude = s->Latitude; station_index = 1; while (s != NULL) { if (s->Xbase != MissingValue) { XX = GetDataValue(s,t,'X'); if (XX != MissingValue) { XX -= (s->Xbase); if (XX > IU) { IU = XX; vU = s; IU_station_index = station_index; } else if (XX == IU) { // Select southern station if (s->Latitude < IU_latitude) { IU_latitude = s->Latitude; vU = s; IU_station_index = station_index; } } if (XX < IL) { IL = XX; vL = s; IL_station_index = station_index; } else if (XX == IL) { // Select southern station if (s->Latitude < IL_latitude) { IL_latitude = s->Latitude; vL = s; IL_station_index = station_index; } } } } s = s->Next; station_index++; } if ( (IU != -999999) && (IL != 999999) ) { /* if (IU < 0) IU = 0; */ /* if (IL > 0) IL = 0; */ SetDataValue(w,t,'X', IL); SetDataValue(w,t,'Y', IU); SetDataValue(w,t,'Z', IU-IL); SetDataValue(w,t,'F', 100*IL_station_index + IU_station_index); SetDataValue(w,t,'H', 10000*vU->Latitude + vU->Longitude); // IU station location SetDataValue(w,t,'D', 10000*vL->Latitude + vL->Longitude); // IL station location vU->Fbase++; /* Fbase is used in computing station weights of IU */ vL->Dbase++; /* Dbase is used in computing station weights of IL */ } else { SetDataValue(w,t,'X', MissingValue); SetDataValue(w,t,'Y', MissingValue); SetDataValue(w,t,'Z', MissingValue); SetDataValue(w,t,'F', 9999); SetDataValue(w,t,'H', 99999999); SetDataValue(w,t,'D', 99999999); } t += w->TimeStep; } } /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; /* Dummy index variable */ long status = 0; /* Status variable for file operations */ long FileCount = 0; /* Number of files in command line */ long HourCount = 0; char FileName[100] = ""; /* Name of IAGA file */ char StationStr[200] = ""; /* List of stations to be processed */ char StartTimeStr[20] = ""; char EndTimeStr[20] = ""; Network_struct IMAGE; /* List of station structures */ Network_struct IL_index; /* Structure containing IL index data */ StationPtr s,w; /* Pointers to station structures */ StationInfoPtr p; Time_sec t,StartTime,EndTime; char TimeStr[20]; char BaselineStr[20] = ""; long ComputeStatistics = 0; /* Compute and write the contributions */ /* of stations into output */ long MinuteOption = 0; /* Write time as minutes from the start */ /* of the current year. Also write only */ /* IL index. */ long FirstSec; /* First second of the current year */ Time_struct MyTime = {0,0,0,1,0,0}; long HeaderInfo = 1; /* Write header info into output file */ long HeaderFormat = 1; /* Header format, either 1 or 2 */ long GeoInfo = 0; /* Write lon and lat of StationStr */ Time_sec BaselineStartTime, BaselineEndTime; long latlon; /*==========================*/ /* 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' : strcpy(StartTimeStr,argv[++params]); break; case 'e' : strcpy(EndTimeStr,argv[++params]); break; case 'h' : HourCount = atol(argv[++params]); break; case 'b' : strcpy(BaselineStr,argv[++params]); break; case 'o' : strcpy(StationStr,argv[++params]); break; case 'c' : ComputeStatistics = 1; break; case 'f' : HeaderFormat = atol(argv[++params]); break; case 'g' : GeoInfo = 1; break; case 'a' : MinuteOption = 1; break; case 'v' : HeaderInfo = 0; break; case 'p' : 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 the values of the parameters */ /*====================================*/ if (FileCount > 1) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if ((HourCount != 0) && (strlen(EndTimeStr) > 0)) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if ((HourCount != 0) && (strlen(StartTimeStr) == 0)) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if (HourCount != 0) { SecsToStrC(StrToSecsC(StartTimeStr,0)+3600*HourCount,EndTimeStr); } /*====================================================*/ /* Read the data from an IAGA format file into memory */ /*====================================================*/ status = ReadIAGA(FileName,&IMAGE,NULL,NULL,StationStr); 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; } /*==================================================*/ /* Create the network structure for IL index data */ /*==================================================*/ s = IMAGE.StationList; InitNetwork(&IL_index); p = FindStationInfo("XXX"); NewStation(&IL_index,p,s->StartTime,s->EndTime,s->TimeStep,s->TempStep); w = IL_index.StationList; strcpy(w->StationID,"IL_"); w->Latitude = 9000; w->Longitude = 0; strcpy(w->Components,s->Components); if (strlen(StartTimeStr) == 0) StartTime = w->StartTime; else StartTime = StrToSecsC(StartTimeStr,0); if (strlen(EndTimeStr) == 0) EndTime = w->EndTime; else EndTime = StrToSecsC(EndTimeStr,0); /*==================================*/ /* Compute the IL index */ /*==================================*/ ComputeBaselines(&IMAGE,StartTime,EndTime,BaselineStr, &BaselineStartTime, &BaselineEndTime); Compute_IL_index(&IL_index,&IMAGE); if (HeaderInfo == 1) { WriteHeader(&IMAGE, BaselineStartTime, BaselineEndTime, HeaderFormat); } /*======================================================*/ /* Write the IL index data into stdout */ /*======================================================*/ t = StartTime; if (MinuteOption == 1) { SecsToTm(t,&MyTime); MyTime.tm_sec = 0; MyTime.tm_min = 0; MyTime.tm_hour = 0; MyTime.tm_mday = 1; MyTime.tm_mon = 0; FirstSec = TmToSecs(&MyTime); /* First second of current year */ } while (t < EndTime) { /* Print time */ SecsToStrC(t,TimeStr); // YYYYMMDD if (MinuteOption == 0) { // if (t >= YEAR2000) // printf("20"); // else // printf("19"); // printf("%.2s %.2s %.2s %.2s %.2s %.2s ",TimeStr,TimeStr+2, // TimeStr+4,TimeStr+6,TimeStr+8,TimeStr+10); printf("%.4s %.2s %.2s %.2s %.2s %.2s ",TimeStr,TimeStr+4, TimeStr+6,TimeStr+8,TimeStr+10,TimeStr+12); /* Print indices IL, IU, IE */ printf(" %7.1f",GetDataValue(w,t,'X')/10.0); printf(" %7.1f",GetDataValue(w,t,'Y')/10.0); printf(" %7.1f",GetDataValue(w,t,'Z')/10.0); if (ComputeStatistics == 1) { printf(" %04d",GetDataValue(w,t,'F')); } if (GeoInfo == 1) { latlon = GetDataValue(w,t,'D'); // IL printf(" %5.2f %5.2f", (double)(latlon/10000)/100.0, (double)(latlon%10000)/100.0); latlon = GetDataValue(w,t,'H'); // IU printf(" %5.2f %5.2f", (double)(latlon/10000)/100.0, (double)(latlon%10000)/100.0); } } else { printf("%8d",1440+(t - FirstSec)/60); printf(" %7.1f",GetDataValue(w,t,'Y')/10.0); printf(" %7.1f",GetDataValue(w,t,'X')/10.0); } printf("\n"); t += w->TimeStep; } FreeNetwork(&IMAGE); FreeNetwork(&IL_index); return OK; }