xrit2grib.cpp

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

  delete [ ] segsindexes;

  return(0);
}

void usage(char *pname)
{
  std::cout << pname << ": Convert HRIT/LRIT files to Grib format."
	    << std::endl;
  std::cout << std::endl << "Usage:" << std::endl << "\t"
            << pname << " [-d directory] -r resol -s prodid1 -c prodid2 -t time"
            << " [-a linestart,nlin,pixstart,npix]" << std::endl << std::endl
            << "Example: " << std::endl << "\t" << pname
            << " -d data/HRIT -r H -s MSG1 -c HRV -t 200402101315"
            << std::endl;
  return;
}

char *chname(char *chdesc, int len)
{
  char *name = new char[len];
  for (int i = 0; i < len; i ++)
    if (chdesc[i] == ' ' || chdesc[i] == '.') name[i] = '_';
    else name[i] = chdesc[i];
  return name;
}

std::string underscoreit(std::string base, int final_len)
{
  char *tmp;
  int mlen;

  tmp = (char *) malloc((final_len+1) * sizeof(char));
  tmp[final_len] = 0;
  memset(tmp, 0x5F, final_len);
  mlen = strlen(base.c_str( ));
  if (mlen > final_len) mlen = final_len;
  memcpy(tmp, base.c_str( ), mlen);
  std::string retval = tmp;
  free(tmp);
  return retval;
}

