Sep 5, 2010

Use GDL to run IDL code and parse GOES AOD binary file

IDL (Interactive Data Language) has a free replacement, named GDL GNU Data Language.
NASA GASP GOES satellite provides AOD data with .aod that needs IDL to parse. Its idl_read code gives some hint about the data format of that binary file.

It is still painful to parse with c code. So GDL is used to function for IDL, which requires non-free licence.
To instal GDL, download it. Download GSL, PLPLOT, CMAKE, etc. packages.
Install all required packages.

After that, download CMSVLIB, and put untared files into <gdl>/share/gnudatalanguage/lib (This is for running 'save' in gdl)

add GDL's bin path to the .bashrc
Possibly, need to launch the program by running:
LD_LIBRARY_PATH=../../plplot/lib:$LD_LIBRARY_PATH gdl
which provides the plplot 's lib path

Besides, for the .aod binary file, its format is guessed to be 10 bytearries of size 2128 x 880. I wrote a c function to parse it also, which is a translation of IDL reader from Appendix A in http://www.ssd.noaa.gov/PS/FIRE/GASP/20090107_GASP_Algorithm_Updates.doc
The c function is also a demo, which is not fully optimized and prettified.

But notice that the archived data may be of size 17000000 bytes, which are before they switched GOES-EAST from GOES-12 to GOES-13. In that case, the array dimension is 2000x850. While there are other associated lat/lon.dat files.

The key points in this function are:
1. there are 10 byte_arrays
2. each of size 2128 x 880 = 1872640 bytes
3. the variable can be declared as uint8_t
4. Matlab can also do this, by
    (1)>>fid=fopen('2010244171519_i18_US.all.aod')
    (2)>>baod=fread(fid, [2128,880],'uint8'); # similar for the rest 9 arrays
5. in c, each array has elements stored as column major. So if loop over with one index from 0 to 1872640, each column is looped first.

#include
#include
#include

int main()
{
    int i = 0;
    FILE *fp;
    long fsize;
    uint8_t *baod;
    uint8_t *bmsk;
    uint8_t *bcls;
    uint8_t *baodstd;
    uint8_t *bsfc;
    uint8_t *bch1;
    uint8_t *bmos;
    uint8_t *bcld;
    uint8_t *bsig;
    uint8_t *bsca;

    fp = fopen("2010244171519_i18_US.all.aod", "rb");
    //2010246174517_i18_US.all.aod","rb");
    //2010244171519_i18_US.all.aod", "rb");

    fseek(fp, 0, SEEK_END);
    fsize = ftell(fp);
    rewind(fp);

    long arrsize = fsize / 10;

    baod = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bmsk = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bcls = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    baodstd = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bsfc = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bch1 = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bmos = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bcld = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bsig = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );
    bsca = (uint8_t *)malloc(sizeof(uint8_t) * arrsize );

    fread(baod, 1, arrsize, fp);
    fread(bmsk, 1, arrsize, fp);
    fread(bcls, 1, arrsize, fp);
    fread(baodstd, 1, arrsize, fp);
    fread(bsfc, 1, arrsize, fp);
    fread(bch1, 1, arrsize, fp);
    fread(bmos, 1, arrsize, fp);
    fread(bcld, 1, arrsize, fp);
    fread(bsig, 1, arrsize, fp);
    fread(bsca, 1, arrsize, fp);

    fclose(fp);

   float *faod = (float *)malloc(sizeof(float) * arrsize );
    float *faodstd = (float *)malloc(sizeof(float) * arrsize );
    float *fsfc = (float *)malloc(sizeof(float) * arrsize );
    float *fch1 = (float *)malloc(sizeof(float) * arrsize );
    float *fmos = (float *)malloc(sizeof(float) * arrsize );
    float *fsig = (float *)malloc(sizeof(float) * arrsize );
    //float *fsca = (float *)malloc(sizeof(float) * arrsize );

    for (i=0; i
    {
        faod[i] = (float)baod[i] / 100.0 - 0.5;
        faodstd[i] = (float)baodstd[i] / 100.0;
        fsfc[i] = (float)bsfc[i] / 500.0 - 0.1;
        fch1[i] = (float)bch1[i] / 600.0;
        fmos[i] = (float)bmos[i] / 600.0;
        fsig[i] = (float)bsig[i] / 250.0 - 0.5;
        //fsca[i] = (float)bsca[i] / 1.0;
    }
   fp = fopen("test0.txt", "w");
   for (i=0; i
    {
        fprintf(fp, "%d %d %d %d %d %d %d %d %d %d %d\n",i,baod[i],bmsk[i],bcls[i],baodstd[i],bsfc[i],bch1[i],bmos[i],bcld[i],bsig[i],bsca[i]);
    }
    fclose(fp);

    for(i=0; i
    {
        if((faodstd[i]>=0.3) || (fsig[i]<=0.01) || (fsfc[i]>=0.15) ||
            (fsfc[i]<=0.005) || (bcls[i]<=15) || (faod[i]>=10.0) ||
            (fch1[i]<=0.0) || (bcld[i]!=1) || (bsca[i]<=70) ||
            (bsca[i]>=170) )
        {
            faod[i] = -9999.0;
        }
    }

    fp = fopen("test.txt","w");
    for (i=0; i
    {
        fprintf(fp, "%d %f\n", i, faod[i]);
    }

    fclose(fp);
   free(baod);
   //free(bmsk);
   free(baodstd);
   free(bsfc);
   free(bch1);
   free(bmos);
   free(bsig);
   //free(bsca);

    free(bcls); free(bsca);
    free(bcld); free(bmsk);

    free(faod);
    free(faodstd);
    free(fsfc);
    free(fch1);
    free(fmos);
    free(fsig);
    //free(fsca);
    return 0;
}

No comments: