diff --git a/data-access/engine/src/common/include/ast4vl.hpp b/data-access/engine/src/common/include/ast4vl.hpp index 46571ab47734b84cbd2969536bd8d14f8c9671de..95d714e3d9089127c2edc0fec2fccf5d0307656d 100644 --- a/data-access/engine/src/common/include/ast4vl.hpp +++ b/data-access/engine/src/common/include/ast4vl.hpp @@ -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 diff --git a/data-access/engine/src/common/src/ast4vl.cpp b/data-access/engine/src/common/src/ast4vl.cpp index a0da6e352297aaa6960d3b101dcd9c3525e6d0b7..c43683500b180186b555288f6329e96fa3bfdf92 100644 --- a/data-access/engine/src/common/src/ast4vl.cpp +++ b/data-access/engine/src/common/src/ast4vl.cpp @@ -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; } diff --git a/data-access/engine/src/common/src/ast_frameset.cpp b/data-access/engine/src/common/src/ast_frameset.cpp index f08299b46fdaa1dfc3ee3fd060c27d7117f826ae..672156126eee9e347a1f5b7dc807affe99819039 100644 --- a/data-access/engine/src/common/src/ast_frameset.cpp +++ b/data-access/engine/src/common/src/ast_frameset.cpp @@ -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); @@ -896,7 +947,7 @@ overlap_ranges ast::frameset::overlap(coordinates coord) /* FIXME ast-overlap computation breaks if bands do not overlap - suspect AST bug ? */ if((( up[m_spec_axis-1] < coord.vl_kmps) && ( up[m_spec_axis-1] < coord.vu_kmps)) || - ((low[m_spec_axis-1] > coord.vl_kmps) && (low[m_spec_axis-1] > coord.vu_kmps))) + ((low[m_spec_axis-1] > coord.vl_kmps) && (low[m_spec_axis-1] > coord.vu_kmps))) { /* set values only to get correct debug print for client */ low[m_spec_axis-1] = coord.vl_kmps; @@ -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; diff --git a/data-access/engine/src/common/src/ast_frameset.hpp b/data-access/engine/src/common/src/ast_frameset.hpp index 23caabec215bba79f5ee0847d7f4847016635de0..413c926a1e17f4a66220de1c9acdd41682ad8650 100644 --- a/data-access/engine/src/common/src/ast_frameset.hpp +++ b/data-access/engine/src/common/src/ast_frameset.hpp @@ -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 }; } diff --git a/data-access/engine/src/common/src/cutout.cpp b/data-access/engine/src/common/src/cutout.cpp index 3651f26f1571b1fb3bc62a4c039f4a1a4b4a94e9..b1875b3a46484cdc3615b8737c4d2e51fc0aea2c 100644 --- a/data-access/engine/src/common/src/cutout.cpp +++ b/data-access/engine/src/common/src/cutout.cpp @@ -147,29 +147,62 @@ int to_v_type(enum specsystem specsys) int to_v_valid(enum specsystem specsys) { - return specsys == specsystem::NONE ? 0 : 1; + return specsys == specsystem::NONE ? 0 : 1; } string to_cfitsio_format(vector<uint_bounds> bounds) { - my_assert(!bounds.empty(),__FILE__,__LINE__,"bounds vector is empty" ); + my_assert(!bounds.empty(),__FILE__,__LINE__,"bounds vector is empty" ); - stringstream ss; - ss << "["; - ss << bounds[0].pix1 << ":" << bounds[0].pix2; - for(unsigned int i = 1; i < bounds.size(); i++) - { - ss << " " << bounds[i].pix1 << ":" << bounds[i].pix2; - } - ss << "]"; - return ss.str(); + stringstream ss; + ss <<"["; + ss << bounds[0].pix1 << ":" << bounds[0].pix2; + for(unsigned int i = 1; i < bounds.size(); i++) + { + ss << " " << bounds[i].pix1 << ":" << bounds[i].pix2; + } + ss << "]"; + 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; + coordinates coord; coord.skysys = pos.sys; if(pos.shape == area::RANGE) diff --git a/data-access/engine/src/vlkb/src/main.cpp b/data-access/engine/src/vlkb/src/main.cpp index af3e0f05f5831c7bd5abc13d62d406a1790d7e88..52b41da1b1216d0016a889de604376f31ea3c38b 100644 --- a/data-access/engine/src/vlkb/src/main.cpp +++ b/data-access/engine/src/vlkb/src/main.cpp @@ -384,7 +384,7 @@ int cmd_overlap(int argc, char * argv[]) if((ov_code >= 2) && (ov_code <= 5)) { rc = EXIT_SUCCESS; - cout << to_cfitsio_format(bnds) << endl; + cout << to_cfitsio_format_with_axis_type(bnds) << endl; } else if((ov_code == 1) || (ov_code == 6)) { diff --git a/data-access/engine/src/vlkb/src/service_string.hpp b/data-access/engine/src/vlkb/src/service_string.hpp index 416dc21c422a0c1f9060585e77be078aa85c05d4..2ce27c9e90b1daba05c13236b3124b673c75ba2e 100644 --- a/data-access/engine/src/vlkb/src/service_string.hpp +++ b/data-access/engine/src/vlkb/src/service_string.hpp @@ -11,5 +11,6 @@ skysystem to_skysystem(std::string str); specsystem to_specsystem(std::string str); specsystem to_specsystem(int i); std::string to_cfitsio_format(std::vector<uint_bounds> bounds); +std::string to_cfitsio_format_with_axis_type(std::vector<uint_bounds> bounds); #endif diff --git a/data-access/servlet/src/main/java/cutout/CutoutImpl.java b/data-access/servlet/src/main/java/cutout/CutoutImpl.java index f22075e057d339e90f0281a1cb2aa0e5cdd96738..dbcb5e3bbb5e13b32c7870b149b7497c7ddebff2 100644 --- a/data-access/servlet/src/main/java/cutout/CutoutImpl.java +++ b/data-access/servlet/src/main/java/cutout/CutoutImpl.java @@ -112,7 +112,6 @@ class CutoutImpl implements Cutout if(bos == null) throw new AssertionError("byte output stream for bounds was not created, is null"); - //String coordString = genRegionForVlkbOverlapCmd(pos, band); JsonEncoder jReq = new JsonEncoder(); jReq.add(pos); jReq.add(band); @@ -139,10 +138,8 @@ class CutoutImpl implements Cutout { boundsString = new String(bos.toByteArray()); - // remove end-of-line (was added by vlkb_ast.cpp: cout << ... << endl) - String lineSeparator = System.lineSeparator(); - boundsString = boundsString.replace(lineSeparator, ""); - LOGGER.info("BOUNDS: " + boundsString); + boundsString = replaceWithGrid(boundsString, pos, band, time, pol); + LOGGER.info("boundsString(with GRID): " + boundsString); has_overlap = !((boundsString != null) && boundsString.trim().isEmpty()); @@ -199,6 +196,78 @@ class CutoutImpl implements Cutout + private String replaceWithGrid(String wcsBounds, Pos pos, Band band, Time time, Pol pol) + { + // remove end-of-line (was added by vlkb_ast.cpp: cout << ... << endl) + String lineSeparator = System.lineSeparator(); + wcsBounds = wcsBounds.replace(lineSeparator, ""); + LOGGER.info("BOUNDS: " + wcsBounds); + + // replace in wcsBounds those bounds where pos,band,time or pol has system=GRID + + String[] substr = wcsBounds.split("(?=AXES)", 2); + for(String ss : substr) LOGGER.info("boundsResult: " + ss); + + String boundsString = substr[0]; + boolean noOverlap = ((boundsString != null) && boundsString.trim().isEmpty()); + if(noOverlap) + { + boundsString = ""; // no overlap + } + else + { + String axesTypes = ""; + if(substr.length > 1) + { + axesTypes = substr[1].replace("AXES"," ").trim(); + LOGGER.info("AXES TYPES: " + axesTypes); + + String[] bnds = normalize(boundsString); + String[] types = normalize(axesTypes); + // assert: bnds.length == types.length + LOGGER.info("boundsCount: " + bnds.length + " typesCount: " + types.length); + + if(bnds.length == types.length) + boundsString = replaceBounds(bnds, types, pos, band); + } + } + return boundsString; + } + + private String replaceBounds(String[] bnds, String[] types, Pos pos, Band band) + { + int ix; + for(ix=0; ix<bnds.length; ix++) + { + if(types[ix].equals("LON") && (pos.system == Pos.System.GRID)) + { + bnds[ix] = pos.lonBoundsString(); + } + else if(types[ix].equals("LAT") && (pos.system == Pos.System.GRID)) + { + bnds[ix] = pos.latBoundsString(); + } + else if(types[ix].equals("BAND") && (band.system == Band.System.GRID)) + { + bnds[ix] = band.boundsString(); + } + } + + LOGGER.info("replaced: " + String.join(" ", bnds)) ; + + return "[" + String.join(" ", bnds) + "]"; + } + + // MAKE SURE vlkb overlap returns space delimited bounds: [ a:b c:d ] + // normalize: strip [,] if any, and split into array by space + private String[] normalize(String spaceStr) + { + String other = spaceStr.replace("[","").replace("]",""); + LOGGER.info("normalize: " + other); + return other.split("\\s+"); + } + + private NullValueCount doCountNullValues(String absPathname, int hdunum) throws IOException, InterruptedException { diff --git a/java-libs/lib/vlkb-volib-0.9.0.jar b/java-libs/lib/vlkb-volib-0.9.0.jar deleted file mode 100644 index 7a469173ef91db8367fb58af748e1f668e405fe3..0000000000000000000000000000000000000000 Binary files a/java-libs/lib/vlkb-volib-0.9.0.jar and /dev/null differ diff --git a/java-libs/lib/vlkb-volib-0.9.1-SNAPSHOT.jar b/java-libs/lib/vlkb-volib-0.9.1-SNAPSHOT.jar new file mode 100644 index 0000000000000000000000000000000000000000..554db77940832d36d03f0cad3addef64dfa3fdc8 Binary files /dev/null and b/java-libs/lib/vlkb-volib-0.9.1-SNAPSHOT.jar differ