package GoMOOS::NetCDF; =head1 NAME GoMOOS::NetCDF - a library of routines for handling NetCDF files =head1 SYNOPSIS Wraps bas ic NetCDF functions for use raw by scripts, or by other classes (e.g. GoM OOS::Buoys::NetCDF). Note that many of the functions in this library may have a tendency to die suddenly. We try to check them, and we mark them if possible. Note that this is primarily a record-oriented library. =cut my $version = '$Revision: 1.29 $'; use strict; use Carp 'cluck'; use NetCDF; my @types; $types[NetCDF::FLOAT] = 'float'; $types[NetCDF::DOUBLE] = 'double'; $types[NetCDF::LONG] = 'long'; $types[NetCDF::SHORT] = 'short'; $types[NetCDF::BYTE] = 'by te'; $types[NetCDF::CHAR] = 'char'; my $ncGLOBAL = NetCDF::GLOBAL; my $_MODE_RO = NetCDF::NOWRITE; sub new { my $proto = shi ft; my $class = ref($proto) || $proto; my $self = {}; if ( my $file = shift){ $self->{_FILE} = $file; } bless($self,$ class); return $self; } sub open { my ($self,$file) = @_; # may pass argument to open if ($file){ $self->{_FILE} = $file; } unless(-e $self->{_FILE}){ cluck "No such file: " . $self->{_FIL E}; return; # NetCDF::open will terminate!! if file can't be opened } # TODO: check to see if it looks like a netcdf file... # it shou ld be a binary file: # file tmp/test.nc = tmp/test.nc: data # and the first three chars are 'CDF' # this next line will crash on fa il, can't eval!! my $ncid = NetCDF::open($self->{_FILE},$_MODE_RO); # read-only mode should be default $self->{_NCID} = $ncid; # Store basic info my ($ndims,$nvars,$ngatts,$recdim); NetCDF::inquire($ncid ,$ndims,$nvars,$ngatts,$recdim); #warn "ndims=$ndims,nvars=$nvars,na tt=$ngatts,recdim=$recdim\n"; $self->{_NDIMS} = $ndims; # $se lf->{_NVARS} = $nvars; $self->{_NGATTS} = $ngatts; $self->{_RECD IM} = $recdim; my ($numrecs,$dim); # STORE OUR DIMENSION SIZ ES for (my $i=0;$i<$ndims;$i++){ my ($dname,$dsize); if (NetCD F::diminq($ncid,$i,\$dname,\$dsize) == 0){ # STORE FOR REFERENCE $self->{_DIM}{$dname} = $dsize; # STORE OUR RECORD DIMENSION: TELLS HOW MANY RECORDS WE HAVE if ($i == $recdim){ $numrecs = $dsize; $self->{_NUMRECS} = $numrecs; } } else { cluck "NetCDF::diminq returned non-zero"; } } # TODO: this could cause other problems... maybe we should yell return $self->{_NCID} unless ($self->{_NUMRECS}); # CREATE LIST OF ALL VARIABLES my @al l_vars = (); for my $i (0 ... ($self->{_NVARS} - 1)){ my ($var,$ty pe,$ndims,@dimids,$natts); NetCDF::varinq($ncid,$i,\$var,\$type,\$ndims ,\@dimids,\$natts); push @all_vars,$var; # also remember the type: used by get_rec_val $self->{_VAR_TYPES}{$var} = $type; } $self- >{_ALLVARS} = \@all_vars; # CREATE LIST OF RECORD VARIABLES # cu rrently only works if if (NetCDF::diminq($ncid,$self->{_RECDIM},\$dim,\ $numrecs) == 0){ # Get record zero and store var names my @recv arids = (); my @recsizes = (); my $nrecvars; NetCDF::recinq ($self->{_NCID}, $nrecvars, \@recvarids, \@recsizes) == 0 || cluck "Couldn't inquire about record variables\n"; my @varnames = (); if ($numrecs > 0){ my @record = (); NetCDF::recget($self-> {_NCID}, 0, \@record) == 0 || cluck "Couldn't get record zero\n"; # So we can immediately call get_next_rec() $self->{_CURRENT_RECORD_NU M} = -1; for my $k (0 ... $#record){ my $varid = $recvarids[ $k]; my ($var,$type,$ndims,@dimids,$natts); NetCDF::varinq($nci d,$varid,\$var,\$type,\$ndims,\@dimids,\$natts); push @varnames, $var ; } } $self->{_RECVARS} = \@varnames; } return $self->{_NCID}; # return 'true' (!=0) value if succeeded } sub clo se { my $self = shift; unless($self->{_NCID}){ cluck "don't know how we can close NetCDF with NCID=" . $self->{_NCID}; return; } my $rs = NetCDF::close($self->{_NCID}); return !$rs; } sub numrecs { # return number of records my $self = shift; r eturn $self->{_NUMRECS}; } =head2 get_rec_vars Get the list o f record variables for this file. =cut sub get_rec_vars { my $s elf = shift; my @vars = @{$self->{_RECVARS}}; return @vars; } =head2 get_all_vars Get the list of all variables for this fi le. =cut sub get_all_vars { my $self = shift; my @vars = @{$self->{_ALLVARS}}; return @vars; } =head2 get_gatt Get the values of a global attribute and return as a list. This is u seful when you just want to extract the value of a given attribute. NOTE : it returns a list, even if there is only one value. =cut sub g et_gatt { my ($self,$attr) = @_; my ($type,$len); # First get the type of the global attribute so we know how to handle the result my $ret = NetCDF::attinq($self->{_NCID}, NetCDF::GLOBAL, $attr, $type, $len ); my @values = (); # NetCDF::attget will core if this is an undef NetCDF::attget($self->{_NCID},NetCDF::GLOBAL,$attr,\@values); @value s = perlify_values($type,\@values); #print "gvalues(" . scalar @val ues . ")=(" . join(",",@values) . ")\n"; return @values; } =h ead2 get_all_gatts Get all global attributes and return as a hash. Sample: my %gatts = $nc->get_all_gatts(); foreach my $gatt ( sort keys %gatts){ my @vals = @{$gatts{$gatt}}; print "global att $gatt = " . join(",",@vals) . "\n"; } =cut sub get_all_gatts { my $self = shift; my %gatts = (); for my $i (0 ... ($ self->{_NGATTS} - 1)){ my ($name,$value,$type,$len); NetCDF::att name($self->{_NCID}, NetCDF::GLOBAL, $i, \$name); # First get the ty pe of the global attribute so we know how to present # the result m y $ret = NetCDF::attinq($self->{_NCID}, NetCDF::GLOBAL, $name, $type, $le n); my @values = (); # NetCDF::attget will core if this is an unde f NetCDF::attget($self->{_NCID},NetCDF::GLOBAL,$name,\@values); @va lues = perlify_values($type,\@values); $gatts{$name} = \@values; } return %gatts; } =head2 get_var_att Pass a name of a variable and the name of its attribute, and get back the value of the at tribute. This also checks to see if there is such an attribute, since th e charming NetCDF library will DIE if you try to look at a non-existent a ttribute. NOTE: it returns a list, even if there is only one value. =cut sub get_var_att { my ($self,$var,$attr) = @_; my @value s = (); # NetCDF::attget will core if this is an undef # GET VARIAB LE ID # which will kill you dead unless you have a valid variable if ( !( grep { $_ eq $var; } @{$self->{_ALLVARS}}) ){ return undef; } my $varid = NetCDF::varid($self->{_NCID}, $var); # Now get the number of attributes it has my ($var,$vtype,$ndims,@dimids,$natts); N etCDF::varinq($self->{_NCID},$varid,\$var,\$vtype,\$ndims,\@dimids,\$natts) ; #warn "varinq = $var,$vtype,$ndims,(" . join(",",@dimids) . "),$natt s" if $GoMOOS::DEBUG; # Find the correct attribute id... you'll see w hy in a minute my $attid = -1; for my $j (0 ... ($natts - 1)){ my $name; NetCDF::attname($self->{_NCID}, $varid, $j, \$name); if ( $name eq $attr){ $attid = $j; #warn "$j:name = $name\n"; last; } } # Return an empty list if there is no such attribute. if ($attid == -1){ return (); } #warn "attid = $attr\n" if $Go MOOS::DEBUG; # Now go after the attribute itself. my ($type,$len); # This one dies with an untrappable error if $attr is bogus... (eval ju st dies) # That's what all the above work is for! NetCDF::attinq($sel f->{_NCID}, $varid, $attr, $type, $len); # Now get the actual values NetCDF::attget($self->{_NCID},$varid,$attr,\@values); @values = perlif y_values($type,\@values); # Return single value as scalar for conveni ence if (@values == 1){ my $retval = $values[0]; return $r etval; } #warn "attr values(" . scalar @values . ")=" . join("," ,@values) if $GoMOOS::DEBUG; return @values; } =head2 get_rec Passed a record number: 0 to numrecs - 1 get the record and store it in _CURRENT_RECORD_NUM. See get_rec_val(). =cut sub get_re c { my ($self,$recnum, %opts) = @_; if( $self->{_NUMRECS} == 0){ cluck "NO RECORDS"; return 0; } if ($recnum < 0 || $r ecnum >= $self->{_NUMRECS}){ if (defined($self->{_NCID})){ unles s( $opts{INTERNAL}){ cluck "Invalid record number: $recnum"; } cluck "Invalid record number: $recnum\n"; $self->{_CURRENT_RECOR D_NUM} = -1; } else { #warn "No NCID: file not open!! (recnu m = $recnum)\n"; cluck "No NCID: file not open!! (recnum = $recnum )\n"; $self->{_CURRENT_RECORD_NUM} = -1; } return 0; } # GET RECORD my @current_record = (); my $status = NetCDF::re cget($self->{_NCID}, $recnum, \@current_record); if($status){ cluck "Could not get record $recnum\n"; $self->{_CURRENT_RECORD_NUM} = -1; return 0; } $self->{_CURRENT_RECORD_NUM} = $recnum; $sel f->{_CURRENT_RECORD} = \@current_record; return 1; } =head 2 get_rec_val Return the values from current record returned by the la st call to get_rec(). Pass a variable name and get a scalar if the variab le reference is a scalar or if it is an array with only one value, also r eturn a scalar. If it is a reference to an array just return the array re ference. Sample: my $sclr = $nc->get_rec_val('temperature_qc'); # time is an array but with only one value. my $sclr_time = $nc->get_r ec_val('time'); Note: @arry = $nc->get_rec_val('time') will also work, arry will have one value. my @arry = @{$nc->get_rec_val('current_u')} ; =cut sub get_rec_val { my ($self,$var) = @_; my @recv ars = @{$self->{_RECVARS}}; my @current_record = @{$self->{_CURRENT_ RECORD}}; my $type = $self->{_VAR_TYPES}{$var}; die "WHAT? NO RECOR D VARIABLES in $self->{_FILE}: (@recvars)" unless @recvars; my $retva l = undef; for my $k (0 ... $#recvars){ if ($recvars[$k] eq $var){ $retval = perlify_values($type,$current_record[$k]); last; } } #if (!defined($retval)){ cluck "did not find $var\n"; } re turn $retval; } sub get_next_rec { my $self = shift; if( $self->{_NUMRECS} == 0){ return 0; } $self->{_CURRENT_RECORD_NUM} ++; if ($self->{_CURRENT_RECORD_NUM} >= $self->{_NUMRECS}){ return 0; } return $self->get_rec($self->{_CURRENT_RECORD_NUM}, INTERNAL=> 1); } sub get_prev_rec { my $self = shift; if( $self->{_ NUMRECS} == 0){ return 0; } $self->{_CURRENT_RECORD_NUM}--; if ($self->{_CURRENT_RECORD_NUM} < 0){ return 0; } return $self->get _rec($self->{_CURRENT_RECORD_NUM}, INTERNAL=>1); } sub get_rewind { # reset get next record pointer my ($self,%opts) = @_; if (def ined($opts{REVERSE})){ $self->{_CURRENT_RECORD_NUM} = $self->{_NUMREC S}; } else { $self->{_CURRENT_RECORD_NUM} = 0; } retur n; } =head2 get_var ($var,[START=>$start,END=>$end,COUNT= >]) Return value(s) corresponding to a NetCDF variable. This accesses variables in an array format, as opposed to the record format. The parameters are: START, which tells which value (record) to start at (0 is t he first); TO, which gives the index of the last value to return; and COU NT, which tells a number of records to return. If none are specified, th e first record is returned. If both COUNT and TO are specified, COUNT ta kes precedence. If TO is -1, then values up to the last one are returned . If this is called in a scalar context, then the first value is retur ned. Otherwise, an array is returned. Caller beware. $start and $co unt are naturally an offset and a number; they default to 0 and 1, which gives the first value only, which is useful when there is only one, which is the only time you wouldn't want to speficy it anyway. If $count is - 1 then it retrieves the entire array. =cut sub get_var { my ( $self,$var,%range) = @_; my @values = (); # NetCDF::varget might cor e if this is an undef # careful... NetCDF::varget will CORE if start or count are undef... my ($start,$count); $start = 0 if not defined ($range{START}); if ($range{COUNT}){ $count = $range{COUNT}; } elsif ($range{TO}){ $count = $range{TO} - $range{START}; } $count = 1 if not $count; # GET VARIABLE ID # which will kill yo u dead unless you have a valid variable, so we check that # it's valid b efore we do anything with it if ( !( grep { $_ eq $var; } @{$self->{_ALL VARS}}) ){ return undef; } my $varid = NetCDF::varid($self->{_N CID}, $var); # GET NUMBER OF RECORDS IF POSSIBLE AND ADJUST COUNT IF NECESSARY if ($self->{_DIM}{$var}){ if ($self->{_DIM}{$var} < ($star t + $count) or ($range{TO} and not $range{COUNT})){ $count = $ self->{_DIM}{$var} - $start; } } # NOW CHECK START, COUNT, ETC . # NOW GET THE VALUES NetCDF::varget($self->{_NCID},$varid,\$start,\ $count,\@values); # And don't forget to perlify them @values = perl ify_values($self->{_VAR_TYPES}{$var},\@values); # Return single value as scalar for convenience if (not wantarray){ my $retval = $value s[0]; # see get_rec_val: here it could also be a string, so we vary the test... if (not $retval =~ /\w/ and $retval == 0){ return 0; } return $retval; } return @values; } =head2 perlify_v alues The purpose of this function is to prettify the values extracted from NetCDF files into something more like what a Perl programmer expect s - strings or numbers, not arrays of numbers representing char values or bytes that don't look like numbers. This expects two arguments: the value's NetCDF type, and a reference to the values. The values can be e ither an array OR scalar, depending on whether we've just gotten the cont ents of a NetCDF variable - which itself can be either - or the content o f a variable's attribute. This function returns either an array or a r eference to an array, e.g. @get_array = perlify_values($type,\@valu es); $get_ref = perlify_values($type,\@values); @use_array = @$ get_ref; =cut sub perlify_values { my ($type,$input) = @_; my @vals; my $result; # the values can be an array reference, a s calar reference, or a scalar. # just don't pass arrays. if (ref($i nput) eq 'ARRAY') { @vals = @$input } elsif (ref($input) eq 'SCALAR' ){ @vals = ($$input) } else { @vals = ($input) } if ($type == NetCDF::CHAR){ # weed out nulls, which can sneak in because NetC DF stores attribute strings # as arrays of chars, not as null-terminate d strings @vals = grep { $_ != 0 } @vals; # we also split it ba ck into an array of lines, which is the perl-ish way to # look at strin gs $result = [split "\n",(pack( 'c' x @vals, @vals))]; } # the numeric types are easy... the input is the same as the output elsif ($t ype == NetCDF::FLOAT){ $result = $input } elsif ($type == NetC DF::DOUBLE){ $result = $input } elsif ($type == NetCDF::SHORT){ $r esult = $input } elsif ($type == NetCDF::LONG){ $result = $inpu t } # except this one doesn't do what you'd expect in perl elsif ($ty pe == NetCDF::BYTE){ $result = [map { unpack('C',$_) } @vals]; } else { die "UNKNOWN NetCDF type='$type'"; } return (wa ntarray ? @$result : (@$result == 1 ? $result->[0] : $result)); } =head1 VERSION Revsion: $Revision: 1.29 $ =cut 1;