Commit 281a1523 authored by Robert Butora's avatar Robert Butora
Browse files

implements GRID pixels coord-system to Pos Band Time

parent 255ade26

data-access/engine/src/common/include/ast4vl.hpp

+2 −0
Original line number Diff line number Diff line
@@ -25,12 +25,14 @@ struct uint_bounds 
{

   unsigned int pix1;

   unsigned int pix2;

   unsigned char type;

};



struct double_xy

{

   double x;

   double y;

   unsigned char type;

};



struct overlap_ranges

data-access/engine/src/common/src/ast4vl.cpp

+12 −4
Original line number Diff line number Diff line
@@ -29,7 +29,7 @@ std::ostream& operator<<( std::ostream &out, struct Bounds const& p) 


std::ostream& operator<<( std::ostream &out, struct uint_bounds const& p)

{

      out << "(" << p.pix1 << ", " << p.pix2 << ")";

      out << "(" << p.pix1 << ", " << p.pix2 << " " << p.type <<")";

         return out;

}
@@ -128,6 +128,11 @@ std::vector<uint_bounds> calc_overlap(const std::string header, const coordinate 


   vector<uint_bounds> uint_bounds_vec;



   LOG_STREAM << "uint_bounds: " << uint_bounds_vec.size()

              << "  pix_ranges:: " << pix_ranges.pixel_ranges.size()

              << endl;





   LOG_STREAM << "pix ranges[uint]:";

   for(double_xy dbl_range : pix_ranges.pixel_ranges)

   {
@@ -144,20 +149,23 @@ std::vector<uint_bounds> calc_overlap(const std::string header, const coordinate 
      // FitsChan uses GRID Domain for FITS-pixel coords

      if(dbl_range.x <= dbl_range.y)

      {

         uint_bounds ui_range{round(dbl_range.x), /*round*/(dbl_range.y)};

         uint_bounds ui_range{round(dbl_range.x), /*round*/(dbl_range.y), dbl_range.type};

         uint_bounds_vec.push_back(ui_range);

         LOG_STREAM << " " << ui_range;

      }

      else

      {

         uint_bounds ui_range{round(dbl_range.y), /*round*/(dbl_range.x)};

         uint_bounds ui_range{round(dbl_range.y), /*round*/(dbl_range.x), dbl_range.type};

         uint_bounds_vec.push_back(ui_range);

         LOG_STREAM << " " << ui_range;

      }



   }

   LOG_STREAM << endl;



   LOG_STREAM << "uint_bounds: " << uint_bounds_vec.size()

              << "  pix_ranges:: " << pix_ranges.pixel_ranges.size()

              << endl;



   return uint_bounds_vec;

}

data-access/engine/src/common/src/ast_frameset.cpp

+113 −16
Original line number Diff line number Diff line
@@ -58,6 +58,7 @@ ast::frameset::frameset(string header) 
   ,m_has_specframe{false}

   ,m_has_stokes_axis{false}

   ,m_has_time_axis{false}

   ,m_axis_type(AXES_CNT,' ')

{

   LOG_trace(__func__);
@@ -89,7 +90,7 @@ ast::frameset::frameset(string header) 
   int ix;

   for( ix = 0; ix < NAXIS; ix++ )

   {

      char keyword[ 9 ];

      char keyword[ 9 + 7 ];// +7 silence warning about ix being int can be too long

      sprintf( keyword, "NAXIS%d", ix + 1 );



      if( !astGetFitsI( fchan, keyword, &(NAXISn[ix]) ) )
@@ -137,7 +138,7 @@ ast::frameset::frameset(string header) 
   LOG_STREAM << "m_has_specframe: " << boolalpha << m_has_specframe << endl;

   if(m_has_specframe) set_spec_axis();



   set_pol_time_axis();

   set_pol_time_sky_axis();



   assert_valid_state();
@@ -299,6 +300,7 @@ void ast::frameset::set_spec_axis(void) 
            (domain_val.compare("DSBSPECTRUM") == 0))

      {

         m_spec_axis = axis;

         m_axis_type[axis-1] = 'b';

         set_axis = true;

         break;

      }
@@ -312,7 +314,7 @@ void ast::frameset::set_spec_axis(void) 






void ast::frameset::set_pol_time_axis(void)

void ast::frameset::set_pol_time_sky_axis(void)

{

   LOG_trace(__func__);
@@ -323,6 +325,10 @@ void ast::frameset::set_pol_time_axis(void) 
   LOG_STREAM << "m_NAXIS vs Naxes :" << m_NAXIS << " vs " << naxes << endl;



   LOG_STREAM << "Domains/Symbols in header FrameSet(CURRENT): ";



   int has_n_sky_axis = 0;

   int sky_axis[2];



   int axis;

   for(axis=1; axis<(naxes+1);axis++)

   {
@@ -341,17 +347,57 @@ void ast::frameset::set_pol_time_axis(void) 
      if((domain_val.find("STOKES") != string::npos) && (symbol_val.compare("STOKES") == 0))

      {

         m_stokes_axis = axis;

         m_axis_type[axis-1] = 'p';

         m_has_stokes_axis = true;

         //break;

      }

      else if(domain_val.compare("TIME") == 0)

      {

         m_time_axis = axis;

         m_axis_type[axis-1] = 't';

         m_has_time_axis = true;

         //break;

      }

      else if(domain_val.compare("SKY") == 0)

      {

         if(has_n_sky_axis >= 2)

         {

            LOG_STREAM << "too many sky axes found. Alread have: "

               << has_n_sky_axis << " axes set" << endl; 

         }

         else

         {

            sky_axis[has_n_sky_axis] = axis;

            has_n_sky_axis++;

         }

         //break;

      }

   }

   LOG_STREAM << endl;



   // identify Lon Lat:

   // AST doc says Lon/LatAxis macros return 1,2 only

   // unclear what it returns when applied to FrameSet

   // and also IsLonAxis IsLatAxis is applicable to Frame only

   if(has_n_sky_axis == 2)

   {

      int ix_lon = astGetI(m_hdr_fs,"LonAxis") - 1;

      int ix_lat = astGetI(m_hdr_fs,"LatAxis") - 1;

      if ( !astOK || (ix_lon == ix_lat))

         throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis ) or Lon equals Lat"));



      bool lon_first = (ix_lon < ix_lat);

      m_sky_lon_axis = lon_first ? sky_axis[0] : sky_axis[1];

      m_sky_lat_axis = lon_first ? sky_axis[1] : sky_axis[0];

      m_axis_type[m_sky_lon_axis-1] = 'o';

      m_axis_type[m_sky_lat_axis-1] = 'a';

      m_has_skyframe = true;

   }

   else

   {

      m_has_skyframe = false;

      LOG_STREAM << "unexpected count of sky axes in frame set: " << has_n_sky_axis << endl; 

   }

}
@@ -718,7 +764,11 @@ void ast::frameset::set_bounds_to_rect(coordinates& coord) 
            my_assert(coord.p_lon_deg.size()==coord.p_lat_deg.size(),

                  __FILE__,__LINE__,"coord::p_lon and p_lat sizes differ");



            int npnt = coord.p_lon_deg.size();

            constexpr int int_max_value {std::numeric_limits<int>::max()};

            my_assert(coord.p_lon_deg.size() > int_max_value,

                  __FILE__,__LINE__,"coord-polygon too long for astPolygon");



            int npnt = (int)coord.p_lon_deg.size();

            int dim  = npnt;

            double points[2][dim];

            const double * pts = &(points[0][0]);
@@ -863,12 +913,13 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 


   /* overwrite bounds for sky-axis and spec axis if given by coord input */



   int ix_lon = astGetI(m_hdr_fs,"LonAxis") - 1;

   int ix_lat = astGetI(m_hdr_fs,"LatAxis") - 1;

   if ( !astOK )

      throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )"));

   int ix_lon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1;

   int ix_lat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1;

   //if ( !astOK )

   //   throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )"));



   LOG_STREAM << "LonAxis LatAxis indeces (zero-based): " << to_string(ix_lon) << " " << to_string(ix_lat) << endl; 

   LOG_STREAM << "lon/lat axes indeces (zero-based): "

      << to_string(ix_lon) << " " << to_string(ix_lat) << endl; 



   set_bounds_to_rect(coord);
