項目要用到matlab中的Vq = interp2(X,Y,V,Xq,Yq)函數，即把一個已知經緯度和對應值的矩陣，插值變換到一個給定經緯度網格中，也就是對給定網格填值，需要用到插值，這里使用雙線性內插法。
*(這只是一個初步完成代碼，僅供參考)
以下是對應C代碼和測試程序：
1 //************************************ 2 // 函數名稱: inter_linear() 3 // 函數說明：計算兩點之間某一給定點的值，（x0,y0）->(x1,y1),已知x(x0<x<x1),求y 4 // 返 回 值: double 5 // 參 數: double x0,y0,x1,y1/*in*/ 兩個點的坐標（x0,y0）（x1,y1） 6 7 // 作 者：WSS 8 // 作成日期：2019/09/05 9 10 //************************************ 11 double inter_linear(double x0, double y0, double x1, double y1, double x) 12 { 13 double a0, a1, y; 14 a0 = (x - x1) / (x0 - x1); 15 a1 = (x - x0) / (x1 - x0); 16 y = a0*y0 + a1*y1; 17 return y; 18 }
1 //************************************ 2 // 函數名稱: interp2d() 3 // 函數說明：二維插值，同matlab的interp2()功能 4 // 返 回 值: double 5 // 參 數: x,y分別為長度為m和n的向量（一維數組），z為矩陣（對應的二維數組(m,n)） 6 // a,b分別為長度為asize和bsize的向量（一維數組），out_result為矩陣（對應的二維數組(asize,bsize)） 7 8 // 作 者：WSS 9 // 作成日期：2019/09/05 10 11 //************************************ 12 int interp2d(double *x, double *y, double *z, int m, int n, double *a, double *b, int asize,int bsize,double *out_result) 13 { 14 //a,b是待插值的矩陣行和列，遍歷a，取出行，遍歷b取出列，將插值計算后的值存儲在out_result中 15 double pointa = 0; 16 double pointb = 0; 17 double w1, w2,w; 18 const int nu_val = -9999; 19 int tempi = 0; 20 int tempj =0; 21 for (int i = 0; i < asize; ++i) 22 { 23 for (int j = 0; j < bsize; ++j) 24 { 25 // 找到網格點位置 26 pointa = a[i]; 27 pointb = b[j]; 28 //確定pointa pointb所在的位置 29 for (int p = 0; p < m-1; ++p) 30 { 31 if (pointa < x[0]||pointa>x[m-1])//范圍外的值不進行插值處理 32 { 33 tempi = nu_val; 34 break; 35 } 36 else if (pointa == x[m - 1]) 37 { 38 tempi = m - 1; 39 break; 40 } 41 else if(pointa >= x[p] && pointa < x[p+1])//x升序 42 { 43 tempi = p; 44 break; 45 } 46 } 47 for (int q = 0; q < n - 1; ++q) 48 { 49 if (pointb< y[0] || pointb>y[n- 1])//范圍外的值不進行插值處理 50 { 51 tempj = nu_val; 52 break; 53 } 54 else if (pointb == y[n - 1]) 55 { 56 tempj = n - 1; 57 break; 58 } 59 else if (pointb >= y[q] && pointb < y[q + 1])//y升序 60 { 61 tempj = q; 62 break; 63 } 64 } 65 if (tempj == nu_val || tempi == nu_val) 66 { 67 out_result[i*bsize + j] = nu_val; 68 } 69 else 70 { 71 //如果四個點，有一個點為-9999，就不進行插值了 72 if(z[tempi*n + tempj]<-900||z[tempi*n + tempj + 1]<-900||z[(tempi + 1)*n + tempj]<-900||z[(tempi + 1)*n + tempj + 1]<-900) 73 { 74 out_result[i*bsize + j] = nu_val; 75 } 76 else 77 { 78 //cout << tempi << " " << tempj << endl; 79 /**************x方向進行插值*************/ 80 if (x[tempi] == pointa) 81 { 82 //取網格節點值 83 w1 = z[tempi*n + tempj]; 84 w2 = z[tempi*n + tempj + 1]; 85 //y方向進行插值 86 if (y[tempj] == pointb) 87 { 88 w = w1; 89 } 90 else 91 { 92 w = inter_linear(y[tempj], w1, y[tempj + 1], w2, pointb); 93 } 94 } 95 else 96 { 97 //x方向進行插值 98 w1 = inter_linear(x[tempi], z[tempi*n + tempj], x[tempi + 1], z[(tempi + 1)*n + tempj], pointa); 99 w2 = inter_linear(x[tempi], z[tempi*n + tempj + 1], x[tempi + 1], z[(tempi + 1)*n + tempj + 1], pointa); 100 /*******y方向進行插值********/ 101 102 if (y[tempj] == pointb) 103 104 { 105 106 //取網格節點值 107 108 w = w1; 109 110 } 111 else 112 { 113 //進行插值(y) 114 115 w = inter_linear(y[tempj], w1, y[tempj + 1], w2, pointb); 116 } 117 } 118 out_result[i*bsize + j] = w; 119 } 120 121 } 122 123 124 } 125 } 126 return 0; 127 }
測試函數：
1 int main() 2 { 3 //插值函數測試 4 int ret = 0; 5 double a1[5] = { 10, 20, 30, 40, 50 }; 6 double a2[6] = {1,2,4,6,8,9}; 7 double z[30] ; 8 for (int i = 0; i < 5; ++i) 9 { 10 for (int j = 0; j < 6; ++j) 11 { 12 z[i*6+j] = a1[i]+a2[j]; 13 } 14 } 15 16 double b1[10] = { 10, 20, 30, 40, 50, 10, 20, 30, 40, 50 }; 17 double b2[18] = { 0, 2, 4, 6, 8, 10, 0, 2, 4, 6, 8, 10, 0, 2, 4, 6, 8, 10 }; 18 double zout[180]; 19 ret = interp2d(a1, a2, z, 5, 6, b1, b2, 10, 18, zout); 20 for (int i = 0; i < 5; ++i) 21 { 22 for (int j = 0; j < 6; j++) 23 { 24 cout << z[i * 6 + j] << " "; 25 } 26 cout << endl; 27 } 28 cout << endl; 29 for (int i = 0; i < 10; ++i) 30 { 31 for (int j = 0; j <18; j++) 32 { 33 cout << zout[i *18 + j] << " "; 34 } 35 cout << endl; 36 } 37 return 0; 38 }