Index: trunk/ppTranslate/src/ppMopsRead.c =================================================================== --- trunk/ppTranslate/src/ppMopsRead.c (revision 28243) +++ trunk/ppTranslate/src/ppMopsRead.c (revision 28623) @@ -4,5 +4,4 @@ #include -#include #include @@ -17,6 +16,5 @@ psArray *detections = psArrayAlloc(num); // Array of detections, to return for (int i = 0; i < num; i++) { - const char *name = inNames->data[i]; // File name - psFits *fits = psFitsOpen(name, "r"); // FITS file + psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file if (!fits) { psError(PS_ERR_IO, false, "Unable to open input %d", i); @@ -26,4 +24,10 @@ if (!header) { psError(PS_ERR_IO, false, "Unable to read header %d", i); + return false; + } + + psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image + if (diffSkyfileId == 0) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i); return false; } @@ -43,8 +47,7 @@ continue; } - ppMopsDetections *det = detections->data[i] = ppMopsDetectionsAlloc(); - det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension - det->header = header; - det->num = size; + ppMopsDetections *det = ppMopsDetectionsAlloc(size); + + psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); @@ -61,55 +64,132 @@ psMetadataLookupF32(NULL, header, "FWHM_MIN")); - det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); - det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); + int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows - det->psfHeader = psFitsReadHeader(NULL, fits); - if (!det->psfHeader) { - psError(psErrorCodeLast(), false, "Unable to read header %d", i); + psFree(header); + + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); return false; } - psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest - if (!table) { - psError(psErrorCodeLast(), false, "Unable to read table %d", i); - return false; + psFitsClose(fits); + + double plateScale = 0.0; // Plate scale + long numGood = 0; // Number of good rows + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + + psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); + if (flags & SOURCE_MASK) { + psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags); + continue; + } + + det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF"); + det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF")); + det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF")); + det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); + det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); + det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ"); + det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF"); + det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA"); + det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); + det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR"); + det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR"); + det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA"); + det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF"); + det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX"); + det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX"); + det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY"); + det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY"); + det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS"); + det->diffSkyfileId->data.S64[numGood] = diffSkyfileId; + det->naxis1->data.S32[numGood] = naxis1; + det->naxis2->data.S32[numGood] = naxis2; + + det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS"); + det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO"); + det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD"); + det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK"); + det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL"); + + // Calculate error in RA, Dec + double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); + double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); + double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); + double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); + + if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || + !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || + !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || + !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { + psTrace("ppMops.read", 10, + "Discarding row %ld from input %d because of non-finite values: " + "%f %f %lf %lf %f %f %f %f %f %f", + j, i, + det->x->data.F32[numGood], det->y->data.F32[numGood], + det->ra->data.F64[numGood], det->dec->data.F64[numGood], + det->mag->data.F32[numGood], det->magErr->data.F32[numGood], + xErr, yErr, scale, angle); + continue; + } + + // XXX Not at all sure I've got the angles around the right way here... + double cosAngle = cos(angle), sinAngle = sin(angle); + double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); + double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); + double errScale = scale / 3600.0; + det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[numGood] = 0; + plateScale += scale; + numGood++; + } + det->seeing *= plateScale / numGood; + + det->x->n = numGood; + det->y->n = numGood; + det->ra->n = numGood; + det->dec->n = numGood; + det->raErr->n = numGood; + det->decErr->n = numGood; + det->mag->n = numGood; + det->magErr->n = numGood; + det->chi2->n = numGood; + det->dof->n = numGood; + det->cr->n = numGood; + det->extended->n = numGood; + det->psfMajor->n = numGood; + det->psfMinor->n = numGood; + det->psfTheta->n = numGood; + det->quality->n = numGood; + det->numPix->n = numGood; + det->xxMoment->n = numGood; + det->xyMoment->n = numGood; + det->yyMoment->n = numGood; + det->flags->n = numGood; + det->diffSkyfileId->n = numGood; + det->naxis1->n = numGood; + det->naxis2->n = numGood; + det->mask->n = numGood; + det->nPos->n = numGood; + det->fPos->n = numGood; + det->nRatioBad->n = numGood; + det->nRatioMask->n = numGood; + det->nRatioAll->n = numGood; + + det->num = numGood; + + if (isfinite(args->zp) && numGood > 0) { + psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); } - psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); - psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); - det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); - det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); - if (!det->ra || !det->dec || !det->x || !det->y) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); - return false; - } - - psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); - - if (!psFitsMoveExtName(fits, "SkyChip.deteff")) { - psWarning("No detection efficiencies included in %s", det->component); - psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component); - psErrorClear(); - continue; - } - - det->deteffHeader = psFitsReadHeader(NULL, fits); - if (!det->deteffHeader) { - psWarning("Unable to read detection efficiency header in %s", det->component); - psErrorClear(); - continue; - } - det->deteffTable = psFitsReadTableAllColumns(fits); - if (!det->deteffTable) { - psFree(det->deteffHeader); - det->deteffHeader = NULL; - psWarning("Unable to read detection efficiency table in %s", det->component); - psErrorClear(); - continue; - } - psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component); - psFitsClose(fits); + psFree(table); + detections->data[i] = det; }