2

我有 2 个数据文件：file01和file02. 在这两个数据集中，字段是： (i) 标识符；(ii) 数字参考；(iii) 经度；(iv) 纬度。对于 中的每一行file01，我想file02使用相同的数字引用来搜索数据，然后找到file02最接近 中的标识符的标识符file01。

如果我file01使用以下代码将值手动传递给 awk 程序，我可以得到这个：

awk 'function acos(x) { return atan2(sqrt(1-x*x), x) }
BEGIN {pi=3.14159;
       ndist=999999999.1;
       date=1001;
       lo1=-1.20; lg1=lo1*(pi/180);
       la1=30.31; lt1=la1*(pi/180)
           }
{if($2==date) {ws=$1;
               lg2=$3*(pi/180);
               lt2=$4*(pi/180);
               dist= 6378.7 * acos( sin(lt1)*sin(lt2) + cos(lt1)*cos(lt2)*cos(lg2-lg1) );
               if(dist < ndist) {ndist=dist; ws0=ws}}}
END {print(ws0,ndist)}' file02

如您所见，date语句中的lo1和la1是BEGIN第 1 行中的值file01（请参阅下面的数据文件）。我的问题是我是否可以一次做到这一点，所以每次我在 中读取一行时file01，我都会获得最近的标识符和距离并附加到 中的行数据file01。我不知道某些 shell 命令是否可以以更简单的方式实现这一点，也许使用管道。

这两个数据文件的示例和所需的输出是：

=== 文件01 ===

A 1001 -1.2 30.31
A 1002 -1.2 30.31
B 1002 -1.8 30.82
B 1003 -1.8 30.82
C 1001 -2.1 28.55

=== 文件02 ===

ws1 1000 -1.3 29.01
ws1 1001 -1.3 29.01
ws1 1002 -1.3 29.01
ws1 1003 -1.3 29.01
ws1 1004 -1.3 29.01
ws1 1005 -1.3 29.01
ws2 1000 -1.5 30.12
ws2 1002 -1.5 30.12
ws2 1003 -1.5 30.12
ws2 1004 -1.5 30.12
ws2 1005 -1.5 30.12
ws3 1000 -1.7 29.55
ws3 1001 -1.7 29.55
ws3 1002 -1.7 29.55
ws3 1003 -1.7 29.55
ws3 1004 -1.7 29.55
ws3 1005 -1.7 29.55
ws4 1000 -1.9 30.33
ws4 1001 -1.9 30.33
ws4 1002 -1.9 30.33
ws4 1003 -1.9 30.33
ws4 1004 -1.9 30.33
ws4 1005 -1.9 30.33

=== 输出文件 ===

A 1001 -1.2 30.31 ws4 67.308
A 1002 -1.2 30.31 ws2 35.783
B 1002 -1.8 30.82 ws4 55.387
B 1003 -1.8 30.82 ws4 55.387
C 1001 -2.1 28.55 ws1 85.369

编辑#1：考虑到@Eran 的建议，我编写了以下代码：

join -j 2 < (sort -k 2,2 file01) < (sort -k 2,2 file02) |
awk 'function acos(x) { return atan2(sqrt(1-x*x), x) }
     BEGIN {pi=3.14159}

     {if (last != $1 $2)
         {print NR, id,r,lon,lat,ws0,ndist;
          last = $1 $2;
          ndist=999999999.1

         } else {

          lg1=$3*(pi/180);
          lt1=$4*(pi/180);
          lg2=$6*(pi/180);
          lt2=$7*(pi/180);
          dist= 6378.7 * acos( sin(lt1)*sin(lt2) + cos(lt1)*cos(lt2)*cos(lg2-lg1) );
          if(dist< ndist) {ndist=dist; ws0=$5}
          id=$2;r=$1;lon=$3;lat=$4

          }
     }'

该脚本的输出是：

1      
4  A 1001 -1.2 30.31 ws4 67.3078
7  C 1001 -2.0 28.55 ws3 115.094
11 A 1002 -1.2 30.31 ws2 35.7827
15 B 1002 -1.8 30.82 ws4 55.387

编辑#2：使用@Dennis 的建议（经过一些修改），我得到了所需的输出。awk 脚本如下：


awk 'function acos(x) { return atan2(sqrt(1-x*x), x) }
     BEGIN {pi=3.14159}
     NR==FNR {c++; a1[c]=$1;a2[c]=$2;a3[c]=$3;a4[c]=$4; next}
             {d++; b1[d]=$1;b2[d]=$2;b3[d]=$3;b4[d]=$4}

     END {
     for(k=1;k<=c;k++) {
         lg1=a3[k]*(pi/180);
         lt1=a4[k]*(pi/180);
         ndist=999999999.1;
         for(l=1;l<=d;l++) {
             if(b2[l]==a2[k]) {kk=b2[l];
                lg2=b3[l]*(pi/180);
                lt2=b4[l]*(pi/180);
                dist= 6378.7 * acos( sin(lt1)*sin(lt2) + cos(lt1)*cos(lt2)*cos(lg2-lg1) );
                if(dist<ndist) {ndist=dist; ws0=b1[l]}
             }
         }
         print a1[k],a2[k],a3[k],a4[k],ws0,ndist
     }
    }' file01 file02
