#
#   ArcInfoAscii.pm
#
#   ARC/INFO ASCII Grid データを利用する Perl のモジュール
# 
#   2011-09-16 by TAKENAKA, A.
#  
#  USAGE
#
#  my $ref = new ArcInfoAscii $file;  # オブジェクトの生成。$file はデータファイル名
#  exit unless($ref);                 # ファイルがない場合。
#
#  my ($x_min, $x_max)  = $ref->get_x_range();  # x 座標の範囲
#  my ($y_min, $y_max)  = $ref->get_y_range();  # y 座標の範囲
#  my $resol = $ref->get_resolution();          # 解像度（点の間隔）
#  my $no_data = $ref->get_NODATA_value();      # 未定義の場合の値
#
#  my $val = $ref->get_value(70, 130);             # (70, 130)の地点での値
#  my $val = $ref->get_value_at_nearest(70, 130);  # (70, 130)の最近地点での値
#
#   (70, 130) から半径 50以内の地点の値の配列
#  my @vals = $ref->get_values_in_buffer(70, 130, 50);  
#

package ArcInfoAscii;

use Exporter;

@ISA = qw(Exporter);

@EXPORT = qw(
new
get_value
get_value_at_nearest
get_values_in_buffer
get_resolution
get_x_range
get_y_range
get_NODATA_value
_load_ascii_grid_file
);


use strict;

my $default_NODATA = -9999;  # デフォルトの未定義値

#####################################
# ArcInfoAscii オブジェクトを生成するコンストラクタ
# パラメータとしてデータファイル名を渡す。
#
# USAGE
# $ref = new ArcInfoAscii $file_name;

sub new
{
    my $class =  shift @_;
    my $data_file = shift @_;

    unless ($data_file && -e $data_file) { 
        print "ArcInfoAscii::new()  File not given or does not exist.\n";
        return "";
    }

    my $self = {};
    &_load_ascii_grid_file($self, $data_file);
    bless $self, 'ArcInfoAscii';

    return $self;
}


#####################################
#
#  ArcInfo Ascii Grid ファイルから読み込んだデータセット(無名ハッシュへの
#  リファレンス）と、 x, y 座標を受けとり、その座標での値を計算して返す。
#
#  受けとった座標を 囲む４点の値から内挿（線形）で計算する。
#  データセットがカバーするエリア外であるか、 囲む４点中に NODATA が
#  あれば、NODATA を返す。
# 
#  USAGE: $refはArcInfoAscii オブジェクト
#  $val = $ref->get_value(x座標, y座標); # $refはArcInfoAscii オブジェクト

sub get_value
{
    my ($_ref, $x, $y) = @_;

    my $w =  $_ref->{cellsize};
    my $NODATA = $_ref->{NODATA_value};

    # 受けとった座標を番地（配列の添え字）に変換
    my $x_index = int(($x - $_ref->{xllcorner}) / $w);
    my $y_index = int(($y - $_ref->{yllcorner}) / $w);

    # エリア外であれば、NODATA値を返す。
    if ($x_index < 0 || $_ref->{ncols} <= $x_index + 1) {
        return  $NODATA;
    }

    if ($y_index < 0 || $_ref->{nrows} <= $y_index + 1) {
        return  $NODATA;
    }

    # エリア内であっても、指定した点を囲む４点のなかに NODATA値があれば、
    # NODATA値を返す。
    if (($_ref->{val}->[$y_index]->[$x_index] == $NODATA) ||
        ($_ref->{val}->[$y_index]->[$x_index + 1] == $NODATA) ||
        ($_ref->{val}->[$y_index + 1]->[$x_index] == $NODATA) ||
        ($_ref->{val}->[$y_index + 1]->[$x_index + 1] == $NODATA) ) {
        return $NODATA;
    }

    # 指定した点を囲む４点の座標を求める。
    my $x_0 = $_ref->{xllcorner} + $x_index * $w;
    my $x_1 = $_ref->{xllcorner} + ($x_index + 1) * $w;
    my $y_0 = $_ref->{yllcorner} + $y_index * $w;
    my $y_1 = $_ref->{yllcorner} + ($y_index + 1) * $w;


    # 囲む４点のうち、原点に近い点からの、x, y それぞれの距離を求める。
    my $x_d = $x - $x_0;
    my $y_d = $y - $y_0;

    # x が原点に近い辺上での、y 座標による内挿
    my $x0_y0_val = $_ref->{val}->[$y_index]->[$x_index];
    my $x0_y1_val = $_ref->{val}->[$y_index + 1]->[$x_index];
    my $x0_val = $x0_y0_val + ($x0_y1_val - $x0_y0_val) * $y_d / $w;

    # x が原点から遠い辺上での、y 座標による内挿
    my $x1_y0_val = $_ref->{val}->[$y_index]->[$x_index + 1];
    my $x1_y1_val = $_ref->{val}->[$y_index + 1]->[$x_index + 1];
    my $x1_val = $x1_y0_val + ($x1_y1_val - $x1_y0_val) * $y_d / $w;

    # 上の２つの内挿値から、さらに x 座標による内挿
    my $val = $x0_val + ($x1_val - $x0_val) * $x_d / $w;
 
    return $val;
}


