/****************************************************************************/ /* */ /* IAGA_to_Matlab.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This is a filter program that reads a IAGA format data file and writes */ /* out the same data into Matlab binary file. */ /* */ /* Usage: */ /* IAGA_to_Matlab [-s] [-e | -h] [-o] [-i] [-w] [<] IAGAFile > MatlabFile */ /* [-s YYMMDDHH] Time of the first record included. If missing */ /* then start of file is assumed. */ /* [-e YYMMDDHH] 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. */ /* [-i] Field values are stored as 32-bit integers. */ /* If missing then 64-bit double precision used. */ /* Integer files are smaller but the marker for */ /* missing field value is 999999 while on double */ /* precision files it is NaN (= Not a number). */ /* The unit of field strength is 0.1 nT in integer */ /* files and 1 nT in double precision files. */ /* [-p] This may be used to prevent the usage text from */ /* appearing when there are no other parameters */ /* and data is read from pipeline. */ /* IAGAFile Name of IAGA-format file. */ /* MatlabFile Name of Matlab binary 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 9.9.1999 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.04 09.09.1999 Fixed a Y2K bug in NewTime.h file which resulted in */ /* incorrect year in date strings if year >= 2000. */ /* 1.03 22.10.1998 */ /* - Internal modifications due to changes in the MissingValue marker. */ /* 1.02 27.05.1998 */ /* - Increased the number of stations allowed in -o option. Now 200 */ /* characters, earlier 50 characters. */ /* 1.01 26.02.1996 */ /* - Cosmetic change in handling the usage text. No apparent change in */ /* program behaviour. */ /* */ /* 1.0 13.11.1995 First official release */ /****************************************************************************/ #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" char *version = "1.04"; char *date = "9.9.1999"; #define HeaderLineCount 13 #define UsageLineCount 23 char *HeaderText[HeaderLineCount] = { "**************************************************************************", "This program converts an IAGA-format magnetometer data file into Matlab ", "binary file. In Matlab the file is read with the load-command. For each ", "station a (N+2)*3 matrix is created whose name is XXXYYMMDD where XXX is ", "the name of the station and YYMMDD is time of the first data point. N is ", "the number of data points. The structure of the matrix is: ", " Year Month Day ", " Hour Minute Samplerate ", " X1 Y1 Z1 ", " . . . ", " . . . ", " XN YN ZN ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-s] [-e | -h] [-o] [-i] [-p] [<] IAGAFile > MatlabFile", " [-s YYMMDDHH] Time of the first record included. If missing ", " then start of file is assumed. ", " [-e YYMMDDHH] 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. ", " [-i] Field values are stored as 32-bit integers. ", " If missing then 64-bit double precision used. ", " Integer files are smaller but the marker for ", " missing field value is 999999 while on double ", " precision files it is NaN (= Not a number). ", " The unit of field strength is 0.1 nT in integer ", " files and 1 nT in double precision files. ", " [-p] This may be used to prevent the usage text from ", " appearing when there are no other parameters ", " and data is read from pipeline. ", " IAGAFile Name of IAGA-format file. ", " MatlabFile Name of Matlab binary file. ", }; /*--------------------------------------------------------------------------*/ /* Here is a constant that determines how your computes handles binary */ /* numbers. The constant consists of four parts: thousands, hundreds, tens */ /* and ones. Change this constant to appropriate value for your computer. */ /* Thousands: 0 for PC, Alpha (Little endian machines) */ /* 1 for Mac,SPARC,Apollo,SGI,HP (Big endian machines) */ /* 2 for VAX D-float */ /* 3 for VAX G-float */ /* Hundreds : 0 always */ /* Tens: 0 if double precision numbers, 1 if single precision, */ /* 2 if signed 32-bit numbers, 3 if signed 16-bit numbers. */ /* Ones: 0 for numeric matrices, 1 for text matrix. */ /*--------------------------------------------------------------------------*/ /* #define ThisMachine 1000 */ #define ThisMachine 0 /*--------------------------------------------------------------------------*/ /* Here is the double precision NaN (=Not a Number). */ /*--------------------------------------------------------------------------*/ #define Double_NaN *(double *)NaN /* static char NaN[10] = {0x7F,0xFF,0xFF,0xFF,0xFF,0xFF,0x00,0x00,0x00,0x00}; */ /* Big endian machines */ static char NaN[10] = {0x00,0x00,0x00,0x00,0xFF,0xFF,0xFF,0xFF,0xFF,0x7F}; /* Little endian machines */ /*--------------------------------------------------------------------------*/ /* Here is the header used in storing a single matrix into a Matlab binary */ /* file. For more details see 'Matlab External Interface Guide' p.1-29. */ /*--------------------------------------------------------------------------*/ struct MLmatrix { int type; /* Type flag */ int mrows; /* row dimension */ int ncols; /* column dimension */ int imagf; /* flag indicating imaginary part */ int namlen; /* name length (including NULL) */ }; typedef struct MLmatrix Fmatrix; /*--------------------------------------------------------------------------*/ /* Copy data from integer matrix 'From' into double precision matrix 'To' */ /* and change the field unit into 1 nT (i.e. divide field values by 10). */ /* Also change the missing value markers from 999999 to NaN (Not a Number) */ /*--------------------------------------------------------------------------*/ static void CopyData(long *From,double *To,long Count) { long i; long *p = From; double *d = To; for (i=0;iEndTime - s->StartTime)/s->TimeStep; for (i=0;iStationList; RowCount = ((long)(s->EndTime - s->StartTime))/(s->TimeStep) + 2; Header.type = ThisMachine+20; /* 20 for integer values */ Header.mrows = RowCount; Header.ncols = 3; Header.imagf = 0; Header.namlen = 10; /* 3 for station + 6 for date + '\0' */ /*--- Allocate space for a temporary matrix containing ---*/ /*--- time information (2 rows) + magnetometer data ---*/ DataMatrix = malloc(3*RowCount*VarSize); if (DataMatrix == NULL) { fprintf(stderr,"!!! Memory full error !!!\n"); return; } /*--- Put time information into the first two rows ---*/ /*--- of DataMatrix. This should be same for all stations. ---*/ /*--- Row 1 : Year Month Day ---*/ /*--- Row 2 : Hour Minute SampleRate ---*/ SecsToTm(s->StartTime,&Tm); DataMatrix[0] = 1900+Tm.tm_year; DataMatrix[RowCount] = Tm.tm_mon+1; DataMatrix[2*RowCount] = Tm.tm_mday; DataMatrix[1] = Tm.tm_hour; DataMatrix[RowCount+1] = Tm.tm_min; DataMatrix[2*RowCount+1] = s->TimeStep; /*--- Go through all stations in NETWORK ---*/ s = NETWORK->StationList; while (s != NULL) { /* Set the name of the matrix */ strncpy(Name,s->StationID,3); /* StationID (three letters) */ SecsToStr(s->StartTime,Name+3); /* Starttime YYDDMM */ Name[9] = (char) 0x00; /* Set missing data marker as 999999 */ SetMissing(s,s->X); SetMissing(s,s->Y); SetMissing(s,s->Z); /* Copy data from Station structure to DataMatrix */ memcpy(&DataMatrix[2],s->X,VarSize*(RowCount-2)); memcpy(&DataMatrix[RowCount+2],s->Y,VarSize*(RowCount-2)); memcpy(&DataMatrix[2*RowCount+2],s->Z,VarSize*(RowCount-2)); /* Write matrix into stdout */ fwrite(&Header,sizeof(Fmatrix),1,stdout); fwrite(Name,sizeof(char),Header.namlen,stdout); fwrite(DataMatrix,VarSize,3*RowCount,stdout); s = s->Next; } } /*--------------------------------------------------------------------------*/ /* This is the same procedure as above but writes the field values with */ /* double precision. The missing data values are written as NaN's (= Not a */ /* number). The unit of field strength is 1 nT. */ /*--------------------------------------------------------------------------*/ void WriteDoublePrecMatlab(NetworkPtr NETWORK) { StationPtr s; /* Pointer to current station structure */ double *DataMatrix; /* Pointer to matrix magnetometer data (double) */ long RowCount; /* Number of data points +2 */ Time_struct Tm; /* For getting the starttime in years, months ..*/ Fmatrix Header; /* Header containing matrix info */ char Name[20]; /* Matrix name */ long VarSize = sizeof(double); /*--- First fill out the elements of the header. ---*/ /*--- These should be the same for all stations. ---*/ s = NETWORK->StationList; RowCount = ((long)(s->EndTime - s->StartTime))/(s->TimeStep) + 2; Header.type = ThisMachine; Header.mrows = RowCount; Header.ncols = 3; Header.imagf = 0; Header.namlen = 10; /* 3 for station + 6 for date + '\0' */ /*--- Allocate space for a temporary matrix containing ---*/ /*--- time information (2 rows) + magnetometer data ---*/ DataMatrix = malloc(3*RowCount*VarSize); if (DataMatrix == NULL) { fprintf(stderr,"!!! Memory full error !!!\n"); return; } /*--- Put time information into the first two rows ---*/ /*--- of DataMatrix. This should be same for all stations. ---*/ /*--- Row 1 : Year Month Day ---*/ /*--- Row 2 : Hour Minute SampleRate ---*/ SecsToTm(s->StartTime,&Tm); DataMatrix[0] = 1900+Tm.tm_year; DataMatrix[RowCount] = Tm.tm_mon+1; DataMatrix[2*RowCount] = Tm.tm_mday; DataMatrix[1] = Tm.tm_hour; DataMatrix[RowCount+1] = Tm.tm_min; DataMatrix[2*RowCount+1] = s->TimeStep; /*--- Go through all stations in NETWORK ---*/ s = NETWORK->StationList; while (s != NULL) { /* Set the name of the matrix */ strncpy(Name,s->StationID,3); /* StationID (three letters) */ SecsToStr(s->StartTime,Name+3); /* Starttime YYDDMM */ Name[9] = (char) 0x00; /* Copy the field values into DataMatrix and change the */ /* unit of field strength into 1 nT. Also convert the */ /* the missing data values into NaN's. */ CopyData(s->X,&DataMatrix[2],RowCount-2); CopyData(s->Y,&DataMatrix[RowCount+2],RowCount-2); CopyData(s->Z,&DataMatrix[2*RowCount+2],RowCount-2); /* Write matrix into stdout */ fwrite(&Header,sizeof(Fmatrix),1,stdout); fwrite(Name,sizeof(char),Header.namlen,stdout); fwrite(DataMatrix,VarSize,3*RowCount,stdout); s = s->Next; } } /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; long status = 0; long FileCount = 0; long HourCount = 0; long DoublePrecision = 1; char FileName[100] = ""; char StartTimeStr[14] = ""; char EndTimeStr[14] = ""; char StationStr[200] = ""; Network_struct NET; /*==========================*/ /* 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 'o' : strcpy(StationStr,argv[++params]); break; case 'i' : DoublePrecision = 0; break; case 'p' : /* Do nothing */ 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) { SecsToStr(StrToSecs(StartTimeStr,0)+3600*HourCount,EndTimeStr); } /*====================================================*/ /* Read the data from an IAGA format file into memory */ /*====================================================*/ status = ReadIAGA(FileName,&NET,StartTimeStr,EndTimeStr,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; } /*======================================================*/ /* Write the data in Matlab binary format into stdout */ /*======================================================*/ if (DoublePrecision) WriteDoublePrecMatlab(&NET); else WriteIntMatlab(&NET); FreeNetwork(&NET); return OK; }