# -*-perl-*-
#
# The contents depend on whether we have bad-value support in PDL.
# (there's only a small amount of code that uses bad values)
#
# - this is called by Makefile.PL to ensure the file exists
# when the makefile is written
#
use strict;
use Config;
use File::Basename qw(&basename &dirname);
# check for bad value support
use vars qw( $bvalflag $usenan );
use File::Spec;
require File::Spec->catfile( File::Spec->updir, File::Spec->updir, "Basic", "Core", "badsupport.p" );
# This forces PL files to create target in same directory as PL file.
# This is so that make depend always knows where to find PL derivatives.
chdir(dirname($0));
my $file;
($file = basename($0)) =~ s/\.PL$//;
$file =~ s/\.pl$//
if ($Config{'osname'} eq 'VMS' or
$Config{'osname'} eq 'OS2'); # "case-forgiving"
if ( $bvalflag ) {
print "Extracting $file (WITH bad value support)\n";
} else {
print "Extracting $file (NO bad value support)\n";
}
open OUT,">$file" or die "Can't create $file: $!";
chmod 0644, $file;
#########################################################
print OUT <<"!WITH!SUBS!";
# This file is automatically created by NDF.pm.PL
# - bad value support = $bvalflag
!WITH!SUBS!
print OUT <<'!NO!SUBS!';
package PDL::IO::NDF;
=head1 NAME
PDL::IO::NDF - PDL Module for reading and writing Starlink
N-dimensional data structures as PDLs.
=head1 SYNOPSIS
use PDL::IO::NDF;
$a = PDL->rndf($file);
$a = rndf('test_image');
$a = rndf('test_image', 1);
$a->wndf($file);
wndf($a, 'out_image');
propndfx($a, 'template', 'out_image');
=head1 DESCRIPTION
This module adds the ability to read and write Starlink N-dimensional data
files as N-dimensional PDLs.
You must have the Starlink NDF library installed to use it. The library
is distributed under the GPL and is available from "http://www.starlink.ac.uk".
=head1 FUNCTIONS
=cut
@EXPORT_OK = qw/rndf wndf/;
%EXPORT_TAGS = (Func=>[@EXPORT_OK]);
@ISA = qw( PDL::Exporter );
use PDL::Core;
use PDL::Types;
use Carp;
use strict;
# Starlink data type conversion
use vars qw/%pdltypes %startypes $ndf_loaded $VERSION $EXTNAME/;
$VERSION = '1.02';
# Set PDL -> Starlink data types
%pdltypes = ("$PDL_B" => "_UBYTE", # was "_BYTE"
"$PDL_S" => "_WORD",
"$PDL_US" => "_UWORD",
"$PDL_L" => "_INTEGER",
"$PDL_F" => "_REAL",
"$PDL_D" => "_DOUBLE"
);
# Set Starlink -> PDL data types
%startypes = ("_BYTE" =>$PDL_B,
"_UBYTE" =>$PDL_B,
"_WORD" =>$PDL_S,
"_UWORD" =>$PDL_US,
"_INTEGER" =>$PDL_L,
"_REAL" =>$PDL_F,
"_DOUBLE" =>$PDL_D
);
# This is the name I use in the PDL header hash that contains the
# NDF extensions so that they can be recreated
$EXTNAME = 'NDF_EXT';
# load NDF library
sub _init_NDF {
return if $ndf_loaded++; # a bit pointless increasing each time?
eval 'use NDF';
croak 'Cannot use NDF library' if $@ ne "";
}
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
use PDL::Bad;
# internal functions for checking bad values
my $starbad = 0; # flag
use vars qw( @starbad );
sub _init_starbad {
return if $starbad;
_init_NDF(); # just to be safe
@starbad = eval
"( &NDF::VAL__BADUB,
&NDF::VAL__BADW,
&NDF::VAL__BADUW,
&NDF::VAL__BADI,
&NDF::VAL__BADF,
&NDF::VAL__BADD );
";
croak 'Unable to initialise bad values from NDF library' if $@ ne "";
$starbad = 1;
}
sub starlink_bad_value ($) {
_init_starbad();
return $starbad[$_[0]];
}
# given a piddle, check if the Starlink bad value matches
# the bad value used for this type
sub compare_bad_values ($) {
_init_starbad();
return ($starbad[$_[0]] == badvalue($_[0]));
}
!NO!SUBS!
} # if: $bvalflag
print OUT <<'!NO!SUBS!';
=head2 rndf()
=for ref
Reads a piddle from a NDF format data file.
=for example
$pdl = rndf('file.sdf');
$pdl = rndf('file.sdf',1);
The '.sdf' suffix is optional. The optional second argument turns off
automatic quality masking and returns a quality array as well.
BUG? C<rndf('bob.sdf',1)> calls ndf_sqmf(1,...), which means that
the quality array is turned into bad pixels - ie the I<opposite> of
above. Or am I confused?
Header information and NDF Extensions are stored in the piddle as a hash
which can be retreived with the C<$pdl-E<gt>gethdr> command.
Array extensions are stored in the header as follows:
$a - the base DATA_ARRAY
If C<$hdr = $a-E<gt>gethdr>;
then:
%{$hdr} contains all the FITS headers plus:
$$hdr{Error} contains the Error/Variance PDL
$$hdr{Quality} The quality byte array (if reqeusted)
@{$$hdr{Axis}} Is an array of piddles containing the information
for axis 0, 1, etc.
$$hdr{NDF_EXT} Contains all the NDF extensions
$$hdr{Hist} Contains the history information
$$hdr{NDF_EXT}{_TYPES}
Data types for non-PDL NDF extensions so that
wndf can reconstruct a NDF.
All extension information is stored in the header hash array.
Extension structures are preserved in hashes, so that the PROJ_PARS
component of the IRAS.ASTROMETRY extension is stored in
C<$$hdr{NDF_EXT}{IRAS}{ASTROMETRY}{'PROJ_PARS'}>. All array structures are
stored as arrays in the Hdr: numeric arrays are stored as PDLs,
logical and character arrays are stored as plain Perl arrays. FITS
arrays are a special case and are expanded as scalars into the header.
PDL does not have a signed byte datatype, so any '_BYTE' data
is read into a C<byte> (unsigned) piddle and a warning is printed
to C<STDOUT>.
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
=for bad
If the starlink bad flag is set, then the bad flag on the output
piddle is also set. Starlink bad values are converted to the
current bad value used by the piddle type (if they are different).
!NO!SUBS!
} # if: $bvalflag
print OUT <<'!NO!SUBS!';
=cut
# This is one form of the new command
sub rndf {PDL->rndf(@_)}
# And this is the real form
# Allows the command to be called in OO form or as a function
sub PDL::rndf { # Read a piddle from a NDF file
my $class = shift;
barf 'Usage: $a = rndf($file,[$nomask]); $a = PDL->rndf(...) '
if ($#_!=0 && $#_!=1);
# ensure we have the NDF module
_init_NDF();
# Setup the Header Hash
my $header = {};
# Read in the filename
# And setup the new PDL
my $file = shift;
my $pdl = $class->new;
my $nomask = shift if $#_ > -1;
my ($infile, $status, $indf, $entry, @info, $value);
# Strip trailing .sdf if one is present
# remove the trailing extension names
# (anything after a `.' UNLESS it's part of a directory
# name, hence the check for ^/) and any trailing (...)
# specifier (NDF sections, a la PDL slices)
#
$file =~ s/\.sdf$//;
$infile = $file;
$infile =~ s!\.[^/]*$!!;
$infile =~ s/\(.*\)$//;
# If file is not there
croak "Cannot find $infile.sdf" unless -e "$infile.sdf";
# Set status
$status = &NDF::SAI__OK;
# Begin NDF and ERR session
err_begin($status);
ndf_begin();
# Open the file
ndf_find(&NDF::DAT__ROOT, $file, $indf, $status);
# unset automatic quality masking if $nomask is defined
# I think the logic here is wrong (ie 0 and 1 should be swapped)
$nomask = 0 unless (defined $nomask);
$nomask = 1 if $nomask != 0;
ndf_sqmf($nomask, $indf, $status);
# Read the data array (need to pass the header as well)
my $badflag;
( $status, $badflag ) = rdata($indf, $pdl, $nomask, $header, $class, $status);
# Read the axes
raxes($indf, $pdl, $header, $class, $status);
# Read the header
rhdr($indf, $pdl, $header, $class, $status);
# Read history information
rhist($indf, $pdl, $header, $status);
# Read labels
@info = ('Label', 'Title', 'Units');
for $entry (@info) {
$value = 'NULL';
ndf_cget($indf, $entry, $value, $status);
$$header{"$entry"} = $value if $value ne 'NULL';
print "$entry is $value\n" if ($PDL::verbose && $value ne 'NULL');
undef $value;
}
# Tidy up
ndf_annul($indf, $status);
ndf_end($status);
if ($status != &NDF::SAI__OK) {
err_flush($status); # Flush the error system
err_end($status);
croak("rndf: obtained NDF error");
}
err_end($status);
# Put the header into the main pdl
$pdl->sethdr($header);
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
# set the bad flag if it was set in the input file
print "Bad flag status is: $badflag\n" if $PDL::verbose;
$pdl->badflag(1) if $badflag;
!NO!SUBS!
} # if; $bvalflag
print OUT <<'!NO!SUBS!';
return $pdl;
}
=head2 wndf()
=for ref
Writes a piddle to a NDF format file:
=for example
$pdl->wndf($file);
wndf($pdl,$file);
wndf can be used for writing PDLs to NDF files.
The '.sdf' suffix is optional. All the extensions
created by rndf are supported by wndf. This means that error, axis
and quality arrays will be written if they exist. Extensions are also
reconstructed by using their name (ie FIGARO.TEST would be expanded as
a FIGARO extension and a TEST component). Hdr keywords Label, Title
and Units are treated as special cases and are written to the label,
title and units fields of the NDF.
Header information is written to corresponding NDF extensions.
NDF extensions can also be created in the {NDF} hash by using a key
containing '.', ie {NDF}{'IRAS.DATA'} would write the information to
an IRAS.DATA extension in the NDF. rndf stores this as
C<$$hdr{NDF}{IRAS}{DATA}> and the two systems are interchangeable.
rndf stores type information in {NDF}{'_TYPES'} and below so that
wndf can reconstruct the data type of non-PDL extensions. If no
entry exists in _TYPES, wndf chooses between characters, integer and
double on a best guess basis. Any perl arrays are written as CHAR
array extensions (on the assumption that numeric arrays will exist as
PDLs).
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
=for bad
The bad flag is written to the file, if set.
!NO!SUBS!
} # if: $bvalflag
print OUT<<'!NO!SUBS!';
=cut
# This is one form of the new command
# OO version.
*wndf = \&PDL::wndf;
# And this is the real form
# Allows the command to be called in OO form or as a function
sub PDL::wndf { # Write a PDL to an NDF format file
barf 'Usage: wndf($pdl,$file)' if $#_!=1;
# ensure we have the NDF module
_init_NDF();
my ($indf, $place, $status, $outndf);
my (@lbnd, @ubnd);
my ($pdl, $outfile) = @_;
# Check that we are dealing with a PDL
croak 'Argument is not a PDL variable' unless $pdl->isa('PDL');
# Set status
$status = &NDF::SAI__OK;
# Strip trailing .sdf if one is present
$outfile =~ s/\.sdf$//;
# Begin context
ndf_begin();
err_begin($status);
# Open the new file
ndf_place(&NDF::DAT__ROOT(), $outfile, $place, $status);
# Need to create data_array component
# Change the bounds to match the PDL
@ubnd = $pdl->dims;
@lbnd = ();
@lbnd = map { 1 } (0..$#ubnd);
ndf_new($pdltypes{$pdl->get_datatype}, $#ubnd+1, \@lbnd, \@ubnd, $place,
$outndf, $status);
# Set the application name
ndf_happn('PDL::wndf', $status);
# Write the data to this file
wdata($outndf, $pdl, $status);
# Write the header and title
whdr($outndf, $pdl, $status);
# Create a history extension
ndf_hcre($outndf, $status);
# Annul the identifier
ndf_annul($outndf, $status);
# End context
ndf_end($status);
if ($status != &NDF::SAI__OK) {
err_flush($status);
croak "Bad STATUS in wndf\n";
}
err_end($status);
return 1;
}
=head2 propndfx()
Routine to write a PDL to an NDF by copying the extension information
from an existing NDF and writing DATA,VARIANCE, QUALITY and AXIS info
from a PDL (if they exist).
Extensions, labels and history are propogated from the old NDF.
No new extension information is written.
This command has been superseded by L<wndf()|/wndf()>.
=cut
*propndfx = \&PDL::propndfx;
sub PDL::propndfx { # Write a PDL to a NDF format file
barf 'Usage: propndfx($pdl, $infile, $outfile)' if $#_!=2;
# ensure we have the NDF module
_init_NDF();
my ($indf, $in_place, $status, $outplace, $outndf);
my ($pdl, $infile, $outfile) = @_;
# Set status
$status = &NDF::SAI__OK;
# Strip trailing .sdf if one is present
$outfile =~ s/\.sdf$//;
$infile =~ s/\.sdf$//;
# File is the first thing before a .
my $file = (split(/\./,$infile))[0];
croak "$file does not exist\n" unless -e "$file.sdf";
# Begin context
ndf_begin();
err_begin($status);
# Open the old file
ndf_find(&NDF::DAT__ROOT(), $infile, $indf, $status);
# Set the application name
ndf_happn('PDL::propndfx', $status);
# Open the new file
ndf_place(&NDF::DAT__ROOT(), $outfile, $outplace, $status);
# Copy the extension information to outfile
ndf_scopy($indf, ' ', $outplace, $outndf, $status);
# Annul the input file identifier
ndf_annul($indf, $status);
# Write the data to this file
wdata($outndf, $pdl, $status);
# Annul the identifier
ndf_annul($outndf, $status);
# End context
ndf_end($status);
if ($status != &NDF::SAI__OK) {
err_flush($status);
croak "Bad STATUS in propndfx\n";
}
err_end($status);
return 1;
}
=head1 NOTES
The perl NDF module must be available. This is available from the
author or from Starlink (http://www.starlink.rl.ac.uk).
If an NDF is read which contains AST World Coordinate information
(a .WCS component) this information is currently ignored. Currently
WCS information can only be written and stored using standard FITS headers.
See http://rlspc5.bnsc.rl.ac.uk/star/docs/sun211.htx/sun211.html#xref_
for more information on AST.
=head1 AUTHOR
This module was written by Tim Jenness E<lt>t.jenness@jach.hawaii.eduE<gt>.
Copyright (C) Tim Jenness 1997-2000. All Rights Reserved.
=head1 SEE ALSO
L<PDL::FAQ> for general information on the Perl Data language,
L<NDF> for information on the NDF module.
=cut
#######################################################################
# These are the generic routines used by the rndf command
# RDATA
# Read Data, Quality and Variance
# returns ( $status, $badflag ), where $badflag is:
#
# 1 - the piddle may contain bad data
# 0 - the piddle contains no bad data
#
# This reports starlink's bad-data flag (ie it's valid
# even in bad value support is not available in PDL)
#
sub rdata {
my ($indf, $pdl, $nomask, $header, $class, $status) = @_;
return ( $status, undef ) if $status != &NDF::SAI__OK;
my ($maxdims, $ndim, @dim, @comps, $dcomp, $tcomp, $exist);
my ($type, $data_pntr, $el, $temppdl, $nbytes, $badbit, $dref);
####################################################################
# D I M S #
####################################################################
# Find the dimensions of the data array
$maxdims = 100; # Maximum number of dimensions allowed
&ndf_dim($indf, $maxdims, \@dim, $ndim, $status);
@dim = @dim[0..$ndim-1];
print "Dims = @dim\n" if $PDL::verbose;
$$header{Extensions} = [];
####################################################################
# D A T A + V A R I A N C E + Q U A L I T Y #
####################################################################
my $badflag = 0;
@comps = ('Data','Error');
push(@comps, 'Quality') if $nomask;
for $dcomp (@comps) {
$tcomp = $dcomp;
$tcomp = 'Variance' if $tcomp eq 'Error';
ndf_state($indf, $tcomp, $exist, $status);
if ($exist && ($status == &NDF::SAI__OK)) {
# Find the type of the data array
ndf_type($indf, $tcomp, $type, $status);
# Map the data array
ndf_map($indf, $dcomp, $type, 'READ', $data_pntr, $el, $status);
# Check the bad-data status of this component
my $this_badflag = 0;
my $ecomp = ( $dcomp =~ /Err/ ? "Variance" : $dcomp);
ndf_bad($indf, $ecomp, 0, $this_badflag, $status );
if ($status == &NDF::SAI__OK) {
print "Reading $dcomp array...\n" if $PDL::verbose;
# combine the bad data flag with the "global" one for
# the file
$badflag |= ($this_badflag != 0);
# Create a temporary PDL to deal with separate components
if ($dcomp eq 'Data') {
$temppdl = $pdl;
} else {
$dcomp = 'Errors' if $dcomp eq 'Error';
$$header{"$dcomp"} = $class->new;
$temppdl = $$header{"$dcomp"};
}
# Set up the PDL
$temppdl->set_datatype($startypes{$type});
$temppdl->setdims([@dim]);
$dref = $temppdl->get_dataref;
# warning about data loss
if ( $type eq "_BYTE" ) {
print "Warning: NDF type (signed byte) converted to PDL's unsigned byte.\n";
}
# How many bytes in this data type?
$type =~ s/^_|^_U//;
$nbytes = &NDF::byte_size($type);
mem2string($data_pntr, $nbytes*$el, $$dref);
push(@{$$header{Extensions}}, $dcomp) if $dcomp ne 'Data';
if ($dcomp eq 'Quality') {
# Quality bad bit mask
ndf_bb($indf, $badbit, $status);
my $qhdr = {};
$$qhdr{BADBIT} = $badbit;
$temppdl->sethdr($qhdr);
}
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
# if badflag is true, convert any STARLINK bad values
# to PDL values
if ( $badflag and
!compare_bad_values($temppdl->get_datatype) ) {
my $star_bad = starlink_bad_value($temppdl->get_datatype);
print "Converting from STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
$temppdl->inplace->setvaltobad( $star_bad );
}
!NO!SUBS!
} # if: $bvalflag
print OUT <<'!NO!SUBS!';
# Flush and Clear temporary PDL
$temppdl->upd_data();
undef $temppdl;
}
ndf_unmap($indf, $tcomp, $status);
}
}
err_flush($status) if $status != &NDF::SAI__OK;
return ( $status, $badflag );
}
####################################################################
# A X I S #
####################################################################
# Axes are stored as follows:
# Array of PDL axes as an array in header of main pdl (called 'Axis')
# Header of this PDL contains label and errors and units etc.
sub raxes {
my ($indf, $pdl, $header, $class, $status) = @_;
return $status if $status != &NDF::SAI__OK;
my ($there, $naxis, @dims, $axcomp, $exist, $axtype, $axpntr, $el);
my ($nbytes, $temppdl, $tcomp, $daxref, $dref);
my ($axhdr, $ndims);
# Read in axis information
ndf_state($indf, 'Axis', $there, $status);
$ndims = $pdl->getndims; # Find number of dimensions
if ($there) {
# Label from 0..$#dims to match array index handling
# Want to put the axes into an array of axes
# And update the extensions field
$$header{Axis} = [];
push(@{$$header{Extensions}}, "Axis");
for $naxis (0..$ndims-1) {
# Does a CENTRE component exist
$axcomp = 'CENTRE';
ndf_astat($indf, $axcomp, $naxis+1, $exist, $status);
if ($exist && ($status == &NDF::SAI__OK)) {
print "Reading Axis $naxis...\n" if $PDL::verbose;
ndf_atype($indf, $axcomp, $naxis+1, $axtype, $status);
ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
$el, $status);
# Set up new PDL for axis info if map was okay
if ($status == &NDF::SAI__OK) {
print "Number of elements: $el\n" if $PDL::verbose;
$$header{Axis}[$naxis] = $class->new;
$temppdl = $$header{Axis}[$naxis];
$temppdl->set_datatype($startypes{$axtype});
$temppdl->setdims([$el]);
$dref = $temppdl->get_dataref;
# How many bytes?
$axtype =~ s/^_|^_U//;
$nbytes = &NDF::byte_size($axtype);
mem2string($axpntr, $nbytes*$el, $$dref);
$temppdl->upd_data();
# Header info for axis
$axhdr = {}; # somewhere to put the pdl
# Read VARIANCE array if there
$axcomp = 'Error';
$tcomp = 'Variance';
ndf_astat($indf, $tcomp, $naxis+1, $exist, $status);
if ($exist && ($status == &NDF::SAI__OK)) {
print "Reading Axis $naxis errors...\n" if $PDL::verbose;
ndf_atype($indf, $tcomp, $naxis+1, $axtype, $status);
ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
$el, $status);
# Set up new PDL for axis info if map was okay
if ($status == &NDF::SAI__OK) {
$$axhdr{Errors} = $class->new;
$$axhdr{Extensions} = ["Errors"];
$$axhdr{Errors}->set_datatype($startypes{$axtype});
$$axhdr{Errors}->setdims([$el]);
$daxref = $$axhdr{Errors}->get_dataref;
# How many bytes?
$axtype =~ s/^_|^_U//;
$nbytes = &NDF::byte_size($axtype);
mem2string($axpntr, $nbytes*$el, $$daxref);
$$axhdr{Errors}->upd_data();
}
}
# Get label and units
for my $entry ('Units', 'Label') {
ndf_astat($indf, $entry, $naxis+1, $exist, $status);
if ($exist) {
my $value = '';
ndf_acget($indf,$entry, $naxis+1, $value, $status);
$axhdr->{"$entry"} = $value;
}
}
}
ndf_aunmp($indf,'*',$naxis+1,$status);
# Now store this header into temppdl
$temppdl->sethdr($axhdr);
} else {
err_annul($status) unless $status == &NDF::SAI__OK;
}
}
}
return $status;
}
####################################################################
# E X T E N S I O N S #
####################################################################
sub rhdr {
my ($indf, $pdl, $header, $class, $status) = @_;
return $status if $status != &NDF::SAI__OK;
my ($nextn, $xname, $xloc, $extn, $size, @fits, $entry, $value);
my ($comment, $item, $el);
# Read in any header information from the extensions
ndf_xnumb($indf, $nextn, $status);
print "Reading in header information...\n"
if (($nextn > 0) && $PDL::verbose);
# Read in extensions one at a time
for $extn (1..$nextn) {
$xname = 'NULL';
ndf_xname($indf, $extn, $xname, $status);
undef $xloc;
ndf_xloc($indf, $xname, 'READ', $xloc, $status);
# FITS is a special case and must be expanded
if ($xname eq 'FITS') {
# Read in the data
dat_size($xloc, $size, $status);
&dat_getvc($xloc, $size, \@fits, $el, $status);
if ($status == &NDF::SAI__OK) {
print "Reading FITS header...\n" if $PDL::verbose;
for $entry (0..$el-1) {
($item, $value, $comment) = &fits_get_nth_item(\@fits,$entry);
$$header{$item} = $value;
$$header{'_COMMENTS'}{$item} = $comment;
}
} else {
print "Error mapping $xname\n";
err_annul($status);
}
undef @fits;
} else {
# Read in extensions to $EXTNAME
$status = &find_prim($xloc, \%{$$header{$EXTNAME}}, $class, $status);
}
dat_annul($xloc, $status);
}
return $status;
}
####################################################################
# H I S T O R Y #
####################################################################
sub rhist {
my ($indf, $pdl, $header, $status) = @_;
return $status if $status != &NDF::SAI__OK;
my ($nrec, $app, $date, $i);
# Read HISTORY information into a Hist header
ndf_hnrec($indf,$nrec,$status);
if ($status == &NDF::SAI__OK && ($nrec > 0)) {
print "Reading history information into 'Hist'...\n" if $PDL::verbose;
$$header{Hist}{'Nrecords'} = $nrec;
# Loop through the history components and find last "APPLICATION"
for $i (1..$nrec) {
ndf_hinfo($indf, 'APPLICATION', $i, $app, $status);
ndf_hinfo($indf, 'DATE', $i, $date, $status);
$$header{Hist}{"Application$i"} = $app;
$$header{Hist}{"Date$i"} = $date;
print " $app on $date\n" if $PDL::verbose;
}
} else {
err_annul($status);
}
return $status;
}
################ Generic routines ###########################
#################### INTERNAL ###############################
# Find primitive components below a given HDS locator
# FITS information is stored in Hdr
# NDF extensions stored in NDF
sub find_prim {
my ($loc, $href, $class, $status) = @_;
my ($prim, $type, $size, @dim, $ndims, $struct, $name, $nloc, $el);
my ($value, @values, $item, $entry, $maxdims);
my ($packtype, $ncomp, $packed, $comp);
my ($pntr, $nbytes, $comment, $dref);
# Return if bad status
return 0 if $status != &NDF::SAI__OK;
# Most dimensions
$maxdims = 100;
# Is it a primitive component
dat_prim($loc, $prim, $status);
# Type, size and shape
dat_type($loc, $type, $status);
dat_size($loc, $size, $status);
dat_shape($loc, $maxdims, \@dim, $ndims, $status);
dat_name($loc, $name, $status);
if ($prim) {
print "\tReading component: $name ($type) " if $PDL::verbose;
# Store type of information
$$href{"_TYPES"}{"$name"} = $type;
# Read it into the header
if (($size == 1) && ($ndims == 0)) {
if ($type eq '_DOUBLE') {
# dat_get0c doesn't return full precision for doubles
dat_get0d($loc, $value, $status);
if ($status != &NDF::SAI__OK) {
print "Error mapping scalar 0d $type $name\n";
err_flush($status);
}
} else {
dat_get0c($loc, $value, $status);
if ($status != &NDF::SAI__OK) {
print "Error mapping scalar 0c $name\n";
err_flush($status);
}
}
$$href{"$name"} = $value if ($status == &NDF::SAI__OK);
undef $value;
} else {
@dim = @dim[0..$ndims-1];
# Inform the user of dimensions
print "[@dim]" if $PDL::verbose;
# Map arrays (don't need to unpack them)
if ($type =~ /CHAR|LOG/) {
# Read in the data
dat_getvc($loc, $size, \@values, $el, $status);
if ($status != &NDF::SAI__OK) {
print "Error mapping $name\n";
} else {
$$href{"$name"} = [@values];
undef @values;
}
} else {
# map the data
dat_mapv($loc, $type, 'READ', $pntr, $el, $status);
if ($status != &NDF::SAI__OK) {
print "Error mapping $name\n";
}
if ($status == &NDF::SAI__OK) {
# Create the pdl
$$href{"$name"} = $class->new;
$$href{"$name"}->set_datatype($startypes{$type});
$$href{"$name"}->setdims([@dim]);
$dref = $$href{"$name"}->get_dataref();
# How many bytes in this data type?
$type =~ s/^_|^_U//;
$nbytes = &NDF::byte_size($type);
mem2string($pntr, $nbytes*$el, $$dref);
$$href{"$name"}->upd_data();
}
dat_unmap($loc,$status);
}
}
print "\n" if $PDL::verbose;
} else {
my ($ndimx, @dim, $ndim);
dat_shape($loc, 20, \@dim, $ndim, $status);
# If we have a single dimension ($ndim=0) then
# proceed. Cant yet cope with multi-dimensional
# extensions
if ($ndim == 0) {
dat_ncomp($loc, $ncomp, $status);
print "$name ($type) has $ncomp components\n" if $PDL::verbose;
# Store type of extension
$$href{$name}{"_TYPES"}{STRUCT}{"$name"} = $type;
for $comp (1..$ncomp) {
dat_index($loc, $comp, $nloc, $status);
$status = &find_prim($nloc, \%{$$href{$name}}, $class, $status);
dat_annul($nloc, $status);
}
} else {
print "Extension $name is an array structure -- skipping\n"
if $PDL::verbose;
}
}
return $status;
}
###### Routines for WRITING data ######################################
sub wdata {
my ($outndf, $pdl, $status) = @_;
my (@bounds, $ndims, @lower, @comps, $dcomp, $temppdl, $type);
my ($pntr, $el, $nbytes, $axis, $axcomp, $axpntr, %hdr, %axhdr);
my ($entry, $value, $i, $axndims, $axtype, $tcomp, @acomps);
# Return if bad status
return 0 if $status != &NDF::SAI__OK;
####################################################################
# D I M S #
####################################################################
# Change the bounds to match the PDL
@bounds = $pdl->dims;
$ndims = $#bounds;
@lower = map { 1 } (0..$ndims);
&ndf_sbnd($ndims+1, \@lower, \@bounds, $outndf, $status);
####################################################################
# D A T A + V A R I A N C E + Q U A L I T Y #
####################################################################
# Map data, variance, quality...
@comps = ('Data','Errors','Quality');
# Retrieve header and check that the header is a hash ref
my $hdr = $pdl->gethdr;
if (ref($hdr) eq 'HASH') {
%hdr = %$hdr;
} else {
%hdr = ();
}
for $dcomp (@comps) {
if ($dcomp eq 'Data') {
$temppdl = $pdl;
} else {
if (exists $hdr{"$dcomp"}) {
$temppdl = $hdr{"$dcomp"};
# Check that we have a PDL
next unless UNIVERSAL::isa($temppdl, 'PDL');
} else {
next; # Skip this component
}
}
# Check that we have some data
if ($temppdl->nelem > 0) {
if (($dcomp eq 'Quality') && ($temppdl->get_datatype != $PDL_B)) {
# Quality must be bytes so convert to BYTE if this is not true
$temppdl = byte ($$pdl{"$dcomp"});
}
$type = $pdltypes{$temppdl->get_datatype};
$type = '_UBYTE' if $dcomp eq 'Quality'; # No choice for QUALITY
# Set the output type of the NDF
$tcomp = $dcomp;
$tcomp = 'Variance' if $tcomp eq 'Errors';
ndf_stype($type, $outndf, $tcomp, $status) unless $dcomp eq 'Quality';
print "Mapping $dcomp , $type\n" if $PDL::verbose;
# Map NDF
$dcomp = 'Error' if $dcomp eq 'Errors';
ndf_map($outndf, $dcomp, $type, 'WRITE', $pntr, $el, $status);
if ($status == &NDF::SAI__OK) {
# Number of bytes per entry
$nbytes = PDL::Core::howbig($temppdl->get_datatype);
# Convert to 1D data stream
my $p1d = $temppdl->clump(-1);
# Total number of bytes
$nbytes *= $p1d->getdim(0);
!NO!SUBS!
if ( $bvalflag ) {
print OUT <<'!NO!SUBS!';
# if badflag is true, convert any PDL bad values to STARLINK ones
# note: we have to create a copy of the data to ensure that the
# bad data value of the original piddle does not get changed.
# per-piddle bad values would obviate the need for this.
if ( $temppdl->badflag() ) {
print "Setting bad flag for component $dcomp ...\n" if $PDL::verbose;
ndf_sbad( 1, $outndf, $dcomp, $status );
if ( !compare_bad_values($temppdl->get_datatype) ) {
my $star_bad = starlink_bad_value($temppdl->get_datatype);
print "Converting from PDL to STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
$temppdl = $temppdl->setbadtoval( $star_bad ); # do NOT do inplace
}
}
!NO!SUBS!
} # if: $bvalflag
print OUT <<'!NO!SUBS!';
# Copy to disk
string2mem(${$temppdl->get_dataref}, $nbytes, $pntr);
# Write badbit mask
if ($dcomp eq 'Quality') {
ndf_sbb($$temppdl{Hdr}{BADBIT}, $outndf, $status)
if defined $$temppdl{Hdr}{BADBIT};
}
}
# Unmap Data
ndf_unmap($outndf, $tcomp, $status);
}
# free the temporary pdl
undef $temppdl;
}
####################################################################
# A X I S #
####################################################################
# A X I S information is expected as an array in the header
# called 'Axis'. These PDLs contain the CENTRE data and any further
# info is stored in their header.
# @bounds is accessed to find expected shape
# Simply look in the header for a Axis
if (exists $hdr{Axis}) {
# Check that we have an array
if (ref($hdr{Axis}) eq 'ARRAY') {
# Now loop over axes
for my $i (0..$#{$hdr{Axis}} ) {
# Loop unless status is bad
last unless $status == &NDF::SAI__OK;
# Extract the ith axis PDL from the array
my $axis = $hdr{Axis}->[$i];
# If we have a PDL
if (UNIVERSAL::isa($axis, 'PDL')) {
# We now want to copy the data and if necessary the
# Error array. Since there are only two I will do it the
# long way by explcitly storing data and then error
# Set data type
$axtype = $pdltypes{$axis->get_datatype};
ndf_astyp($axtype, $outndf, 'Centre', $i+1, $status);
# Okay we can now map this pdl
ndf_amap($outndf, 'Centre', $i+1, $axtype, 'WRITE', $axpntr,
$el, $status);
# Check that we have the correct number of entries
my $nelem = $axis->nelem;
if ($el == $nelem) {
print "Mapping axis " , $i+1 , "\n" if $PDL::verbose;
# Number of bytes per entry
$nbytes = PDL::Core::howbig($axis->get_datatype) * $el;
# Copy to disk
string2mem( $ { $axis->get_dataref }, $nbytes, $axpntr)
if ($status == &NDF::SAI__OK);
} else {
carp "Axis ",$i+1 .
" is the wrong size ($el values required but got ".
$nelem . ")- ignoring";
}
# Unmap
ndf_aunmp($outndf, '*', $i+1, $status);
# Errors
my $axhdr = $axis->gethdr; # Retrieve and check header
if (ref($axhdr) eq 'HASH') {
%axhdr = %$axhdr;
} else {
%axhdr = ();
}
# Look for an Errors component in the header hash
if (exists $axhdr{Errors}) {
my $axerr = $axhdr{Errors};
if (UNIVERSAL::isa($axerr, 'PDL')) {
# Set data type
$axtype = $pdltypes{$axerr->get_datatype};
ndf_astyp($axtype, $outndf, 'Variance', $i+1, $status);
# Okay we can now map this pdl
ndf_amap($outndf, 'Errors', $i+1, $axtype, 'WRITE',
$axpntr, $el, $status);
# Check that we have the correct number of entries
my $nelem = $axerr->nelem;
print "Nelem: $nelem and $el\n";
if ($nelem == $el) {
print "Mapping errors for axis " . $i+1 . "\n"
if $PDL::verbose;
# Number of bytes per entry
$nbytes = PDL::Core::howbig($axerr->get_datatype) * $el;
# Copy to disk
string2mem($ {$axerr->get_dataref}, $nbytes, $axpntr)
if ($status == &NDF::SAI__OK);
} else {
carp "Error PDL for Axis ",$i+1,
" is the wrong size ($el values required but got ".
$axerr->nelem . ")- ignoring";
}
# Unmap
ndf_aunmp($outndf, '*', $i+1, $status);
}
}
# Now character components
foreach my $char (keys %axhdr) {
if ($char =~ /LABEL|UNITS/i) {
$value = $axhdr{"$char"};
ndf_acput($value, $outndf, $char, $i+1, $status)
if length($value) > 0;
}
}
}
}
}
}
}
# Write header information to NDF extension
# This routine does not reconstruct an NDF identically to one that
# was read in. FITS extensions are made. All non-integer numeric types
# are written as doubles (apart from PDLs which carry there own type
# information) unless _TYPES information exists in {NDF}.
sub whdr {
my ($outndf, $pdl, $status) = @_;
my (%header, @fitsdim, $fitsloc, $value);
my (%unused, @fits, $hashref, $hdr);
# Return if bad status
return 0 if $status != &NDF::SAI__OK;
print "Writing header information...\n" if $PDL::verbose;
# Write FITS header from {Hdr}
# Retrieve and check header from PDL
$hdr = $pdl->gethdr;
if (ref($hdr) eq 'HASH') {
%header = %$hdr;
} else {
%header = ();
}
%unused = ();
@fits = ();
foreach my $key (sort keys %header) {
next if $key eq '_COMMENTS';
if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
# This is not extension info
ndf_cput($header{$key}, $outndf, $key, $status)
if length($header{$key}) > 0;
}
# Only write scalars
if (not ref($header{"$key"})) {
push(@fits, fits_construct_string($key,$header{$key},
$header{'_COMMENTS'}{$key}) );
} else {
# Deal with later
$unused{$key} = $header{$key};
}
}
# Write the FITS extension
if ($#fits > -1) {
push(@fits, "END");
# Write it out
@fitsdim = ($#fits+1);
ndf_xnew($outndf, 'FITS', '_CHAR*80', 1, @fitsdim, $fitsloc, $status);
dat_put1c($fitsloc, $#fits+1, \@fits, $status);
dat_annul($fitsloc, $status);
}
# Loop through all NDF header info and array Hdr information
# Note that %unused still contains our NDF hash so
# we must remove that before proceeding (else we end up with
# an 'NDF_EXT' extension!
delete $unused{$EXTNAME};
foreach $hashref (\%unused, $header{$EXTNAME}) {
whash($hashref, $outndf, '', '', $status);
}
}
#-----------------------------------------------------
# This routine writes a hash array to a NDF extension
# Keys that have leading underscores are skipped
# $hashname contains the extension name to use by default.
# This is extended by any '.' in the key.
sub whash {
my ($hash, $outndf, $hashname, $stypes, $status) = @_;
my ($key, $comp, @structures, $extension, $xtype, @dim, $ndims, $loc);
my (@locs, $there, $type, %header, $packed, $length);
my (@bounds, $nbytes, $pntr, $maxsize, $which, $el, $path);
my ($oldhash, $oldtypes, $structs);
if (defined $hash) {
%header = %$hash;
} else {
%header = ();
}
my $root = 1 if length($hashname) == 0; # Mark a ROOT structure
# Return if bad status
return 0 if $status != &NDF::SAI__OK;
foreach $key (sort keys %header) {
next if $key =~ /^_/;
# If the key matches the Extensions (Axis or errors)
# then skip. Should probably use the Extensions array
# itself to work out which components to skip
next if $key eq 'Axis';
next if $key eq 'Errors';
next if $key eq 'Extensions';
next if $key eq 'Hist';
# Deal with HASH arrays early
ref($header{"$key"}) eq 'HASH' && do {
$oldhash = $hashname;
$oldtypes = $stypes;
$hashname .= ".$key";
$hashname =~ s/^\.//; # remove leading dot
$stypes .= ".".$header{"$key"}{'_TYPES'}{STRUCT}{"$key"};
$stypes =~ s/^\.//; # remove leading dot
whash(\%{$header{"$key"}}, $outndf, $hashname, $stypes, $status);
$hashname = $oldhash;
$stypes = $oldtypes;
next;
};
$path = $hashname . ".$key";
$path =~ s/^\.//; # remove leading dot
if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
# This is not extension info
ndf_cput($header{$key}, $outndf, $key, $status);
} else {
# Find nested structures
@structures = split(/\./, uc($path));
# Last entry in structure is the important name
$comp = pop(@structures);
# Find list of structures
$structs = join(".", @structures);
$structs = 'PERLDL' unless (defined $structs && $structs =~ /./);
$stypes = 'PERLDL_HDR' unless (defined $stypes && $stypes =~ /./);
$loc = mkstruct($outndf, $structs, $stypes, $status);
undef $stypes;
# Now write the component
# $loc is a locator to the last structure component
$type = $header{"_TYPES"}{"$key"};
# What is the data type?
unless ($type =~ /./) {
# number or character or array or pdl
# All arrays are chars - should be PDLs if nums
ref($header{"$key"}) eq 'ARRAY' && ($type = "_CHAR");
ref($header{"$key"}) eq 'PDL' &&
($type = $pdltypes{$header{"$key"}{"Datatype"}});
ref($header{"$key"}) || do {
if ($header{"$key"} =~ /^(-?)(\d*)(\.?)(\d*)([Ee][-\+]?\d+)?$/) {
$header{"$key"} =~ /\.|[eE]/
&& ($type = "_DOUBLE")
|| ($type = "_INTEGER");
} else {
($type = '_CHAR') # Character
}
};
}
# Create and write component
# PDLs
if (ref($header{"$key"}) eq 'PDL') {
my $pdl = $header{"$key"};
@bounds = $pdl->dims;
my $n = $#bounds + 1;
dat_new($loc, $comp, $type, $n, \@bounds, $status);
$nbytes = PDL::Core::howbig($pdl->get_datatype) *
$pdl->nelem;
cmp_mapv($loc, $comp, $type, 'WRITE', $pntr, $el, $status);
string2mem(${$pdl->get_dataref}, $nbytes, $pntr)
if ($status == &NDF::SAI__OK);
cmp_unmap($loc, $comp, $status);
}
# SCALARS
ref($header{"$key"}) || do {
if ($type =~ /^_CHAR/) {
$length = length($header{"$key"});
dat_new0c($loc, $comp, $length, $status);
cmp_put0c($loc, $comp, $header{"$key"}, $status)
if $length > 0;
} elsif ($type =~ /^_LOG/) {
# In this case, add a check for FALSE or TRUE as strings
dat_new0l($loc, $comp, $status);
my $value = $header{$key};
if ($value eq 'FALSE') {
$value = 0;
} elsif ($value eq 'TRUE') {
$value = 1;
}
cmp_put0l($loc, $comp, $value, $status);
} else {
no strict "refs"; # Turn off strict
$which = lc(substr($type, 1,1));
&{"dat_new0$which"}($loc, $comp, $status);
&{"cmp_put0$which"}($loc, $comp, $header{"$key"}, $status);
}
};
# CHARACTER ARRAYS
ref($header{"$key"}) eq 'ARRAY' && do {
# First go through the array and remove any
# entries that are references. I draw the line at storing
# nested arrays
my @norefs = grep { not ref($_) } @{$header{"$key"}};
($maxsize, $packed) = &NDF::pack1Dchar(\@norefs);
$ndims = 1;
@dim = ($#norefs+1);
dat_newc($loc, $comp, $maxsize, $ndims, @dim, $status);
cmp_putvc($loc, $comp, $dim[0], @norefs, $status);
};
# Annul all the locators
dat_annul($loc, $status);
}
}
}
# Subroutine to make extension structures of arbritrary depth
# Just pass dot separated string to this function
# A locator to the final structure is returned
sub mkstruct {
my ($outndf, $structs, $types, $status) = @_;
my ($extension, @locs, @dim, $ndims, $xtype, $loc, $retloc);
my ($set, $prmry, $there);
my (@structures) = split('\.', $structs);
my (@types) = split('\.', $types);
return &NDF::DAT__NOLOC if $#structures < 0;
return &NDF::DAT__NOLOC if $status != &NDF::SAI__OK;
# Does extension exist
$extension = shift(@structures);
$xtype = shift(@types);
ndf_xstat($outndf, $extension, $there, $status);
# Make it if necessary
unless ($there) {
print "Writing $extension extension ($xtype)...\n" if $PDL::verbose;
@dim = (0); # No pdl extensions!
$ndims = 0;
ndf_xnew($outndf, $extension, $xtype, $ndims, @dim, $loc, $status);
} else { # Get the locator to the extension
ndf_xloc($outndf, $extension, 'UPDATE', $loc, $status);
}
@locs = (); # store the locators so that they can be annulled later
push(@locs, $loc);
# Make the structures, end up with a locator to the correct component
for (0..$#structures) {
dat_there($loc, $structures[$_], $there, $status);
unless ($there) { # Make it if it isnt there
my $type = $types[$_];
@dim = (0);
dat_new($loc, $structures[$_], $type, 0, @dim, $status);
}
dat_find($loc, $structures[$_], $loc, $status);
push(@locs, $loc); # Store the new locator
}
# Clone the locator
dat_clone($loc, $retloc, $status);
$set = 1;
$prmry = 1;
dat_prmry($set, $retloc, $prmry, $status);
# Annul all except the last locator
foreach (@locs) {
dat_annul($_, $status);
}
return &NDF::DAT__NOLOC unless $status == &NDF::SAI__OK;
return $retloc;
}
# end of module
1;
!NO!SUBS!