1

我有一个巨大的ESRI Grid文件，其中包含空格分隔的数据。elevation.asc为了简化这个问题，我将使用一个具有 5x5 值的示例文件。

my 的标题elevation.asc包含一些关于数据的附加信息，包括第一个值的起始坐标（纬度、经度）。完整的文件如下所示：

ncols         5
nrows         5
xllcorner     3356385.137
yllcorner     5800799.818
cellsize      1.0
NODATA_value  -9999
31.11266 31.03987 31.15038 30.98865 30.96297
29.65054 29.65345 29.65598 29.60781 29.61685
29.70712 29.66978 29.73194 29.83858 29.87868
29.54893 29.60815 29.62812 29.66953 29.70786
29.55878 29.55927 29.58562 29.66112 29.79232

现在我的问题是，如何使用bash脚本通过给定坐标访问此文件中的高程数据？

我想调用我的脚本并产生这样的第一个数据值usage: myscript.sh {file} {x} {y}：

 $ ./myscript.sh elevation.asc 3356385.137 5800799.818
31.11266

或者：

 $ ./myscript.sh elevation.asc 3356387.137 5800803.818
29.58562

现在，到目前为止我尝试了什么？我正在玩while循环遍历 x 和 y 坐标并awk解析标题并bc进行一些浮点计算。但我现在有点迷失了如何继续。这就是我得到的：

#!/bin/bash
E_BADARGS=65
E_NOINPUT=66
N_ARGS=3

# Checks for proper number of command line args.
if [ $# -ne $N_ARGS ] ; then
  echo "Usage: `basename $0` {input.file} {x} {y}"
  exit $E_BADARGS
fi

# Checks for proper input file.
INPUT=$1
[ ! -f $INPUT ] && { echo "$INPUT file not found"; exit $E_NOINPUT; }

# Parses file header info.
cols=$(awk '$1 == "ncols" { print $2 }' $INPUT)
rows=$(awk '$1 == "nrows" { print $2 }' $INPUT)
x11=$(awk '$1 == "xllcorner" { print $2 }' $INPUT)
y11=$(awk '$1 == "yllcorner" { print $2 }' $INPUT)
size=$(awk '$1 == "cellsize" { print $2 }' $INPUT)
nodata=$(awk '$1 == "NODATA_value" { print $2 }' $INPUT)

# Calculates maximum coordinates.
xpp=$(bc <<< "$x11+$cols-1")
ypp=$(bc <<< "$y11+$rows-1")

# Gets requested coordinates from args.
X=$2
Y=$3

### What now?

但是，是的，现在呢？我可以使用循环遍历整个栅格while以找出坐标的位置，但后来我注意到我所做的只是找到输入坐标而不是存储的数据值。

# Iterates through the whole raster.
while [ $(echo "$y11 < $ypp" | bc) == 1 ] ; do
  while [ $(echo "$x11 < $xpp" | bc) == 1 ] ; do
    if [ $X == $x11 ] && [ $Y == $y11 ] ; then
      ### What now?
    fi
    x11=$(bc <<< "$x11+$size")
  done
  y11=$(bc <<< "$y11+$size")
done

我不知道该怎么做。如何使用 bash 脚本访问高程数据？

更新：澄清

此问题顶部给定的 5x5 矩阵是表示数字高程模型的数据图。每个值表示海拔米。

输入和输出：我用三个参数调用我的脚本：文件名、纬度坐标 (x) 和经度坐标 (y)。像那样：

 $ ./myscript.sh elevation.asc 3356385.137 5800799.818

现在在标题中定义了数据集从左上角（第一个坐标）的坐标xllcorner=3356385.137和开始yllcorner=5800799.818。因此，使用这两个坐标调用脚本应该会在 5x5 矩阵的左上角产生第一个海拔数据 (z)，即z=31.11266. size是 x 和 y 方向上 2 个数据字段之间的步长，以米为单位。所以调用这个：

 $ ./myscript.sh elevation.asc 3356387.137 5800803.818

... 意味着在 x 方向上走两步，在 y 方向上走 4 步，这会产生z=29.58562. 只需在矩阵中计算出来。

如果xllcorner并且yllcorner会更简单，0但事实并非如此。

4

2 回答 2

1

这是一个 perl 解决方案。

#!/usr/bin/perl
use strict;
use warnings;
use v5.10;
use autodie;
use List::MoreUtils qw(any);

my $data_file = shift;
my %metadata;
my @data;

open my $fh, '<', $data_file;
while (<$fh>) {
    chomp;
    my @F = split;
    if (any {$F[0] eq $_} qw(ncols nrows xllcorner yllcorner cellsize NODATA_value)) {
        $metadata{$F[0]} = $F[1];
    }
    else {
        push @data, \@F;
    }
}
close $fh;

while (@ARGV) {
    my $x = shift;
    my $y = shift;
    my $x_delta = int(($x - $metadata{xllcorner}) / $metadata{cellsize});
    my $y_delta = int(($y - $metadata{yllcorner}) / $metadata{cellsize});
    if ($x_delta < 0 or $y_delta < 0 or not defined $data[$y_delta][$x_delta]) {
        say $metadata{NODATA_value};
    }
    else {
        say $data[$y_delta][$x_delta];
    }
}

由于读取所有数据可能会很昂贵，因此您应该能够一次传递几对坐标：这就是使用 @ARGV 的 while 循环。这给了你：

$ perl esri.pl elevation.asc 0 0 3356385.137 5800799.818 3356387.137 5800803.818
-9999
31.11266
29.58562
于 2013-06-16T18:24:30.723 回答
1

好吧，想通了我的自我。在上面的代码中只添加 3 行代码非常简单：

# Calculates the line number and column.
LINE=$(bc <<< "scale=0;($Y-$y11+7)/1")
COLN=$(bc <<< "scale=0;($X-$x11+1)/1")

# Prints the Z Coordinates at X=COLN and Y=LINE
awk -v line=$LINE 'NR == line { print $0 }' $INPUT  | cut -f $COLN -d " "

只需通过减去输入$X和$Y起始坐标x11等来计算高程数据值的位置，y11就像上面写的那样。

整个工作代码是：

#!/bin/bash
INPUT=$1

# Parses starting coordinates.
x11=$(awk '$1 == "xllcorner" { print $2 }' $INPUT)
y11=$(awk '$1 == "yllcorner" { print $2 }' $INPUT)

# Gets requested coordinates from args.
X=$2
Y=$3

# Calculates the line number and column.
LINE=$(bc <<< "scale=0;($Y-$y11+7)/1")
COLN=$(bc <<< "scale=0;($X-$x11+1)/1")

# Prints the Z Coordinates at X=COLN and Y=LINE
awk -v line=$LINE 'NR == line { print $0 }' $INPUT  | cut -f $COLN -d " "

参数是：{file} {x} {y} 示例用法：./myscript.sh elevation.asc 3356387.137 5800803.818

于 2013-06-16T16:03:15.363 回答