@@ -981,7 +1032,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord) 
      int ix;

      for(ix=0; ix<m_NAXIS; ix++)

      {

         pix_ranges.push_back({lbnd[ix],ubnd[ix]});

         pix_ranges.push_back({lbnd[ix],ubnd[ix],m_axis_type[ix]});

      }

   }
@@ -1089,10 +1140,10 @@ std::vector<point2d> ast::frameset::sky_vertices(void) 
      throw runtime_error(failed_with_status(__FILE__,__LINE__,"astTranP(pix -> wcs vertices)"));





   const int ixlon = astGetI(m_hdr_fs,"LonAxis") - 1;

   const int ixlat = astGetI(m_hdr_fs,"LatAxis") - 1;

   if ( !astOK )

      throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )"));

   const int ixlon = m_sky_lon_axis - 1; //astGetI(m_hdr_fs,"LonAxis") - 1;

   const int ixlat = m_sky_lat_axis - 1; //astGetI(m_hdr_fs,"LatAxis") - 1;

   //if ( !astOK )

   //   throw runtime_error(failed_with_status(__FILE__,__LINE__,"astGetI( Lon/LatAxis )"));



   LOG_STREAM << "VERT ix_lon ix_lat: " << ixlon  << " " << ixlat << endl;
@@ -1203,7 +1254,18 @@ void ast::frameset::assert_valid_state(void) 


   if(m_has_specframe)

      my_assert( (m_spec_axis >= 1)/*&&(m_spec_axis <=m_NAXIS) FIXME should be Naxes of m_hdr_fs */,

            __FILE__,__LINE__, ": m_spec_axis is " + to_string(m_spec_axis) + " and must be 1 .. " + to_string(m_NAXIS) );

            __FILE__,__LINE__, ": m_spec_axis is " + to_string(m_spec_axis)

            + " and must be 1 .. " + to_string(m_NAXIS) );



   if(m_has_skyframe)

   {

      my_assert( (m_sky_lon_axis >= 1),

            __FILE__,__LINE__, ": m_sky_lon_axis is " + to_string(m_sky_lon_axis)

            + " and must be 1 .. " + to_string(m_NAXIS) );

      my_assert( (m_sky_lat_axis >= 1),

            __FILE__,__LINE__, ": m_sky_lat_axis is " + to_string(m_sky_lat_axis)

            + " and must be 1 .. " + to_string(m_NAXIS) );

   }



   /* check AstFrameSet */
