grd2vtk/lonlatzvevn2vtk.awk at master · thwbecker/grd2vtk · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#
#
# convert lon lat depth(m) v_east v_north to VTK vectors
#
# Thorsten Becker, UT Austin
# thorstinski@gmail.com
#
# $Id: gmtpoly2vtk.awk,v 1.2 2008/10/16 15:51:22 becker Exp becker $
#

BEGIN{
  # counters
  nsc=0;
  nsi=0;
  np=0;
  R = 1;
  # constants
  f=0.017453292519943296;
  #
  use_attribute=0

  vert_ex = 25;			# vertical exageration
}
{
  if(substr($1,1)!="#"){


      np++;

      lambda=$2*f;
      theta=(90.0-$2)*f;
      phi=$1*f;

      r = (1+$3*vert_ex/6371e3);


# polar velocity components
      vr=0;
      vtheta=-$5;
      vphi=$4;
#
# base vecs
      ct=cos(theta);
      cp=cos(phi);
      st=sin(theta);
      sp=sin(phi);
#
      polar_base_r[1]= st * cp;
      polar_base_r[2]= st * sp;
      polar_base_r[3]= ct;
#
      polar_base_theta[1]= ct * cp;
      polar_base_theta[2]= ct * sp;
      polar_base_theta[3]= -st;
#
      polar_base_phi[1]= -sp;
      polar_base_phi[2]=  cp;
      polar_base_phi[3]= 0.0;
# convert vector

      for(i=1;i<=3;i++){
	  cart_vec[np*3+i]  = polar_base_r[i]    * vr ;
	  cart_vec[np*3+i] += polar_base_theta[i]* vtheta;
	  cart_vec[np*3+i] += polar_base_phi[i]  * vphi;
      }


      tmp=cos(lambda) * r;
      x[np]=tmp * cos(phi);
      y[np]=tmp * sin(phi);
      z[np]=sin(lambda)*r;

    }

}
END{
  print("# vtk DataFile Version 4.0");
  print("converted from GMT -M file");
  print("ASCII");
  print("DATASET POLYDATA");
  print("POINTS",np,"float")
  for(i=1;i<=np;i++)
    printf("%.6e %.6e %.6e\n",x[i],y[i],z[i]);
  print("");

  print("POINT_DATA ",np)
  print("VECTORS velocity float");
  for(i=1;i<=np;i++){
      printf("%g %g %g\n",cart_vec[i*3+1],cart_vec[i*3+2],cart_vec[i*3+3]);
  }


}