#####################################
#
#  ArcInfo Ascii Grid ファイルから読み込んだデータセット(無名ハッシュへの
#  リファレンス）と、 x, y 座標を受けとり、その座標にもっとも近い点の値を返す。
#
#  データセットがカバーするエリア外であるか、 もっとも近い点の値が NODATA で
#  あれば、NODATA を返す。
# 
#
#  USAGE: $refはArcInfoAscii オブジェクト
#  $val = $ref->get_value_at_nearest(x座標, y座標); 

sub get_value_at_nearest
{
    my ($_ref, $x, $y) = @_;

    my $w =  $_ref->{cellsize};
    my $NODATA = $_ref->{NODATA_value};

    # 受けとった座標を番地（配列の添え字）に変換
    my $x_index = int(($x - $_ref->{xllcorner}) / $w);
    my $y_index = int(($y - $_ref->{yllcorner}) / $w);

    # エリア外であれば、NODATA値を返す。
    if ($x_index < 0 || $_ref->{ncols} <= $x_index + 1) {
        return  $NODATA;
    }

    if ($y_index < 0 || $_ref->{nrows} <= $y_index + 1) {
        return  $NODATA;
    }

    # 囲む４点のうち一番原点よりの点の座標を求める。
    my $x_0 = $_ref->{xllcorner} + $x_index * $w;
    my $y_0 = $_ref->{yllcorner} + $y_index * $w;

    # この点からの距離が、グリッド間隔の 1/2 より大きいかどうかで最近点を決める。
    # 中間の場合は原点に近いほうを優先。
    
    my $x_nearest_index = ($x - $x_0 <= $w / 2) ? $x_index : $x_index + 1;
    my $y_nearest_index = ($y - $y_0 <= $w / 2) ? $y_index : $y_index + 1;

    return  ($_ref->{val}->[$y_nearest_index]->[$x_nearest_index]);
}


#####################################
#
#  ArcInfo Ascii Grid ファイルから読み込んだデータセット(無名ハッシュへの
#  リファレンス）と、 x, y 座標、およびこれを中心とする円形バッファの半径を
#  受けとり、バッファに含まれる点の値の配列を返す。
#  ※円周上の点は含める
#
#  データセットのエリア外の点、および NODATA の点の値は含まない。
# 
#  USAGE: $refはArcInfoAscii オブジェクト
#  @values = $ref->get_values_in_buffer(x座標, y座標, 半径);

sub get_values_in_buffer
{
    my ($_ref, $x, $y, $radius) = @_;

    my $w =  $_ref->{cellsize};
    my $NODATA = $_ref->{NODATA_value};

    my @values = ();  #  値をためる配列

    my ($y_range_min, $y_range_max) = &get_y_range($_ref);
    my ($x_range_min, $x_range_max) = &get_x_range($_ref);

    # バッファの y 座標の最大・最小値
    my $y_min =  $y - $radius;
    my $y_max =  $y + $radius ;

    # バッファが完全に範囲外なら、すぐに終了
    if (($y_max < $y_range_min) || ($y_range_max < $y_min)) {
        return @values;
    }

    # 定義域からはみ出した部分の切り捨て
    $y_min = ($y_min < $y_range_min) ? $y_range_min : $y_min;
    $y_max = ($y_range_max < $y_max) ? $y_range_max : $y_max;

    # $y_min から $y_max に含まれる点の添え字
    my $y_index_min = int(($y_min - $_ref->{yllcorner}) / $w);
    my $y_index_max = int(($y_max - $_ref->{yllcorner}) / $w);

    # ちょうど円周上に点がある場合を除き、添え字を１増やしてバッファ内へ。
    my $y_std = ($y_min- $_ref->{yllcorner}) / $w;
    $y_index_min += 1 unless (int($y_std)  == $y_std); 

    # y 軸方向の添え字が小さいほうから大きいほうへ
    for (my $i_y = $y_index_min; $i_y <= $y_index_max; ++$i_y) {

        my $yy = $i_y * $w + $_ref->{yllcorner};        #  $i_y に対応する y 座標
        my $rt = sqrt($radius ** 2 - ($yy - $y) ** 2);  # x 軸と並行で $yy を通る弦の長さの半分

        my $x_min = $x - $rt;  # $yy を通る弦の両端の x 座標
        my $x_max = $x + $rt;

        if (($x_max < $x_range_min) || ($x_range_max < $x_min)) {
            next; # この yy に対しては、x 座標はすべて定義域外
        }

        # 定義域からはみ出した部分の切り捨て
        $x_min = ($x_min < $x_range_min) ? $x_range_min : $x_min;
        $x_max = ($x_range_max < $x_max) ? $x_range_max : $x_max;

        # $x_min から $x_max に含まれる点の添え字
        my $x_index_min = int(($x_min - $_ref->{xllcorner}) / $w);
        my $x_index_max = int(($x_max - $_ref->{xllcorner}) / $w);

        # ちょうど円周上に点がある場合を除き、添え字を１増やしてバッファ内へ。
        my $x_std = ($x_min- $_ref->{xllcorner}) / $w;
        $x_index_min += 1 unless (int($x_std) == $x_std); 

        # 各点の値を配列に加える。
        for (my $i_x = $x_index_min; $i_x <= $x_index_max; ++$i_x) {
            my $val = $_ref->{val}->[$i_y]->[$i_x];
            next if ($val == $NODATA);    # 値が NODATA ならスキップ
            push @values, $val;
        }
    }

    return @values;
}