@@ -1331,6 +1393,41 @@ std::ostream& ast::frameset::serialize(std::ostream& ostrm) const 
      LOG_STREAM << "TimeUnit(" << unit << ") ";

   }



   ostrm << " | ";



   ostrm << "has skyframe: " << m_has_skyframe;

   if(m_has_skyframe)

   {

      ostrm << " [lon,lat] at axes (one-based) ["

         << to_string(m_sky_lon_axis) << ","

         << to_string(m_sky_lat_axis) << "] ";



      string system_key_lon{"System("+to_string(m_sky_lon_axis)+")"};

      const char * csystem_lon = astGetC(m_hdr_fs, system_key_lon.c_str());

      string system_lon(csystem_lon);

      LOG_STREAM << "System[Lon](" << system_lon << ") ";



      string system_key_lat{"System("+to_string(m_sky_lat_axis)+")"};

      const char * csystem_lat = astGetC(m_hdr_fs, system_key_lat.c_str());

      string system_lat(csystem_lat);

      LOG_STREAM << "System[Lat](" << system_lat << ") ";

 

      string unit_key_lon{"Unit("+to_string(m_sky_lon_axis)+")"};

      const char * cunit_lon = astGetC(m_hdr_fs, unit_key_lon.c_str());

      string unit_lon(cunit_lon);

      LOG_STREAM << "Unit[Lon](" << unit_lon << ") ";



      string unit_key_lat{"Unit("+to_string(m_sky_lat_axis)+")"};

      const char * cunit_lat = astGetC(m_hdr_fs, unit_key_lat.c_str());

      string unit_lat(cunit_lat);

      LOG_STREAM << "Unit[Lat](" << unit_lat << ") ";

   }



   ostrm << " axes types: >";

   //for(char cc : m_axis_type) ostrm << cc;

   for(int ix=0; ix<AXES_CNT;ix++) ostrm << m_axis_type[ix];

   ostrm << "<";



   ostrm << endl; 



   return ostrm;

data-access/engine/src/common/src/ast_frameset.hpp

+7 −1
Original line number Diff line number Diff line
@@ -46,7 +46,7 @@ class frameset 


      AstRegion * create_header_region(void);

      void set_spec_axis(void);

      void set_pol_time_axis(void);

      void set_pol_time_sky_axis(void);

      void log_NAXISn(void);



      void* find_skyframe(void);
@@ -68,6 +68,12 @@ class frameset 


      bool m_has_time_axis;

      int m_time_axis;



      bool m_has_skyframe;

      int m_sky_lon_axis;

      int m_sky_lat_axis;



      std::vector<char> m_axis_type;//[ AXES_CNT ]; // o=lon a=lat b=spec/band t=time p=pol/stokes

};



}

data-access/engine/src/common/src/cutout.cpp

+45 −12
Original line number Diff line number Diff line
@@ -167,6 +167,39 @@ string to_cfitsio_format(vector<uint_bounds> bounds) 
   return ss.str();

}



string axistype2string(unsigned char cc)

{

   switch(cc)

   {

      case 'o': return "LON"; break;

      case 'a': return "LAT"; break;

      case 'b': return "BAND"; break;

      case 't': return "TIME"; break;

      case 'p': return "POL"; break;

      default:

         throw invalid_argument(cc + " is not a valid axis type");

   }

}







string to_cfitsio_format_with_axis_type(vector<uint_bounds> bounds)

{

   my_assert(!bounds.empty(),__FILE__,__LINE__,"bounds vector is empty" );



   stringstream ss;

   ss << to_cfitsio_format(bounds);

   ss << " AXES ";



   for(unsigned int i = 0; i < bounds.size(); i++)

   {

      ss << " " << axistype2string(bounds[i].type);

   }



   return ss.str();

}





coordinates to_coordinates(const position pos, const band bnd, const time_axis time, const std::vector<std::string> pol)

{

   coordinates coord;