Skip to content

Commit

Permalink
smurf: Handle HPX discontinuity at RA=12h
Browse files Browse the repository at this point in the history
  • Loading branch information
David Berry committed Oct 2, 2014
1 parent 56bf8fa commit 7c1fbc0
Show file tree
Hide file tree
Showing 22 changed files with 981 additions and 446 deletions.
3 changes: 3 additions & 0 deletions applications/smurf/libsmf/Makefile.am
View file
Expand Up @@ -281,6 +281,9 @@ smf_isfft.c \
smf_iteratemap.c \
smf_jsainstrument.c \
smf_jsadicer.c \
smf_jsaproj.c \
smf_jsaproj_fromstr.c \
smf_jsaproj_tostr.c \
smf_jsatileheader.c \
smf_jsatile.c \
smf_jsatiling.c \
Expand Down
37 changes: 27 additions & 10 deletions applications/smurf/libsmf/jsatiles.h
View file
Expand Up @@ -15,6 +15,19 @@ typedef enum {
} smf_inst_t;


/* The four HPX-related projections available for mosaics of JSA
tiles. The projection for a mosaic is chosen in order to avoid any
discontinuities being present within the mosaic. Individual JSA tiles
always use SMF__JSA_HPX. */
typedef enum {
SMF__JSA_NULL, /* Not an JSA projection */
SMF__JSA_HPX, /* HPX projection centred on RA=0, Dec=0 */
SMF__JSA_HPX12, /* HPX projection centred on RA=12h, Dec=0 */
SMF__JSA_XPHN, /* XPH projection centred on RA=0, Dec=90 */
SMF__JSA_XPHS /* XPH projection centred on RA=0, Dec=-90 */
} smf_jsaproj_t;


/* Struct to store information about the JSA tiling scheme for a given
JCMT instrument. */

