xrit2grib.cpp

HRIT读取,用于在LINUX下显示高束数据图像

//---------------------------------------------------------------------------
//
//  File        :   XRIT2Grib.cpp
//  Description :   Export MSG HRIT format in Grib format
//  Project     :   Lamma 2004
//  Author      :   Graziano Giuliani (Lamma Regione Toscana)
//  Source      :   n/a
//
//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//---------------------------------------------------------------------------

#include <config.h>

#include <cstdio>
#include <cmath>
#include <cstdlib>

// Grib Library

#include <GRIB.h>

#define FILL_VALUE         9.999E-20
#define GRIB_DECIMAL_SCALE 1
#define GRIB_CENTER        201 // HIMET
#define GRIB_SUBCENTER     0
#define GRIB_TABLE         1
#define GRIB_PROCESS       254

#define LOCALDEF3LEN       12
#define LOCALDEF3ID        3    // Satellite image data
#define LOCALDEF3CLASS     1    // Operations
#define LOCALDEF3TYPE      40   // Image Data
#define LOCALDEF3STREAM    55   // Meteosat 8 = 55
#define LOCALDEF3FUNC      1    // T in K - 145

#define LOCALDEF24LEN      16
#define LOCALDEF24ID       24   // Satellite image simulation (?)
#define LOCALDEF24CLASS    1    // Operations
#define LOCALDEF24TYPE     40   // Image Data
#define LOCALDEF24STREAM   55   // Meteosat 8 = 55
#define LOCALDEF24SATID    55   // Meteosat 8 = 55
#define LOCALDEF24INSTR    207  // Seviri = 207
#define LOCALDEF24FUNC     1    // T in K - 145

#define SEVIRI_SCANSTEP    (0.00025153/3.0) // scan step angle (rd) for 3 lines
#define SEVIRI_ORIENTATION 180.0
#define SEVIRI_CAMERA_H    6610.6839590101218
#define SEVIRI_DX          2.0*asin(1.0/SEVIRI_CAMERA_H)/SEVIRI_SCANSTEP*1000
#define SEVIRI_DY          SEVIRI_DX

// HRI format interface

#include <hrit/MSG_HRIT.h>

#include <getopt.h>
#include <glob.h>

#define PATH_SEPARATOR "/"
// For windows use #define PATH_SEPARATOR "\"

void usage(char *pname);
std::string underscoreit(std::string base, int final_len);
char *chname(char *chdesc, int len);
bool GribProduct(MSG_header *PRO_head, MSG_data* PRO_data,
                   int totalsegs, int *segsindexes,
		   MSG_header *header, MSG_data *msgdat);

bool is_subarea = false;
int AreaLinStart = 1, AreaNlin, AreaPixStart = 1, AreaNpix;

/* ************************************************************************* */
/* Reads Data Images, performs calibration and store output in Grib format */
/* ************************************************************************* */

