3
# The contents depend on whether we have bad-value support in PDL.
4
# (there's only a small amount of code that uses bad values)
6
# - this is called by Makefile.PL to ensure the file exists
7
# when the makefile is written
13
use File::Basename qw(&basename &dirname);
15
# check for bad value support
16
use vars qw( $bvalflag $usenan );
18
require File::Spec->catfile( File::Spec->updir, File::Spec->updir, "Basic", "Core", "badsupport.p" );
20
# This forces PL files to create target in same directory as PL file.
21
# This is so that make depend always knows where to find PL derivatives.
24
($file = basename($0)) =~ s/\.PL$//;
26
if ($Config{'osname'} eq 'VMS' or
27
$Config{'osname'} eq 'OS2'); # "case-forgiving"
30
print "Extracting $file (WITH bad value support)\n";
32
print "Extracting $file (NO bad value support)\n";
35
open OUT,">$file" or die "Can't create $file: $!";
38
#########################################################
40
print OUT <<"!WITH!SUBS!";
42
# This file is automatically created by NDF.pm.PL
43
# - bad value support = $bvalflag
46
print OUT <<'!NO!SUBS!';
52
PDL::IO::NDF - PDL Module for reading and writing Starlink
53
N-dimensional data structures as PDLs.
59
$a = PDL->rndf($file);
61
$a = rndf('test_image');
62
$a = rndf('test_image', 1);
65
wndf($a, 'out_image');
67
propndfx($a, 'template', 'out_image');
71
This module adds the ability to read and write Starlink N-dimensional data
72
files as N-dimensional PDLs.
78
@EXPORT_OK = qw/rndf wndf/;
79
%EXPORT_TAGS = (Func=>[@EXPORT_OK]);
81
@ISA = qw( PDL::Exporter );
88
# Starlink data type conversion
89
use vars qw/%pdltypes %startypes $ndf_loaded $VERSION $EXTNAME/;
93
# Set PDL -> Starlink data types
94
%pdltypes = ("$PDL_B" => "_UBYTE", # was "_BYTE"
96
"$PDL_US" => "_UWORD",
97
"$PDL_L" => "_INTEGER",
102
# Set Starlink -> PDL data types
103
%startypes = ("_BYTE" =>$PDL_B,
112
# This is the name I use in the PDL header hash that contains the
113
# NDF extensions so that they can be recreated
115
$EXTNAME = 'NDF_EXT';
119
return if $ndf_loaded++; # a bit pointless increasing each time?
121
croak 'Cannot use NDF library' if $@ ne "";
127
print OUT <<'!NO!SUBS!';
131
# internal functions for checking bad values
133
my $starbad = 0; # flag
135
use vars qw( @starbad );
139
_init_NDF(); # just to be safe
148
croak 'Unable to initialise bad values from NDF library' if $@ ne "";
152
sub starlink_bad_value ($) {
154
return $starbad[$_[0]];
157
# given a piddle, check if the Starlink bad value matches
158
# the bad value used for this type
159
sub compare_bad_values ($) {
161
return ($starbad[$_[0]] == badvalue($_[0]));
168
print OUT <<'!NO!SUBS!';
174
Reads a piddle from a NDF format data file.
178
$pdl = rndf('file.sdf');
179
$pdl = rndf('file.sdf',1);
181
The '.sdf' suffix is optional. The optional second argument turns off
182
automatic quality masking and returns a quality array as well.
184
BUG? C<rndf('bob.sdf',1)> calls ndf_sqmf(1,...), which means that
185
the quality array is turned into bad pixels - ie the I<opposite> of
186
above. Or am I confused?
188
Header information and NDF Extensions are stored in the piddle as a hash
189
which can be retreived with the C<$pdl-E<gt>gethdr> command.
190
Array extensions are stored in the header as follows:
192
$a - the base DATA_ARRAY
194
If C<$hdr = $a-E<gt>gethdr>;
198
%{$hdr} contains all the FITS headers plus:
199
$$hdr{Error} contains the Error/Variance PDL
200
$$hdr{Quality} The quality byte array (if reqeusted)
201
@{$$hdr{Axis}} Is an array of piddles containing the information
203
$$hdr{NDF_EXT} Contains all the NDF extensions
204
$$hdr{Hist} Contains the history information
205
$$hdr{NDF_EXT}{_TYPES}
206
Data types for non-PDL NDF extensions so that
207
wndf can reconstruct a NDF.
209
All extension information is stored in the header hash array.
210
Extension structures are preserved in hashes, so that the PROJ_PARS
211
component of the IRAS.ASTROMETRY extension is stored in
212
C<$$hdr{NDF_EXT}{IRAS}{ASTROMETRY}{'PROJ_PARS'}>. All array structures are
213
stored as arrays in the Hdr: numeric arrays are stored as PDLs,
214
logical and character arrays are stored as plain Perl arrays. FITS
215
arrays are a special case and are expanded as scalars into the header.
217
PDL does not have a signed byte datatype, so any '_BYTE' data
218
is read into a C<byte> (unsigned) piddle and a warning is printed
224
print OUT <<'!NO!SUBS!';
227
If the starlink bad flag is set, then the bad flag on the output
228
piddle is also set. Starlink bad values are converted to the
229
current bad value used by the piddle type (if they are different).
234
print OUT <<'!NO!SUBS!';
237
# This is one form of the new command
239
sub rndf {PDL->rndf(@_)}
241
# And this is the real form
242
# Allows the command to be called in OO form or as a function
243
sub PDL::rndf { # Read a piddle from a NDF file
246
barf 'Usage: $a = rndf($file,[$nomask]); $a = PDL->rndf(...) '
247
if ($#_!=0 && $#_!=1);
249
# ensure we have the NDF module
252
# Setup the Header Hash
255
# Read in the filename
256
# And setup the new PDL
258
my $pdl = $class->new;
259
my $nomask = shift if $#_ > -1;
261
my ($infile, $status, $indf, $entry, @info, $value);
263
# Strip trailing .sdf if one is present
264
# remove the trailing extension names
265
# (anything after a `.' UNLESS it's part of a directory
266
# name, hence the check for ^/)
269
$infile =~ s!\.[^/]*$!!;
271
# If file is not there
272
croak "Cannot find $infile.sdf" unless -e "$infile.sdf";
275
$status = &NDF::SAI__OK;
277
# Begin NDF and ERR session
282
ndf_find(&NDF::DAT__ROOT, $file, $indf, $status);
284
# unset automatic quality masking if $nomask is defined
285
# I think the logic here is wrong (ie 0 and 1 should be swapped)
286
$nomask = 0 unless (defined $nomask);
287
$nomask = 1 if $nomask != 0;
288
ndf_sqmf($nomask, $indf, $status);
290
# Read the data array (need to pass the header as well)
292
( $status, $badflag ) = rdata($indf, $pdl, $nomask, $header, $class, $status);
295
raxes($indf, $pdl, $header, $class, $status);
298
rhdr($indf, $pdl, $header, $class, $status);
300
# Read history information
301
rhist($indf, $pdl, $header, $status);
304
@info = ('Label', 'Title', 'Units');
307
ndf_cget($indf, $entry, $value, $status);
308
$$header{"$entry"} = $value if $value ne 'NULL';
309
print "$entry is $value\n" if ($PDL::verbose && $value ne 'NULL');
314
ndf_annul($indf, $status);
317
if ($status != &NDF::SAI__OK) {
318
err_flush($status); # Flush the error system
320
croak("rndf: obtained NDF error");
325
# Put the header into the main pdl
326
$pdl->sethdr($header);
331
print OUT <<'!NO!SUBS!';
332
# set the bad flag if it was set in the input file
333
print "Bad flag status is: $badflag\n" if $PDL::verbose;
334
$pdl->badflag(1) if $badflag;
340
print OUT <<'!NO!SUBS!';
348
Writes a piddle to a NDF format file:
355
wndf can be used for writing PDLs to NDF files.
356
The '.sdf' suffix is optional. All the extensions
357
created by rndf are supported by wndf. This means that error, axis
358
and quality arrays will be written if they exist. Extensions are also
359
reconstructed by using their name (ie FIGARO.TEST would be expanded as
360
a FIGARO extension and a TEST component). Hdr keywords Label, Title
361
and Units are treated as special cases and are written to the label,
362
title and units fields of the NDF.
364
Header information is written to corresponding NDF extensions.
365
NDF extensions can also be created in the {NDF} hash by using a key
366
containing '.', ie {NDF}{'IRAS.DATA'} would write the information to
367
an IRAS.DATA extension in the NDF. rndf stores this as
368
C<$$hdr{NDF}{IRAS}{DATA}> and the two systems are interchangeable.
370
rndf stores type information in {NDF}{'_TYPES'} and below so that
371
wndf can reconstruct the data type of non-PDL extensions. If no
372
entry exists in _TYPES, wndf chooses between characters, integer and
373
double on a best guess basis. Any perl arrays are written as CHAR
374
array extensions (on the assumption that numeric arrays will exist as
380
print OUT <<'!NO!SUBS!';
383
The bad flag is written to the file, if set.
389
print OUT<<'!NO!SUBS!';
392
# This is one form of the new command
397
# And this is the real form
398
# Allows the command to be called in OO form or as a function
399
sub PDL::wndf { # Write a PDL to an NDF format file
401
barf 'Usage: wndf($pdl,$file)' if $#_!=1;
403
# ensure we have the NDF module
406
my ($indf, $place, $status, $outndf);
409
my ($pdl, $outfile) = @_;
411
# Check that we are dealing with a PDL
412
croak 'Argument is not a PDL variable' unless $pdl->isa('PDL');
415
$status = &NDF::SAI__OK;
417
# Strip trailing .sdf if one is present
418
$outfile =~ s/\.sdf$//;
425
ndf_place(&NDF::DAT__ROOT(), $outfile, $place, $status);
427
# Need to create data_array component
428
# Change the bounds to match the PDL
432
@lbnd = map { 1 } (0..$#ubnd);
434
ndf_new($pdltypes{$pdl->get_datatype}, $#ubnd+1, \@lbnd, \@ubnd, $place,
437
# Set the application name
438
ndf_happn('PDL::wndf', $status);
440
# Write the data to this file
441
wdata($outndf, $pdl, $status);
443
# Write the header and title
444
whdr($outndf, $pdl, $status);
446
# Create a history extension
447
ndf_hcre($outndf, $status);
449
# Annul the identifier
450
ndf_annul($outndf, $status);
455
if ($status != &NDF::SAI__OK) {
457
croak "Bad STATUS in wndf\n";
467
Routine to write a PDL to an NDF by copying the extension information
468
from an existing NDF and writing DATA,VARIANCE, QUALITY and AXIS info
469
from a PDL (if they exist).
471
Extensions, labels and history are propogated from the old NDF.
472
No new extension information is written.
474
This command has been superseded by L<wndf()|/wndf()>.
478
*propndfx = \&PDL::propndfx;
480
sub PDL::propndfx { # Write a PDL to a NDF format file
482
barf 'Usage: propndfx($pdl, $infile, $outfile)' if $#_!=2;
484
# ensure we have the NDF module
487
my ($indf, $in_place, $status, $outplace, $outndf);
489
my ($pdl, $infile, $outfile) = @_;
492
$status = &NDF::SAI__OK;
494
# Strip trailing .sdf if one is present
495
$outfile =~ s/\.sdf$//;
496
$infile =~ s/\.sdf$//;
497
# File is the first thing before a .
498
my $file = (split(/\./,$infile))[0];
500
croak "$file does not exist\n" unless -e "$file.sdf";
507
ndf_find(&NDF::DAT__ROOT(), $infile, $indf, $status);
509
# Set the application name
510
ndf_happn('PDL::propndfx', $status);
513
ndf_place(&NDF::DAT__ROOT(), $outfile, $outplace, $status);
515
# Copy the extension information to outfile
516
ndf_scopy($indf, ' ', $outplace, $outndf, $status);
518
# Annul the input file identifier
519
ndf_annul($indf, $status);
521
# Write the data to this file
522
wdata($outndf, $pdl, $status);
524
# Annul the identifier
525
ndf_annul($outndf, $status);
530
if ($status != &NDF::SAI__OK) {
532
croak "Bad STATUS in propndfx\n";
542
The perl NDF module must be available. This is available from the
543
author or from Starlink (http://www.starlink.rl.ac.uk).
545
If an NDF is read which contains AST World Coordinate information
546
(a .WCS component) this information is currently ignored. Currently
547
WCS information can only be written and stored using standard FITS headers.
548
See http://rlspc5.bnsc.rl.ac.uk/star/docs/sun211.htx/sun211.html#xref_
549
for more information on AST.
553
This module was written by Tim Jenness E<lt>t.jenness@jach.hawaii.eduE<gt>.
554
Copyright (C) Tim Jenness 1997-2000. All Rights Reserved.
558
L<PDL::FAQ> for general information on the Perl Data language,
559
L<NDF> for information on the NDF module.
564
#######################################################################
566
# These are the generic routines used by the rndf command
569
# Read Data, Quality and Variance
571
# returns ( $status, $badflag ), where $badflag is:
573
# 1 - the piddle may contain bad data
574
# 0 - the piddle contains no bad data
576
# This reports starlink's bad-data flag (ie it's valid
577
# even in bad value support is not available in PDL)
580
my ($indf, $pdl, $nomask, $header, $class, $status) = @_;
582
return ( $status, undef ) if $status != &NDF::SAI__OK;
584
my ($maxdims, $ndim, @dim, @comps, $dcomp, $tcomp, $exist);
585
my ($type, $data_pntr, $el, $temppdl, $nbytes, $badbit, $dref);
587
####################################################################
589
####################################################################
591
# Find the dimensions of the data array
592
$maxdims = 100; # Maximum number of dimensions allowed
594
&ndf_dim($indf, $maxdims, \@dim, $ndim, $status);
595
@dim = @dim[0..$ndim-1];
597
print "Dims = @dim\n" if $PDL::verbose;
599
$$header{Extensions} = [];
601
####################################################################
602
# D A T A + V A R I A N C E + Q U A L I T Y #
603
####################################################################
607
@comps = ('Data','Error');
608
push(@comps, 'Quality') if $nomask;
610
for $dcomp (@comps) {
613
$tcomp = 'Variance' if $tcomp eq 'Error';
615
ndf_state($indf, $tcomp, $exist, $status);
617
if ($exist && ($status == &NDF::SAI__OK)) {
619
# Find the type of the data array
620
ndf_type($indf, $tcomp, $type, $status);
623
ndf_map($indf, $dcomp, $type, 'READ', $data_pntr, $el, $status);
625
# Check the bad-data status of this component
626
my $this_badflag = 0;
627
ndf_bad($indf, $dcomp, 0, $this_badflag, $status );
629
if ($status == &NDF::SAI__OK) {
631
print "Reading $dcomp array...\n" if $PDL::verbose;
633
# combine the bad data flag with the "global" one for
635
$badflag |= ($this_badflag != 0);
637
# Create a temporary PDL to deal with separate components
638
if ($dcomp eq 'Data') {
641
$dcomp = 'Errors' if $dcomp eq 'Error';
642
$$header{"$dcomp"} = $class->new;
643
$temppdl = $$header{"$dcomp"};
647
$temppdl->set_datatype($startypes{$type});
648
$temppdl->setdims([@dim]);
649
$dref = $temppdl->get_dataref;
651
# warning about data loss
652
if ( $type eq "_BYTE" ) {
653
print "Warning: NDF type (signed byte) converted to PDL's unsigned byte.\n";
656
# How many bytes in this data type?
658
$nbytes = &NDF::byte_size($type);
660
mem2string($data_pntr, $nbytes*$el, $$dref);
662
push(@{$$header{Extensions}}, $dcomp) if $dcomp ne 'Data';
664
if ($dcomp eq 'Quality') {
665
# Quality bad bit mask
666
ndf_bb($indf, $badbit, $status);
669
$$qhdr{BADBIT} = $badbit;
670
$temppdl->sethdr($qhdr);
676
print OUT <<'!NO!SUBS!';
677
# if badflag is true, convert any STARLINK bad values
680
!compare_bad_values($temppdl->get_datatype) ) {
681
my $star_bad = starlink_bad_value($temppdl->get_datatype);
682
print "Converting from STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
683
$temppdl->inplace->setvaltobad( $star_bad );
690
print OUT <<'!NO!SUBS!';
691
# Flush and Clear temporary PDL
692
$temppdl->upd_data();
696
ndf_unmap($indf, $tcomp, $status);
699
err_flush($status) if $status != &NDF::SAI__OK;
701
return ( $status, $badflag );
706
####################################################################
708
####################################################################
710
# Axes are stored as follows:
711
# Array of PDL axes as an array in header of main pdl (called 'Axis')
712
# Header of this PDL contains label and errors and units etc.
715
my ($indf, $pdl, $header, $class, $status) = @_;
717
return $status if $status != &NDF::SAI__OK;
719
my ($there, $naxis, @dims, $axcomp, $exist, $axtype, $axpntr, $el);
720
my ($nbytes, $temppdl, $tcomp, $daxref, $dref);
723
# Read in axis information
724
ndf_state($indf, 'Axis', $there, $status);
726
$ndims = $pdl->getndims; # Find number of dimensions
730
# Label from 0..$#dims to match array index handling
731
# Want to put the axes into an array of axes
732
# And update the extensions field
734
push(@{$$header{Extensions}}, "Axis");
736
for $naxis (0..$ndims-1) {
737
# Does a CENTRE component exist
739
ndf_astat($indf, $axcomp, $naxis+1, $exist, $status);
741
if ($exist && ($status == &NDF::SAI__OK)) {
743
print "Reading Axis $naxis...\n" if $PDL::verbose;
744
ndf_atype($indf, $axcomp, $naxis+1, $axtype, $status);
746
ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
749
# Set up new PDL for axis info if map was okay
750
if ($status == &NDF::SAI__OK) {
751
print "Number of elements: $el\n" if $PDL::verbose;
752
$$header{Axis}[$naxis] = $class->new;
753
$temppdl = $$header{Axis}[$naxis];
754
$temppdl->set_datatype($startypes{$axtype});
755
$temppdl->setdims([$el]);
756
$dref = $temppdl->get_dataref;
760
$axtype =~ s/^_|^_U//;
761
$nbytes = &NDF::byte_size($axtype);
763
mem2string($axpntr, $nbytes*$el, $$dref);
764
$temppdl->upd_data();
766
# Header info for axis
767
$axhdr = {}; # somewhere to put the pdl
769
# Read VARIANCE array if there
772
ndf_astat($indf, $tcomp, $naxis+1, $exist, $status);
774
if ($exist && ($status == &NDF::SAI__OK)) {
775
print "Reading Axis $naxis errors...\n" if $PDL::verbose;
776
ndf_atype($indf, $tcomp, $naxis+1, $axtype, $status);
778
ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
781
# Set up new PDL for axis info if map was okay
782
if ($status == &NDF::SAI__OK) {
784
$$axhdr{Errors} = $class->new;
785
$$axhdr{Extensions} = ["Errors"];
787
$$axhdr{Errors}->set_datatype($startypes{$axtype});
788
$$axhdr{Errors}->setdims([$el]);
789
$daxref = $$axhdr{Errors}->get_dataref;
792
$axtype =~ s/^_|^_U//;
793
$nbytes = &NDF::byte_size($axtype);
795
mem2string($axpntr, $nbytes*$el, $$daxref);
796
$$axhdr{Errors}->upd_data();
801
# Get label and units
802
for my $entry ('Units', 'Label') {
803
ndf_astat($indf, $entry, $naxis+1, $exist, $status);
806
ndf_acget($indf,$entry, $naxis+1, $value, $status);
807
$axhdr->{"$entry"} = $value;
811
ndf_aunmp($indf,'*',$naxis+1,$status);
813
# Now store this header into temppdl
814
$temppdl->sethdr($axhdr);
817
err_annul($status) unless $status == &NDF::SAI__OK;
828
####################################################################
829
# E X T E N S I O N S #
830
####################################################################
833
my ($indf, $pdl, $header, $class, $status) = @_;
835
return $status if $status != &NDF::SAI__OK;
837
my ($nextn, $xname, $xloc, $extn, $size, @fits, $entry, $value);
838
my ($comment, $item, $el);
840
# Read in any header information from the extensions
842
ndf_xnumb($indf, $nextn, $status);
843
print "Reading in header information...\n"
844
if (($nextn > 0) && $PDL::verbose);
846
# Read in extensions one at a time
847
for $extn (1..$nextn) {
850
ndf_xname($indf, $extn, $xname, $status);
853
ndf_xloc($indf, $xname, 'READ', $xloc, $status);
855
# FITS is a special case and must be expanded
856
if ($xname eq 'FITS') {
858
dat_size($xloc, $size, $status);
859
&dat_getvc($xloc, $size, \@fits, $el, $status);
861
if ($status == &NDF::SAI__OK) {
862
print "Reading FITS header...\n" if $PDL::verbose;
863
for $entry (0..$el-1) {
864
($item, $value, $comment) = &fits_get_nth_item(\@fits,$entry);
865
$$header{$item} = $value;
866
$$header{'_COMMENTS'}{$item} = $comment;
869
print "Error mapping $xname\n";
875
# Read in extensions to $EXTNAME
876
$status = &find_prim($xloc, \%{$$header{$EXTNAME}}, $class, $status);
879
dat_annul($xloc, $status);
884
####################################################################
886
####################################################################
889
my ($indf, $pdl, $header, $status) = @_;
891
return $status if $status != &NDF::SAI__OK;
893
my ($nrec, $app, $date, $i);
895
# Read HISTORY information into a Hist header
897
ndf_hnrec($indf,$nrec,$status);
899
if ($status == &NDF::SAI__OK && ($nrec > 0)) {
900
print "Reading history information into 'Hist'...\n" if $PDL::verbose;
902
$$header{Hist}{'Nrecords'} = $nrec;
904
# Loop through the history components and find last "APPLICATION"
907
ndf_hinfo($indf, 'APPLICATION', $i, $app, $status);
908
ndf_hinfo($indf, 'DATE', $i, $date, $status);
910
$$header{Hist}{"Application$i"} = $app;
911
$$header{Hist}{"Date$i"} = $date;
913
print " $app on $date\n" if $PDL::verbose;
923
################ Generic routines ###########################
924
#################### INTERNAL ###############################
926
# Find primitive components below a given HDS locator
927
# FITS information is stored in Hdr
928
# NDF extensions stored in NDF
932
my ($loc, $href, $class, $status) = @_;
933
my ($prim, $type, $size, @dim, $ndims, $struct, $name, $nloc, $el);
934
my ($value, @values, $item, $entry, $maxdims);
935
my ($packtype, $ncomp, $packed, $comp);
936
my ($pntr, $nbytes, $comment, $dref);
938
# Return if bad status
939
return 0 if $status != &NDF::SAI__OK;
944
# Is it a primitive component
945
dat_prim($loc, $prim, $status);
947
# Type, size and shape
948
dat_type($loc, $type, $status);
949
dat_size($loc, $size, $status);
950
dat_shape($loc, $maxdims, \@dim, $ndims, $status);
951
dat_name($loc, $name, $status);
954
print "\tReading component: $name ($type) " if $PDL::verbose;
956
# Store type of information
957
$$href{"_TYPES"}{"$name"} = $type;
959
# Read it into the header
960
if (($size == 1) && ($ndims == 0)) {
961
if ($type eq '_DOUBLE') {
962
# dat_get0c doesn't return full precision for doubles
963
dat_get0d($loc, $value, $status);
964
if ($status != &NDF::SAI__OK) {
965
print "Error mapping scalar 0d $type $name\n";
970
dat_get0c($loc, $value, $status);
971
if ($status != &NDF::SAI__OK) {
972
print "Error mapping scalar 0c $name\n";
978
$$href{"$name"} = $value if ($status == &NDF::SAI__OK);
982
@dim = @dim[0..$ndims-1];
984
# Inform the user of dimensions
985
print "[@dim]" if $PDL::verbose;
987
# Map arrays (don't need to unpack them)
988
if ($type =~ /CHAR|LOG/) {
991
dat_getvc($loc, $size, \@values, $el, $status);
993
if ($status != &NDF::SAI__OK) {
994
print "Error mapping $name\n";
996
$$href{"$name"} = [@values];
1003
dat_mapv($loc, $type, 'READ', $pntr, $el, $status);
1005
if ($status != &NDF::SAI__OK) {
1006
print "Error mapping $name\n";
1010
if ($status == &NDF::SAI__OK) {
1012
$$href{"$name"} = $class->new;
1013
$$href{"$name"}->set_datatype($startypes{$type});
1014
$$href{"$name"}->setdims([@dim]);
1015
$dref = $$href{"$name"}->get_dataref();
1017
# How many bytes in this data type?
1018
$type =~ s/^_|^_U//;
1019
$nbytes = &NDF::byte_size($type);
1021
mem2string($pntr, $nbytes*$el, $$dref);
1022
$$href{"$name"}->upd_data();
1025
dat_unmap($loc,$status);
1028
print "\n" if $PDL::verbose;
1032
my ($ndimx, @dim, $ndim);
1033
dat_shape($loc, 20, \@dim, $ndim, $status);
1035
# If we have a single dimension ($ndim=0) then
1036
# proceed. Cant yet cope with multi-dimensional
1039
dat_ncomp($loc, $ncomp, $status);
1040
print "$name ($type) has $ncomp components\n" if $PDL::verbose;
1042
# Store type of extension
1043
$$href{$name}{"_TYPES"}{STRUCT}{"$name"} = $type;
1045
for $comp (1..$ncomp) {
1046
dat_index($loc, $comp, $nloc, $status);
1047
$status = &find_prim($nloc, \%{$$href{$name}}, $class, $status);
1048
dat_annul($nloc, $status);
1051
print "Extension $name is an array structure -- skipping\n"
1060
###### Routines for WRITING data ######################################
1064
my ($outndf, $pdl, $status) = @_;
1065
my (@bounds, $ndims, @lower, @comps, $dcomp, $temppdl, $type);
1066
my ($pntr, $el, $nbytes, $axis, $axcomp, $axpntr, %hdr, %axhdr);
1067
my ($entry, $value, $i, $axndims, $axtype, $tcomp, @acomps);
1070
# Return if bad status
1071
return 0 if $status != &NDF::SAI__OK;
1073
####################################################################
1075
####################################################################
1076
# Change the bounds to match the PDL
1077
@bounds = $pdl->dims;
1079
@lower = map { 1 } (0..$ndims);
1081
&ndf_sbnd($ndims+1, \@lower, \@bounds, $outndf, $status);
1084
####################################################################
1085
# D A T A + V A R I A N C E + Q U A L I T Y #
1086
####################################################################
1087
# Map data, variance, quality...
1088
@comps = ('Data','Errors','Quality');
1090
# Retrieve header and check that the header is a hash ref
1091
my $hdr = $pdl->gethdr;
1092
if (ref($hdr) eq 'HASH') {
1098
for $dcomp (@comps) {
1100
if ($dcomp eq 'Data') {
1103
if (exists $hdr{"$dcomp"}) {
1104
$temppdl = $hdr{"$dcomp"};
1105
# Check that we have a PDL
1106
next unless UNIVERSAL::isa($temppdl, 'PDL');
1108
next; # Skip this component
1112
# Check that we have some data
1113
if ($temppdl->nelem > 0) {
1115
if (($dcomp eq 'Quality') && ($temppdl->get_datatype != $PDL_B)) {
1116
# Quality must be bytes so convert to BYTE if this is not true
1117
$temppdl = byte ($$pdl{"$dcomp"});
1120
$type = $pdltypes{$temppdl->get_datatype};
1121
$type = '_UBYTE' if $dcomp eq 'Quality'; # No choice for QUALITY
1123
# Set the output type of the NDF
1125
$tcomp = 'Variance' if $tcomp eq 'Errors';
1126
ndf_stype($type, $outndf, $tcomp, $status) unless $dcomp eq 'Quality';
1128
print "Mapping $dcomp , $type\n" if $PDL::verbose;
1130
$dcomp = 'Error' if $dcomp eq 'Errors';
1131
ndf_map($outndf, $dcomp, $type, 'WRITE', $pntr, $el, $status);
1133
if ($status == &NDF::SAI__OK) {
1134
# Number of bytes per entry
1135
$nbytes = PDL::Core::howbig($temppdl->get_datatype);
1137
# Convert to 1D data stream
1138
my $p1d = $temppdl->clump(-1);
1139
# Total number of bytes
1140
$nbytes *= $p1d->getdim(0);
1145
print OUT <<'!NO!SUBS!';
1146
# if badflag is true, convert any PDL bad values to STARLINK ones
1147
# note: we have to create a copy of the data to ensure that the
1148
# bad data value of the original piddle does not get changed.
1149
# per-piddle bad values would obviate the need for this.
1150
if ( $temppdl->badflag() ) {
1151
print "Setting bad flag for component $dcomp ...\n" if $PDL::verbose;
1152
ndf_sbad( 1, $outndf, $dcomp, $status );
1153
if ( !compare_bad_values($temppdl->get_datatype) ) {
1154
my $star_bad = starlink_bad_value($temppdl->get_datatype);
1155
print "Converting from PDL to STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
1156
$temppdl = $temppdl->setbadtoval( $star_bad ); # do NOT do inplace
1164
print OUT <<'!NO!SUBS!';
1166
string2mem(${$temppdl->get_dataref}, $nbytes, $pntr);
1169
if ($dcomp eq 'Quality') {
1170
ndf_sbb($$temppdl{Hdr}{BADBIT}, $outndf, $status)
1171
if defined $$temppdl{Hdr}{BADBIT};
1175
ndf_unmap($outndf, $tcomp, $status);
1179
# free the temporary pdl
1185
####################################################################
1187
####################################################################
1188
# A X I S information is expected as an array in the header
1189
# called 'Axis'. These PDLs contain the CENTRE data and any further
1190
# info is stored in their header.
1191
# @bounds is accessed to find expected shape
1194
# Simply look in the header for a Axis
1195
if (exists $hdr{Axis}) {
1196
# Check that we have an array
1197
if (ref($hdr{Axis}) eq 'ARRAY') {
1198
# Now loop over axes
1199
for my $i (0..$#{$hdr{Axis}} ) {
1200
# Loop unless status is bad
1201
last unless $status == &NDF::SAI__OK;
1203
# Extract the ith axis PDL from the array
1204
my $axis = $hdr{Axis}->[$i];
1207
if (UNIVERSAL::isa($axis, 'PDL')) {
1208
# We now want to copy the data and if necessary the
1209
# Error array. Since there are only two I will do it the
1210
# long way by explcitly storing data and then error
1213
$axtype = $pdltypes{$axis->get_datatype};
1214
ndf_astyp($axtype, $outndf, 'Centre', $i+1, $status);
1216
# Okay we can now map this pdl
1217
ndf_amap($outndf, 'Centre', $i+1, $axtype, 'WRITE', $axpntr,
1220
# Check that we have the correct number of entries
1221
my $nelem = $axis->nelem;
1222
if ($el == $nelem) {
1223
print "Mapping axis " , $i+1 , "\n" if $PDL::verbose;
1225
# Number of bytes per entry
1226
$nbytes = PDL::Core::howbig($axis->get_datatype) * $el;
1229
string2mem( $ { $axis->get_dataref }, $nbytes, $axpntr)
1230
if ($status == &NDF::SAI__OK);
1234
" is the wrong size ($el values required but got ".
1235
$nelem . ")- ignoring";
1238
ndf_aunmp($outndf, '*', $i+1, $status);
1241
my $axhdr = $axis->gethdr; # Retrieve and check header
1242
if (ref($axhdr) eq 'HASH') {
1248
# Look for an Errors component in the header hash
1249
if (exists $axhdr{Errors}) {
1250
my $axerr = $axhdr{Errors};
1251
if (UNIVERSAL::isa($axerr, 'PDL')) {
1254
$axtype = $pdltypes{$axerr->get_datatype};
1255
ndf_astyp($axtype, $outndf, 'Variance', $i+1, $status);
1257
# Okay we can now map this pdl
1258
ndf_amap($outndf, 'Errors', $i+1, $axtype, 'WRITE',
1259
$axpntr, $el, $status);
1261
# Check that we have the correct number of entries
1262
my $nelem = $axerr->nelem;
1263
print "Nelem: $nelem and $el\n";
1264
if ($nelem == $el) {
1265
print "Mapping errors for axis " . $i+1 . "\n"
1268
# Number of bytes per entry
1269
$nbytes = PDL::Core::howbig($axerr->get_datatype) * $el;
1272
string2mem($ {$axerr->get_dataref}, $nbytes, $axpntr)
1273
if ($status == &NDF::SAI__OK);
1276
carp "Error PDL for Axis ",$i+1,
1277
" is the wrong size ($el values required but got ".
1278
$axerr->nelem . ")- ignoring";
1281
ndf_aunmp($outndf, '*', $i+1, $status);
1287
# Now character components
1288
foreach my $char (keys %axhdr) {
1289
if ($char =~ /LABEL|UNITS/i) {
1290
$value = $axhdr{"$char"};
1291
ndf_acput($value, $outndf, $char, $i+1, $status)
1292
if length($value) > 0;
1308
# Write header information to NDF extension
1309
# This routine does not reconstruct an NDF identically to one that
1310
# was read in. FITS extensions are made. All non-integer numeric types
1311
# are written as doubles (apart from PDLs which carry there own type
1312
# information) unless _TYPES information exists in {NDF}.
1317
my ($outndf, $pdl, $status) = @_;
1319
my (%header, @fitsdim, $fitsloc, $value);
1320
my (%unused, @fits, $hashref, $hdr);
1322
# Return if bad status
1323
return 0 if $status != &NDF::SAI__OK;
1325
print "Writing header information...\n" if $PDL::verbose;
1327
# Write FITS header from {Hdr}
1329
# Retrieve and check header from PDL
1330
$hdr = $pdl->gethdr;
1331
if (ref($hdr) eq 'HASH') {
1340
foreach my $key (sort keys %header) {
1342
next if $key eq '_COMMENTS';
1344
if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
1345
# This is not extension info
1346
ndf_cput($header{$key}, $outndf, $key, $status)
1347
if length($header{$key}) > 0;
1350
# Only write scalars
1351
if (not ref($header{"$key"})) {
1353
push(@fits, fits_construct_string($key,$header{$key},
1354
$header{'_COMMENTS'}{$key}) );
1358
$unused{$key} = $header{$key};
1363
# Write the FITS extension
1367
@fitsdim = ($#fits+1);
1368
ndf_xnew($outndf, 'FITS', '_CHAR*80', 1, @fitsdim, $fitsloc, $status);
1369
dat_put1c($fitsloc, $#fits+1, \@fits, $status);
1370
dat_annul($fitsloc, $status);
1373
# Loop through all NDF header info and array Hdr information
1374
# Note that %unused still contains our NDF hash so
1375
# we must remove that before proceeding (else we end up with
1376
# an 'NDF_EXT' extension!
1378
delete $unused{$EXTNAME};
1380
foreach $hashref (\%unused, $header{$EXTNAME}) {
1381
whash($hashref, $outndf, '', '', $status);
1386
#-----------------------------------------------------
1387
# This routine writes a hash array to a NDF extension
1388
# Keys that have leading underscores are skipped
1389
# $hashname contains the extension name to use by default.
1390
# This is extended by any '.' in the key.
1394
my ($hash, $outndf, $hashname, $stypes, $status) = @_;
1396
my ($key, $comp, @structures, $extension, $xtype, @dim, $ndims, $loc);
1397
my (@locs, $there, $type, %header, $packed, $length);
1398
my (@bounds, $nbytes, $pntr, $maxsize, $which, $el, $path);
1399
my ($oldhash, $oldtypes, $structs);
1401
if (defined $hash) {
1408
my $root = 1 if length($hashname) == 0; # Mark a ROOT structure
1410
# Return if bad status
1411
return 0 if $status != &NDF::SAI__OK;
1413
foreach $key (sort keys %header) {
1415
next if $key =~ /^_/;
1417
# If the key matches the Extensions (Axis or errors)
1418
# then skip. Should probably use the Extensions array
1419
# itself to work out which components to skip
1421
next if $key eq 'Axis';
1422
next if $key eq 'Errors';
1423
next if $key eq 'Extensions';
1424
next if $key eq 'Hist';
1426
# Deal with HASH arrays early
1427
ref($header{"$key"}) eq 'HASH' && do {
1428
$oldhash = $hashname;
1429
$oldtypes = $stypes;
1430
$hashname .= ".$key";
1431
$hashname =~ s/^\.//; # remove leading dot
1432
$stypes .= ".".$header{"$key"}{'_TYPES'}{STRUCT}{"$key"};
1433
$stypes =~ s/^\.//; # remove leading dot
1434
whash(\%{$header{"$key"}}, $outndf, $hashname, $stypes, $status);
1435
$hashname = $oldhash;
1436
$stypes = $oldtypes;
1442
$path = $hashname . ".$key";
1443
$path =~ s/^\.//; # remove leading dot
1445
if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
1446
# This is not extension info
1447
ndf_cput($header{$key}, $outndf, $key, $status);
1450
# Find nested structures
1451
@structures = split(/\./, uc($path));
1453
# Last entry in structure is the important name
1454
$comp = pop(@structures);
1456
# Find list of structures
1457
$structs = join(".", @structures);
1458
$structs = 'PERLDL' unless (defined $structs && $structs =~ /./);
1459
$stypes = 'PERLDL_HDR' unless (defined $stypes && $stypes =~ /./);
1461
$loc = mkstruct($outndf, $structs, $stypes, $status);
1464
# Now write the component
1465
# $loc is a locator to the last structure component
1466
$type = $header{"_TYPES"}{"$key"};
1468
# What is the data type?
1469
unless ($type =~ /./) {
1470
# number or character or array or pdl
1471
# All arrays are chars - should be PDLs if nums
1472
ref($header{"$key"}) eq 'ARRAY' && ($type = "_CHAR");
1474
ref($header{"$key"}) eq 'PDL' &&
1475
($type = $pdltypes{$header{"$key"}{"Datatype"}});
1477
ref($header{"$key"}) || do {
1478
if ($header{"$key"} =~ /^(-?)(\d*)(\.?)(\d*)([Ee][-\+]?\d+)?$/) {
1479
$header{"$key"} =~ /\.|[eE]/
1480
&& ($type = "_DOUBLE")
1481
|| ($type = "_INTEGER");
1483
($type = '_CHAR') # Character
1489
# Create and write component
1492
if (ref($header{"$key"}) eq 'PDL') {
1493
my $pdl = $header{"$key"};
1494
@bounds = $pdl->dims;
1495
my $n = $#bounds + 1;
1496
dat_new($loc, $comp, $type, $n, \@bounds, $status);
1497
$nbytes = PDL::Core::howbig($pdl->get_datatype) *
1499
cmp_mapv($loc, $comp, $type, 'WRITE', $pntr, $el, $status);
1500
string2mem(${$pdl->get_dataref}, $nbytes, $pntr)
1501
if ($status == &NDF::SAI__OK);
1502
cmp_unmap($loc, $comp, $status);
1506
ref($header{"$key"}) || do {
1508
if ($type =~ /^_CHAR/) {
1509
$length = length($header{"$key"});
1510
dat_new0c($loc, $comp, $length, $status);
1511
cmp_put0c($loc, $comp, $header{"$key"}, $status)
1513
} elsif ($type =~ /^_LOG/) {
1514
# In this case, add a check for FALSE or TRUE as strings
1515
dat_new0l($loc, $comp, $status);
1516
my $value = $header{$key};
1517
if ($value eq 'FALSE') {
1519
} elsif ($value eq 'TRUE') {
1522
cmp_put0l($loc, $comp, $value, $status);
1524
no strict "refs"; # Turn off strict
1525
$which = lc(substr($type, 1,1));
1526
&{"dat_new0$which"}($loc, $comp, $status);
1527
&{"cmp_put0$which"}($loc, $comp, $header{"$key"}, $status);
1532
ref($header{"$key"}) eq 'ARRAY' && do {
1533
# First go through the array and remove any
1534
# entries that are references. I draw the line at storing
1536
my @norefs = grep { not ref($_) } @{$header{"$key"}};
1538
($maxsize, $packed) = &NDF::pack1Dchar(\@norefs);
1540
@dim = ($#norefs+1);
1541
dat_newc($loc, $comp, $maxsize, $ndims, @dim, $status);
1542
cmp_putvc($loc, $comp, $dim[0], @norefs, $status);
1545
# Annul all the locators
1546
dat_annul($loc, $status);
1553
# Subroutine to make extension structures of arbritrary depth
1554
# Just pass dot separated string to this function
1555
# A locator to the final structure is returned
1559
my ($outndf, $structs, $types, $status) = @_;
1560
my ($extension, @locs, @dim, $ndims, $xtype, $loc, $retloc);
1561
my ($set, $prmry, $there);
1563
my (@structures) = split('\.', $structs);
1564
my (@types) = split('\.', $types);
1566
return &NDF::DAT__NOLOC if $#structures < 0;
1567
return &NDF::DAT__NOLOC if $status != &NDF::SAI__OK;
1569
# Does extension exist
1571
$extension = shift(@structures);
1572
$xtype = shift(@types);
1574
ndf_xstat($outndf, $extension, $there, $status);
1576
# Make it if necessary
1578
print "Writing $extension extension ($xtype)...\n" if $PDL::verbose;
1579
@dim = (0); # No pdl extensions!
1581
ndf_xnew($outndf, $extension, $xtype, $ndims, @dim, $loc, $status);
1582
} else { # Get the locator to the extension
1583
ndf_xloc($outndf, $extension, 'UPDATE', $loc, $status);
1587
@locs = (); # store the locators so that they can be annulled later
1590
# Make the structures, end up with a locator to the correct component
1591
for (0..$#structures) {
1592
dat_there($loc, $structures[$_], $there, $status);
1593
unless ($there) { # Make it if it isnt there
1594
my $type = $types[$_];
1596
dat_new($loc, $structures[$_], $type, 0, @dim, $status);
1598
dat_find($loc, $structures[$_], $loc, $status);
1599
push(@locs, $loc); # Store the new locator
1603
dat_clone($loc, $retloc, $status);
1606
dat_prmry($set, $retloc, $prmry, $status);
1608
# Annul all except the last locator
1610
dat_annul($_, $status);
1613
return &NDF::DAT__NOLOC unless $status == &NDF::SAI__OK;
added
removed
IO/NDF/NDF.pm.PL
