awk运用三维数组进行插值获得任意经纬度处的水层沉积层地壳厚度

awk三维数组与插值

目的：给定经纬度，获得该点地下的冰层水层沉积层和地壳的厚度
实现：awk一行命令

下载Crust1.0模型

该数据集的详细介绍见官网.
解压后有几个文件:crust1.vp,crust1.vs,crust1.rho,crust1.bnds，分别代表P波速度，S波速度，密度和深度.

数据格式

每个文件共有64800行，9列，每行经纬度如下表所列。crust1.bnds的每行代表某个经纬度所对应的水层，冰层，上中下沉积层和上中下地壳的深度。

行号	latitude	lat_index	longitude	lon_index
1	89.5	180	-179.5	0
2	89.5	180	-178.5	1
---	---	---	---	---
360	89.5	180	179.5	359
361	88.5	179	-179.5	0
362	89.5	179	-178.5	1
---	---	---	---	---
64800	-89.5	0	179.5	359

awk三维数组存储

在BEGIN{}里将三维数组初始化，在{}里将数据存储到三维数组中，在END{}里线性插值

lat=-89.5lon=-179.5cat /opt/crust10/crust1.bnds | awk 'BEGIN{  for (lat_index=0; lat_index<180; lat_index++) {    for (lon_index=0; lon_index<360; lon_index++) {      for (i=1; i<=9; i++) {        value[lat_index,lon_index,i]=0.0      }    }  } }{  lat_index=180-int(NR/360);   lon_index=NR%360;   for (i=1; i<=9; i++) {    value[lat_index,lon_index,i]=$i  }}END{  lat_index1=int('"$lat"'+89.5)  lat_index2=lat_index+1  lon_index1=int('"$lon"'+179.5)  lon_index2=lon_index+1  for (i=1; i<9; i++) {    lat1='"$lat"'+89.5-lat_index1    lat2=1.0-lat1    lon1='"$lon"'+179.5-lon_index1    lon2=1.0-lon1    print i,value[lat_index1,lon_index1,i]*lon2*lat2+value[lat_index1,lon_index2,i]*lon1*lat2+value[lat_index2,lon_index1,i]*lon2*lat1+value[lat_index2,lon_index2,i]*lon1*lat1  }}' | cat