int main( int argc, char* argv[] )
{
  char *directory = NULL;
  char *resolution = NULL;
  char *productid1 = NULL;
  char *productid2 = NULL;
  char *timing = NULL;
  char *areastring = NULL;

  char *pname = strdup(argv[0]);

  while (1) {
    int option_index = 0;
    int c;

    static struct option long_options[] = {
      {"directory", 1, 0, 'd'},
      {"resolution", 1, 0, 'r'},
      {"productid1", 1, 0, 's'},
      {"productid2", 1, 0, 'c'},
      {"time", 1, 0, 't'},
      {"subarea", 1, 0, 'a'},
      {"version", 0, 0, 'V'},
      {"help", 0, 0, 'h'},
      {0, 0, 0, 0}
    };

    c = getopt_long (argc, argv, "d:r:s:c:t:a:Vh?",
                     long_options, &option_index);
    if (c == -1)
      break;

    switch (c)
    {
      case 'd':
       directory = strdup(optarg);
       break;
      case 'r':
       resolution = strdup(optarg);
       break;
      case 's':
       productid1 = strdup(optarg);
       break;
      case 'c':
       productid2 = strdup(optarg);
       break;
      case 't':
       timing = strdup(optarg);
       break;
      case 'a':
       is_subarea = true;
       areastring = strdup(optarg);
       sscanf(areastring, "%d,%d,%d,%d",
              &AreaLinStart, &AreaNlin, &AreaPixStart, &AreaNpix);
       AreaPixStart --;
       AreaLinStart --;
       break;
      case 'V':
       std::cout << pname << " " << PACKAGE_STRING << std::endl;
       return 0;
      case '?':
      case 'h':
       usage(pname);
       return(0);
       break;
      default:
       std::cerr << "?? getopt returned character code"
                 << std::oct << c << "??" << std::endl;
       usage(pname);
       return(1);
    }
  }

  if (resolution == NULL || productid1 == NULL ||
      productid2 == NULL || timing == NULL)
  {
    usage(pname);
    return(1);
  }

  std::string filename;

  MSG_header PRO_head;
  MSG_data PRO_data;

  if (directory) filename = directory;
  else filename = ".";

  filename = filename + PATH_SEPARATOR + resolution;
  filename = filename + "-???-??????-";
  filename = filename + underscoreit(productid1, 12) + "-";
  filename = filename + underscoreit("_", 9) + "-";
  filename = filename + "PRO______-";
  filename = filename + timing + "-__";

  glob_t globbuf;
  globbuf.gl_offs = 1;

  if ((glob(filename.c_str( ), GLOB_DOOFFS, NULL, &globbuf)) != 0)
  {
    std::cerr << "No such file(s)." << std::endl;
    return 1;
  }

  if (globbuf.gl_pathc > 1)
  {
    std::cerr << "Non univoque prologue file.... Do not trust calibration."
              << std::endl;
    return 1;
  }

  std::ifstream hrit(globbuf.gl_pathv[1], (std::ios::binary | std::ios::in));
  if (hrit.fail())
  {
    std::cerr << "Cannot open input hrit file "
              << globbuf.gl_pathv[1] << std::endl;
      return 1;
  }
  PRO_head.read_from(hrit);
  PRO_data.read_from(hrit, PRO_head);
  hrit.close( );
  std::cout << PRO_head;

  globfree(&globbuf);

  if (directory) filename = directory;
  else filename = ".";

  filename = filename + PATH_SEPARATOR + resolution;
  filename = filename + "-???-??????-";
  filename = filename + underscoreit(productid1, 12) + "-";
  filename = filename + underscoreit(productid2, 9) + "-";
  filename = filename + "0?????___" + "-";
  filename = filename + timing + "-" + "?_";

  globbuf.gl_offs = 1;

  if ((glob(filename.c_str( ), GLOB_DOOFFS, NULL, &globbuf)) != 0)
  {
    std::cerr << "No such file(s)." << std::endl;
    return 1;
  }

  int nsegments = globbuf.gl_pathc;

  MSG_header *header;
  MSG_data *msgdat;

  header = new MSG_header[nsegments];
  msgdat = new MSG_data[nsegments];

  for (int i = 0; i < nsegments; i ++)
  {
    std::ifstream hrit(globbuf.gl_pathv[i+1],
                       (std::ios::binary | std::ios::in));
    if (hrit.fail())
    {
      std::cerr << "Cannot open input hrit file "
                << globbuf.gl_pathv[i+1] << std::endl;
      return 1;
    }
    header[i].read_from(hrit);
    msgdat[i].read_from(hrit, header[i]);
    hrit.close( );
    // std::cout << header[i];
  }

  globfree(&globbuf);

  if (header[0].segment_id->data_field_format == MSG_NO_FORMAT)
  {
    std::cout << "Product dumped in binary format." << std::endl;
    return 0;
  }

  int totalsegs = header[0].segment_id->planned_end_segment_sequence_number;
  int *segsindexes = new int[totalsegs]; 

  for (int i = 0; i < totalsegs; i ++)
    segsindexes[i] = -1;

  for (int i = 0; i < nsegments; i ++)
    segsindexes[header[i].segment_id->sequence_number-1] = i;

  std::cout << "Writing Grib output..." << std::endl;
  if(!GribProduct(&PRO_head, &PRO_data,
                    totalsegs, segsindexes, header, msgdat))
    throw "Created Grib product NOT OK";

  delete [ ] header;
  delete [ ] msgdat;