//
// Creates Grib product
//
bool GribProduct(MSG_header *PRO_head, MSG_data* PRO_data,
                   int totalsegs, int *segsindexes,
		   MSG_header *header, MSG_data *msgdat)
{
  struct tm *tmtime;
  char GribName[1024];
  int bpp;
  int ncal;
  float *cal;
  GRIB_MSG_PDS pds;

  // std::cout << PRO_data->prologue->sat_status;
  // std::cout << PRO_data->prologue->image_acquisition;
  // std::cout << PRO_data->prologue->celestial_events;
  // std::cout << PRO_data->prologue->image_description;
  // std::cout << PRO_data->prologue->radiometric_proc;
  // std::cout << PRO_data->prologue->geometric_proc;

  int npix = header[0].image_structure->number_of_columns;
  int nlin = header[0].image_structure->number_of_lines;
  size_t npixperseg = npix*nlin;
  nlin *= totalsegs;
  int SP_X0 = npix/2;
  int SP_Y0 = nlin/2;

  size_t total_size = totalsegs*npixperseg;
  MSG_SAMPLE *pixels = new MSG_SAMPLE[total_size];
  MSG_SAMPLE *rotated = new MSG_SAMPLE[total_size];
  memset(pixels, 0, total_size*sizeof(MSG_SAMPLE));
  size_t pos = 0;
  for (int i = 0; i < totalsegs; i ++)
  {
    if (segsindexes[i] >= 0)
      memcpy(pixels+pos, msgdat[segsindexes[i]].image->data,
             npixperseg*sizeof(MSG_SAMPLE));
    pos += npixperseg;
  }

  // Rotate 180deg the image!!!
  for (int i = 0; i < nlin; i ++)
    for (int j = 0; j < npix; j ++)
      rotated[i*npix+j] = pixels[((nlin-1)-i)*npix+(npix-1-j)];

  delete [ ] pixels;
  pixels = rotated;

  // Manage subarea
  if (is_subarea)
  {
    if (AreaLinStart < 0                             ||
        AreaLinStart > nlin - AreaNlin ||
        AreaNlin > nlin - AreaLinStart)
    {
      std::cerr << "Wrong Subarea in lines...." << std::endl;
      throw;
    }
    if (AreaPixStart < 0               ||
        AreaPixStart > npix - AreaNpix ||
        AreaNpix > npix - AreaPixStart)
    {
      std::cerr << "Wrong Subarea in Pixels...." << std::endl;
      throw;
    }
    size_t newsize = AreaNpix * AreaNlin;
    MSG_SAMPLE *newpix = new MSG_SAMPLE[newsize];
    memset(newpix, 0, newsize*sizeof(MSG_SAMPLE));
    for (int i = 0; i < AreaNlin; i ++)
      memcpy(newpix + i * AreaNpix,
             pixels + (AreaLinStart + i) * npix + AreaPixStart,
             AreaNpix * sizeof(MSG_SAMPLE));
    delete [ ] pixels;
    pixels = newpix;
    total_size = newsize;
  }
  else
  {
    AreaNpix = npix;
    AreaNlin = nlin;
  }

  tmtime = PRO_data->prologue->image_acquisition.PlannedAquisitionTime.TrueRepeatCycleStart.get_timestruct( );

  // Fill pds extra section

  pds.spc = (unsigned short) header[0].segment_id->spacecraft_id;
  pds.chn = (unsigned char ) header[0].segment_id->spectral_channel_id;
  pds.sublon = header[0].image_navigation->subsatellite_longitude;
  pds.npix = npix;
  pds.nlin = nlin;
  pds.cfac = header[0].image_navigation->column_scaling_factor;
  pds.lfac = header[0].image_navigation->line_scaling_factor;
  pds.coff = header[0].image_navigation->column_offset;
  pds.loff = header[0].image_navigation->line_offset;
  pds.sh = header[0].image_navigation->satellite_h;
  MSG_data_level_15_header *p = PRO_data->prologue;
  pds.cal_offset = p->radiometric_proc.ImageCalibration[pds.chn-1].Cal_Offset;
  pds.cal_slope = p->radiometric_proc.ImageCalibration[pds.chn-1].Cal_Slope;

  char *channelstring = strdup(MSG_channel_name((t_enum_MSG_spacecraft) pds.spc,
                               pds.chn).c_str( ));
  char *channel = chname(channelstring, strlen(channelstring) + 1);

  // Build up output Grib file name and open it
  sprintf( GribName, "%s_%4d%02d%02d_%02d%02d.grb", channel,
           tmtime->tm_year + 1900, tmtime->tm_mon + 1, tmtime->tm_mday,
	   tmtime->tm_hour, tmtime->tm_min );

  GRIB_FILE gf;
  int ret = gf.OpenWrite(GribName);
  if (ret != 0) return false;

  GRIB_MESSAGE m;
  GRIB_GRID g;
  GRIB_LEVEL l;
  GRIB_TIME t;
  GRIB_FIELD f;

  int fakechan = 0;
  float abase = 145.0;
  if (pds.chn == MSG_SEVIRI_1_5_VIS_0_6) { fakechan = 10; abase = 0.0; }
  if (pds.chn == MSG_SEVIRI_1_5_VIS_0_8) { fakechan = 11; abase = 0.0; }
  if (pds.chn == MSG_SEVIRI_1_5_HRV)     { fakechan = 12; abase = 0.0; }
  if (pds.chn == MSG_SEVIRI_1_5_IR_1_6)   fakechan = 0;
  if (pds.chn == MSG_SEVIRI_1_5_IR_3_9)   fakechan = 1;
  if (pds.chn == MSG_SEVIRI_1_5_IR_8_7)   fakechan = 2;
  if (pds.chn == MSG_SEVIRI_1_5_IR_9_7)   fakechan = 3;
  if (pds.chn == MSG_SEVIRI_1_5_IR_10_8)  fakechan = 4;
  if (pds.chn == MSG_SEVIRI_1_5_IR_12_0)  fakechan = 5;
  if (pds.chn == MSG_SEVIRI_1_5_IR_13_4)  fakechan = 6;
  if (pds.chn == MSG_SEVIRI_1_5_WV_6_2)   fakechan = 20;
  if (pds.chn == MSG_SEVIRI_1_5_WV_7_3)   fakechan = 21;
  fakechan = pds.chn;

  t.set(tmtime->tm_year+1900, tmtime->tm_mon+1,
        tmtime->tm_mday, tmtime->tm_hour, tmtime->tm_min,
        GRIB_TIMEUNIT_MINUTE, GRIB_TIMERANGE_FORECAST_AT_REFTIME_PLUS_P1,
        0, 0, 0, 0);

  // l.set(GRIB_LEVEL_SATELLITE_METEOSAT8, pds.chn, LOCALDEF24FUNC);
 
  l.set(GRIB_LEVEL_SATELLITE_METEOSAT8, fakechan, LOCALDEF3FUNC);

  // Dimensions
  bpp = header[0].image_structure->number_of_bits_per_pixel;
  ncal = (int) pow(2.0, bpp);

  g.set_size(AreaNpix, AreaNlin);
  g.set_earth_spheroid( );
  g.is_dirincgiven = true;
  // Earth Equatorial Radius is 6378.160 Km (IAU 1965)
  g.set_spaceview(0.0, pds.sublon, SEVIRI_DX, SEVIRI_DY,
                  SP_X0, SP_Y0, SEVIRI_ORIENTATION,
                  SEVIRI_CAMERA_H, AreaPixStart+1, AreaLinStart+1);

  // Get calibration values
  cal = PRO_data->prologue->radiometric_proc.get_calibration((int) pds.chn,
                  bpp);

  float *fvals = new float[total_size];
  for (int i = 0; i < (int) total_size; i ++)
  {
    if (pixels[i] > 0) fvals[i] = cal[pixels[i]] - abase;
    else fvals[i] = 0.0;
  }

  f.set_table(GRIB_CENTER, GRIB_SUBCENTER, GRIB_TABLE, GRIB_PROCESS);
  f.set_field(GRIB_PARAMETER_IMG_D, fvals, total_size, FILL_VALUE, FILL_VALUE);
  f.set_scale(GRIB_DECIMAL_SCALE);

  unsigned char localdefinition3[LOCALDEF3LEN];

  localdefinition3[0] = LOCALDEF3ID;
  localdefinition3[1] = LOCALDEF3CLASS;
  localdefinition3[2] = LOCALDEF3TYPE;
  localdefinition3[3] = (LOCALDEF3STREAM >> 8) & 255;
  localdefinition3[4] = LOCALDEF3STREAM & 255;
  localdefinition3[5] = 0x30; // 0
  localdefinition3[6] = 0x30; // 0
  localdefinition3[7] = 0x30; // 0
  localdefinition3[8] = 0x31; // 1
  localdefinition3[9] = fakechan;
  localdefinition3[10] = LOCALDEF3FUNC;
  localdefinition3[11] = 0;

  // unsigned char localdefinition24[LOCALDEF24LEN];

  // localdefinition24[0] = LOCALDEF24ID;
  // localdefinition24[1] = LOCALDEF24CLASS;
  // localdefinition24[2] = LOCALDEF24TYPE;
  // localdefinition24[3] = (LOCALDEF24STREAM >> 8) & 255;
  // localdefinition24[4] = LOCALDEF24STREAM & 255;
  // localdefinition24[5] = 0x30; // 0
  // localdefinition24[6] = 0x30; // 0
  // localdefinition24[7] = 0x30; // 0
  // localdefinition24[8] = 0x31; // 1
  // localdefinition24[9] = (LOCALDEF24SATID >> 8) & 255;
  // localdefinition24[10] = LOCALDEF24SATID & 255;
  // localdefinition24[11] = (LOCALDEF24INSTR >> 8) & 255;
  // localdefinition24[12] = LOCALDEF24INSTR & 255;
  // localdefinition24[13] = (pds.chn >> 8) & 255;
  // localdefinition24[14] = pds.chn & 255;
  // localdefinition24[15] = LOCALDEF24FUNC;

  std::cout << "Adding local use PDS of " << sizeof(GRIB_MSG_PDS) << " bytes"
            << std::endl;
  size_t add_total = LOCALDEF3LEN + sizeof(GRIB_MSG_PDS);
  unsigned char *pdsadd = new unsigned char[add_total];
  memcpy(pdsadd, localdefinition3, LOCALDEF3LEN);
  memcpy(pdsadd+LOCALDEF3LEN, &pds, sizeof(GRIB_MSG_PDS));
  f.add_local_def(add_total, pdsadd);

  // Write output values
  m.set_time(t);
  m.set_grid(g);
  m.set_level(l);
  m.set_field(f);

  gf.WriteMessage(m);

  // Close Grib output
  ret = gf.Close( );
  if (ret != 0) return -1;

  delete [ ] pixels;
  delete [ ] cal;
  delete [ ] fvals;

  return( true );
}

//---------------------------------------------------------------------------