#! /usr/bin/env perl
#
#    This example reads some surface pressure and temperatures. It is
#    intended to illustrate the use of the netCDF perl API. The companion
#    program sfc_pres_temp_wr.pl shows how to write the netCDF data file
#    needed by this program.
#
# Heiko Klein <Heiko.Klein at met.no> 2007-02-06
#

use strict;
use warnings;

use PDL::Lite;
use PDL::NetCDF;
use Fcntl;

my $NLAT = 6;
my $NLON = 12;
my $SAMPLE_PRESSURE = 900.0;
my $SAMPLE_TEMP = 9.0;
my $START_LAT = 25.0;
my $START_LON = -125.0;

my $latitude = (PDL::Basic::sequence($NLAT) * 5 + $START_LAT)->float;
my $longitude = (PDL::Basic::sequence($NLON) * 5 + $START_LON)->float;

my $pressure = PDL::Core::zeroes($NLON, $NLAT)->float;
my $temp = PDL::Core::zeroes($NLON, $NLAT)->float;

for (my $i = 0; $i < $NLAT; $i++) {
    for (my $j = 0; $j < $NLON; $j++) {
	$pressure->mslice($j, $i) .= $SAMPLE_PRESSURE + $i + $j * 5;
	$temp->mslice($j, $i) .= $SAMPLE_TEMP + $i * .25 + $j * .5;
    }
}

eval {
    my $ncfile = new PDL::NetCDF("sfc_pres_temp.nc", {REVERSE_DIMS => 1,
						      MODE => O_RDONLY});

    # checking data
    if (not equalPdls($latitude, $ncfile->get('latitude'))) {
	die "different latitude values";
    }
    if (not equalPdls($longitude, $ncfile->get('longitude'))) {
	die "different longitude values";
    }
    if (not equalPdls($pressure, $ncfile->get('pressure'))) {
	die "different pressure values";
    }
    if (not equalPdls($temp, $ncfile->get('temperature'))) {
	die "different temperature values";
    }

    # checking attributes
    if ('degrees_north' ne $ncfile->getatt('units', 'latitude')) {
	die "wrong units for latitude";
    }
    if ('degrees_east' ne $ncfile->getatt('units', 'longitude')) {
	die "wrong units for longtude";
    }
    if ('hPa' ne $ncfile->getatt('units', 'pressure')) {
	die "wrong units for pressure";
    }
    if ('celcius' ne $ncfile->getatt('units', 'temperature')) {
	die "wrong units for temperature";
    }
    $ncfile->close;
};

if ($@) {
    print STDERR "$@\n";
} else {
    print "SUCCESS reading file\n";
}


sub equalPdls {
    my ($pdlA, $pdlB) = @_;
    return $pdlA->where($pdlA != $pdlB)->isempty;
}