如何通过 bash 脚本中的坐标 args 访问网格文件？答案

【问题标题】：How to access grid files by coordinate args in a bash script?如何通过 bash 脚本中的坐标 args 访问网格文件？
【发布时间】：2016-01-25 04:27:19
【问题描述】：

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

我的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 会更简单，但事实并非如此。

【问题讨论】：

不要使用 bash，它只会让你感到沮丧和烦扰计算机。既然您已经了解 awk 的基本原理，请尝试 Perl；这就是它的设计目的。
我所说的“惹恼计算机”是指“在您的“巨大文件”上“将花费真正很长时间”。
+1 用于示范问题。一切都在那里，数据、预期输出和喘气，代码！一件事..，所以数据格式在一个文件中一遍又一遍地重复？你无法控制这个文件是如何呈现给你的？如果我真的受限于命令行技术，我认为这可以在 awk 中完成，并且会考虑将数据展平，以便每组都在一条线上。假设这不是一次性的一次性任务，我会将数据放入数据库中。一个基本说明：任何$(... anything ..) 的使用都会创建另一个进程并会显着减慢速度。祝你好运。
+1 关于这个问题。你能解释一下3356385.137 5800799.818如何返回31.11266和3356387.137 5800803.818的值返回29.58562
@Jaypal 将数据矩阵视为坐标系。左上角的第一个数据具有 (x, y, z) 值 (3356385.137, 5800799.818, 31.11266) 和模拟： (3356387.137, 5800803.818, 29.58562) 是中间最后一行的数据三元组。知道了？ xllcorner 和 yllcorner 是起始坐标，在左上角，size 是从 x0 到 x1 到 x2 和 y0 到 y1 到 y2 的步长大小，依此类推，在本例中为 +1.0。

标签： bash matrix geolocation grid

【解决方案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

【讨论】：

这是一个非常好用的脚本，我刚刚测试过了。但是内存似乎存在很大的问题。在一个 2.6GB 的网格文件上，我需要大约 26GB 主内存。怎么样？
另外，我的 shell 脚本需要 7 秒才能产生一个数据值，而您的 perl 脚本需要大约 2 分钟。从 18 个或更多请求值开始，您的脚本将比我的更省时，因为您的脚本时间与输入参数的数量无关。如果在可用容量较低的机器上使用，巨大的内存消耗将是一个问题。有什么想法吗？
有了这么多数据，使用真实的数据库会让您受益匪浅。
顺便说一句，您可以通过将文件保存为行数组而不是数组数组来节省一些内存，并且仅在需要访问x 坐标。

【解决方案2】：

好的，了解我自己。在上面的代码中只添加 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

【讨论】：