Expand All @@ -37,14 +50,15 @@ void smf_jsainstrument( const char *param, AstFitsChan *fc,
smf_inst_t def, smfJSATiling *tiling,
int *status );
void smf_jsadicer( int indf, const char *base, int trim,
smf_inst_t instrument, int usexph, size_t *ntile,
Grp *grp, int *status );
void smf_jsatile( int itile, smfJSATiling *jsatiling,
int local_origin, int usexph, AstFitsChan **fc,
smf_inst_t instrument, smf_jsaproj_t proj,
size_t *ntile, Grp *grp, int *status );
void smf_jsatile( int itile, smfJSATiling *jsatiling, int local_origin,
smf_jsaproj_t proj, AstFitsChan **fc,
AstFrameSet **fs, AstRegion **region, int lbnd[2],
int ubnd[2], int *status );
AstFitsChan *smf_jsatileheader( int itile, smfJSATiling *jsatiling,
int local_origin, int usexph, int *move, int *status );
AstFitsChan *smf_jsatileheader( int itile, smfJSATiling *skytiling,
int local_origin, smf_jsaproj_t proj,
int *move, int *status );
void smf_jsatilei2xy( int itile, smfJSATiling *jsatiling, int *xt,
int *yt, int *fi, int *status );
int * smf_jsatiles_region( AstRegion *region, smfJSATiling *tiling,
Expand All @@ -53,9 +67,12 @@ int * smf_jsatiles_data( ThrWorkForce *wf, Grp *igrp, size_t size,
smfJSATiling *tiling, int *ntile, int *status );
int smf_jsatilexy2i( int xt, int yt, smfJSATiling *jsatiling,
int *status );
void smf_jsatilexyconv( smfJSATiling *skytiling, int usexph, int xt_hpx,
int yt_hpx, int *xt_out, int *yt_out, int *status );
void smf_jsatilexyconv( smfJSATiling *skytiling, smf_jsaproj_t proj,
int xt_hpx, int yt_hpx, int raw, int *xt_out,
int *yt_out, int *status );
void smf_jsatiling( smf_inst_t instrument, smfJSATiling *jsatiling,
int *status );


smf_jsaproj_t smf_jsaproj( int ntile, const int *tiles,
smfJSATiling *skytiling, int *status );
smf_jsaproj_t smf_jsaproj_fromstr( const char *str, int rep, int *status );
const char * smf_jsaproj_tostr( smf_jsaproj_t proj, int *status );
127 changes: 40 additions & 87 deletions applications/smurf/libsmf/smf_getrefwcs.c
View file
Expand Up @@ -20,10 +20,10 @@
* param = const char * (Given)
* The name of the ADAM parameter used to get the reference NDF. A
* value of "JSA" may also be provided for the parameter, in which
* case the WCS defined by the JSA all-sky pixel grid is returned.
* The choice of JSA projection (HPX, XPHS or XPHN) is based on
* whether the map crosses an HPX discontinuity. Alternatively,
* "JSA-HPX", "JSA-XPHS" or "JSA-XPHN" can be supplied if it is
* case the WCS defined by the SMF__JSA_HPX all-sky projection is
* returned. The suitable JSA projection is chosen to ensure that
* the map does not cross any discontinuities. Alternatively,
* "HPX", "HPX12", "XPHS" or "XPHN" can be supplied if it is
* necessary to specifiy the JSA projection explicitly.
* igrp = Grp * (Given)
* A group holding the paths to the input science files. The first
Expand Down Expand Up @@ -91,10 +91,12 @@
* 26-AUG-2014 (DSB):
* Allow some tolerance for rounding errors when checking whether the
* observation instersects an HPX discontinuity.
* 1-OCT-2014 (DSB):
* Added support for HPX projections centred on RA=12h.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 2007,2013 Science & Technology Facilities Council.
* Copyright (C) 2007,2013-2014 Science & Technology Facilities Council.
* All Rights Reserved.
* Licence:
Expand Down Expand Up @@ -146,22 +148,20 @@ void smf_getrefwcs( const char *param, Grp *igrp, AstFrameSet **specwcs,
AstMapping *splitmap = NULL; /* Mapping from required wcs to PIXEL coords */
AstRegion *circle;
char text[ 255 ]; /* Parameter value */
double centre[ 2 ];
double hpxlat[] = {0.0, 90.0, 180.0, 270.0 };
double p0[ 2 ];
double radius;
int *tiles;
int i;
int inax[ 2 ]; /* Indices of required WCS axes */
int jsatiles;
int lbnd[2]; /* Lower pixel index bounds of mid tile */
int ntile;
int outax[ 7 ]; /* Indices of corresponding PIXEL axes */
int refndf; /* NDF identifier for the refence NDF */
int ubnd[2]; /* Upper pixel index bounds of mid tile */
int usexph;
size_t code;
smfData *data = NULL; /* Structure describing 1st input file */
smfJSATiling skytiling;
smf_inst_t inst = SMF__INST_NONE;
smf_jsaproj_t proj; /* Specific JSA projection to use */
smf_subinst_t subinst;

/* Initialise the returned values. */
Expand Down Expand Up @@ -193,9 +193,11 @@ void smf_getrefwcs( const char *param, Grp *igrp, AstFrameSet **specwcs,
if( *status == PAR__NULL ) {
errAnnul( status );

/* If it is "JSA", we return WCS that describes the JSA all-sky pixel grid. */
/* If it is "JSA", or one of the JSA projection codes, we return WCS that
describes one of the the JSA all-sky pixel grids. */
} else if( *status == SAI__OK ) {
if( astChrMatchN( text, "JSA", 3 ) ) {
proj = smf_jsaproj_fromstr( text, 0, status );
if( astChrMatch( text, "JSA" ) || proj != SMF__JSA_NULL ) {
*isjsa = 1;

/* Report an error if the instrument cannot be determined. */
Expand Down Expand Up @@ -239,9 +241,6 @@ void smf_getrefwcs( const char *param, Grp *igrp, AstFrameSet **specwcs,
instrument. */
smf_jsatiling( inst, &skytiling, status );

/* Assume the HPX projection is to be used. */
usexph = 0;

/* For "JSA" - choose the best projection. */
if( astChrMatch( text, "JSA" ) ) {

Expand All @@ -253,88 +252,42 @@ void smf_getrefwcs( const char *param, Grp *igrp, AstFrameSet **specwcs,
/* Convert the circle to ICRS (as used by the JSA all-sky grid). */
astSetC( circle, "System", "ICRS" );

/* Get the centre and radius of this circle. */
astCirclePars( circle, centre, &radius, NULL );
if( ( centre[ 0 ] == AST__BAD || centre[ 1 ] == AST__BAD ||
radius == AST__BAD ) && *status == SAI__OK ) {
*status = SAI__ERROR;
errRep( "", "smf_getrefwcs: the centre and/or radius of "
"the map is undefined.", status );
}

/* See if the observation may straddle the gap between two polar cusps of
the standard JSA HPX-projected all-sky grid. If it does, we need to
swap to use an XPH (polar-HEALPix) projection instead since otherwise
the map will probably be too big (in pixels) to create. The gaps are
defined by the sections of meridians at longitude=0, 90, 180, 270 that
have absolute latitude greater than 41.8103 degrees (the latitude of
the transition between the equatorial and polar regions of the HPX
projection). If the offset in longitude between the circle centre and
any of these meridians is less than the radius, then we may have an
intersection. Check each meridian in turn. */
usexph = 0;
for( i = 0; i < 4; i++ ) {

/* Get the coords of the point on the meridian which is at the same
latitude as the circle centre. */
p0[ 0 ] = AST__DD2R*hpxlat[ i ];
p0[ 1 ] = centre[ 1 ];

/* If the arc-distance between this point and the circle centre is less
than the radius, the circle intersects the meridian. Use twice the
radius to allow some tolerance (there's no real harm in using XPH
rather than HPX in the polar regions so it's better to be safe). */
if( fabs( astDistance( circle, centre, p0 ) <= 2*radius ) ) {

/* If the circle centre is above the northern transitional latitude, we
need to use the XPH projection centred on the north pole. */
if( centre[ 1 ] > AST__DD2R*SMF__HPX_TRANS ) {
msgOutif( MSG__VERB, "", "The Map will be created "
"using an XPH projection centred on the "
"north pole.", status );
usexph = 1;
break;

/* If the circle centre is below the southern transitional latitude, we
need to use the XPH projection centred on the south pole. */
} else if( centre[ 1 ] < -AST__DD2R*SMF__HPX_TRANS ) {
msgOutif( MSG__VERB, "", "The Map will be created "
"using an XPH projection centred on the "
"south pole.", status );
usexph = -1;
break;
}
}
}

/* If a projection was specified, use it. */
} else if( astChrMatch( text, "JSA-HPX" ) ) {
usexph = 0;
/* Get a list of the tiles that touch this circle. */
tiles = smf_jsatiles_region( circle, &skytiling,
&ntile, status );

} else if( astChrMatch( text, "JSA-XPHN" ) ) {
usexph = 1;
/* Choose the best projection (i.e. the projection that puts the circle
furthest away from any singularities). */
proj = smf_jsaproj( ntile, tiles, &skytiling, status);

} else if( astChrMatch( text, "JSA-XPHS" ) ) {
usexph = -1;
/* Free resources. */
tiles = astFree( tiles );
circle = astAnnul( circle );

} else if( *status == SAI__OK ) {
/* If a good projection was specified, use it. Otherwise report an error. */
} else if( proj == SMF__JSA_NULL && *status == SAI__OK ) {
*status = SAI__ERROR;
errRepf( "", "Bad value '%s' supplied for parameter %s.",
status, text, param );
}

/* Get the WCS FrameSet for what was the central tile in the old
raster-based indexing scheme - since we do not use the "local-origin"
option, the WCS for every tile is identical. The base Frame will be
/* Report the projection type. */
msgOutf( " ", "The %s will be created on the JSA %s "
"pixel grid.", status,
(data->hdr->instrument==INST__ACSIS)?"cube":"map",
smf_jsaproj_tostr( proj, status ) );

/* All tiles within the same JSA projection use the same WCS, so we get
the WCS FrameSet for an arbitrary central tile, and use it for the
full map. The exception is that tiles within the HPX facet that is
split between bottom-left and top-right, use a different WCS (they
have different reference points). But our choice of projection should
mean that the map never falls in that facet. The base Frame will be
GRID coords within the tile, and the current Frame will be ICRS
(RA,Dec). This used to be the central tile in numerical order, not
in spatial position. In other words, it was the lower of the two middle
tiles in numerical order. This tile is in the north west corner of
HPX facet 0, so in the new nested scheme, it has number 0 * ntpf^2 + all
the odd bits of the index within the tile. */
(RA,Dec). */
smf_jsatile( ((skytiling.ntpf * skytiling.ntpf - 1) * 2) / 3,
&skytiling, 0, usexph,
NULL, spacewcs, NULL, lbnd, ubnd, status );
&skytiling, 0, proj, NULL, spacewcs, NULL, lbnd,
ubnd, status );

/* Change the base Frame to be PIXEL. */
for( i = 1; i <= astGetI( *spacewcs, "NFrame" ); i++ ) {
Expand Down

0 comments on commit 7c1fbc0

Please sign in to comment.