#################################
#  ArcInfo Ascii Grid ファイルから読み込んだデータセットを受けとり、
#  解像度（点の間隔）を返す。
#
#  USAGE: $refはArcInfoAscii オブジェクト
#  $val = $ref->get_resolution(x座標, y座標); 

sub get_resolution
{
    my $_ref = shift @_;
    return  $_ref->{cellsize};
}


#################################
#  ArcInfo Ascii Grid ファイルから読み込んだデータセットを受けとり、
#  x 軸方向の座標の最小・最大を返す。
#
# USAGE: $refはArcInfoAscii オブジェクト
# ($x_min, $x_max) = $ref->get_x_range(); 

sub get_x_range
{
    my $_ref = shift @_;

    my $x_min = $_ref->{xllcorner};
    my $x_max = $x_min + ($_ref->{ncols} - 1) * $_ref->{cellsize};
    
    return ($x_min, $x_max);
}


#################################
#  ArcInfo Ascii Grid ファイルから読み込んだデータセットを受けとり、
# y 軸方向の座標の最小・最大を返す。
#
#  USAGE: $refはArcInfoAscii オブジェクト
#  ($y_min, $y_max) = $ref->get_y_range(); 

sub get_y_range
{
    my $_ref = shift @_;

    my $y_min = $_ref->{yllcorner};
    my $y_max = $y_min + ($_ref->{nrows} - 1) * $_ref->{cellsize};

    return ($y_min, $y_max);
}


#################################
#
#  NODATA の値を返す。
#    ファイルを読んだあとは、そのファイルで定義されている値。
#    そうでなければ、デフォルト値（-9999)
#
#  USAGE: $refはArcInfoAscii オブジェクト
#  $no_data = $ref->get_NODATA_value(); 

sub get_NODATA_value
{
    my $_ref = shift @_;
    return  $_ref->{NODATA_value};
}

#####################################
# Arc/Info ASCII Grid ファイルの読み込み。
# データを読み込んだ無名ハッシュへのリファレンスを返す。
# new から呼び出される。原則として直接は使わない。

sub _load_ascii_grid_file
{
    my ($_ref)    = shift @_;
    my $data_file = shift @_;

    unless (-e $data_file) {
        die ("$data_file not found");
    }

    open (my $fh, "<", $data_file);

    my $record;
    my ($item, $val);
    my $array_ref = [];  # 無名配列のレファレンス

    for (my $i = 0; $i < 6; ++$i) {
        $record = <$fh>;
        chomp $record;
        ($item, $val) = split(/\s+/, $record);
        $_ref->{$item} = $val;
    }

    unless ($_ref->{nrows} && $_ref->{ncols} & $_ref->{cellsize}) {
        die ("$data_file is not ArcINFO ASCII grid file.");
    }

 #  NODATA_value の行がない場合、最後に読んだ６番めの行はデータ行
    unless ($_ref->{NODATA_value}) {
        my @data = split(/\s+/, $record);
        push @$array_ref, \@data;

        $_ref->{NODATA_value} = $default_NODATA; # デフォルトの値を設定
    }

    while ($record = <$fh>) {  # データ読み込み
        chomp $record;
        my @data = split(/\s+/, $record);
        push @$array_ref, \@data;
    }

    $_ref->{val} = $array_ref;  # 
    close($fh);

    return $_ref;
}


1;