4

4 回答 4

2

将您的值从 file01 读入一个或多个数组。您可以getline在BEGIN块中使用，或者规范的方式是使用FNR == NR循环作为主要块之一。

FNR == NR {array[$1] = $1; ...; next } # read file01 into some arrays
{ for (item in array) { ... }     # process each item in the array(s) against each line in file02

您的脚本将被调用为awk '...' file01 file02

array1[c] = $1; array2[c] = $2; c++您可以使用计数器索引它们并使用计数器迭代而不是使用in: ，而不是按字段值索引数组for (i=0; i<c; i++)。

当然，您应该选择有意义的数组名称。

于 2012-05-09T16:04:23.843 回答
1

有趣的挑战。由于您必须先读取 file02 并将信息存储在数据结构中，所以我首先倾向于 Perl。

#!/usr/bin/perl
use strict;
use warnings;

# see http://perldoc.perl.org/Math/Trig.html
use Math::Trig qw(great_circle_distance deg2rad);
sub NESW {deg2rad($_[0]), deg2rad(90-$_[1])}

# read file02
my %data;
my $file2 = 'file02';
open my $fid, '<', $file2 or die "can't read $file2: $!\n";
while (<$fid>) {
    my ($id, $ref, $long, $lat) = split;
    push @{$data{$ref}}, [$id, $long, $lat];
}
close $fid;

$, = " ";

# process file01
my $file1 = 'file01';
open $fid, '<', $file1 or die "can't read $file1: $!\n";
while (<$fid>) {
    my ($id, $ref, $long, $lat) = split;
    my @here = NESW($long, $lat);
    my $min = 99_999_999;
    my (@min_id, $dist);

    while (my ($key, $listref) = each %data) {
        next unless $key == $ref;

        foreach my $trioref (@$listref) {
            my ($other_id, $other_long, $other_lat) = @$trioref;
            my @there = NESW($other_long, $other_lat);
            $dist = great_circle_distance(@here, @there, 6378.7);
            if ($dist < $min) {
                $min = $dist;
                @min_id = @$trioref;
            }
        }
    }

    printf "%s %d %s %s %s %6.3f\n", $id, $ref, $long, $lat, $min_id, $min;
}
close $fid;

这输出

A 1001 -1.2 30.31 ws4 67.308
A 1002 -1.2 30.31 ws2 35.783
B 1002 -1.8 30.82 ws4 55.387
B 1003 -1.8 30.82 ws4 55.387
C 1001 -2.1 28.55 ws1 93.361

我注意到“C”距离与您建议的距离不同。

于 2012-05-09T16:47:06.057 回答
1

要立即执行，请运行

join -j 2 <(sort -k 2,2 file01) <(sort -k 2,2 file02)

并将其通过管道传输到 awk，在每次更改引用时执行计算：

gawk '{if (last != $1 $2) {calc_nearest_on_array; last=$1 $2; add_point_to_array} else {add_point_to_array}}'
于 2012-05-10T05:13:14.577 回答
0

发射器：

@(do
   (defvar pi 3.1415926535)
   (defvar earth-radius 6378.7)
   (defun rad (deg) (/ (* deg pi) 180))
   (defun sphere-distance (lat0 lon0 lat1 lon1)
     (let ((lat0 (rad lat0)) (lat1 (rad lat1))
           (lon0 (rad lon0)) (lon1 (rad lon1)))
       (* earth-radius (acos (+ (* (sin lat0) (sin lat1))
                                (* (cos lat0) (cos lat1) (cos (- lon1 lon0)))))))))
@(next "file01")
@(collect)
@id @ref @lon0 @lat0
@  (filter :tonumber lon0 lat0)
@  (next "file02")
@  (bind (min-dist ws) (1e99 nil))
@  (collect)
@ws1 @ref @lon1 @lat1
@    (filter :tonumber lon1 lat1)
@    (do (let ((d (sphere-distance lat0 lon0 lat1 lon1)))
           (cond ((< d min-dist)
                  (set min-dist d)
                  (set ws ws1)))))
@  (end)
@  (do (format t "~a ~a ~0,2f ~0,2f ~a ~0,3f\n" id ref lon0 lat0 ws min-dist))
@(end)

跑：

$ txr dist.txr
A 1001 -1.20 30.31 ws4 67.308
A 1002 -1.20 30.31 ws2 35.783
B 1002 -1.80 30.82 ws4 55.387
B 1003 -1.80 30.82 ws4 55.387
C 1001 -2.10 28.55 ws1 93.361
于 2012-05-11T07:09:27.663 回答