ARTS  2.4.0(git:4fb77825)
special_interp.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
62 #include <cmath>
63 #include <iostream>
64 #include <stdexcept>
65 #include "auto_md.h"
66 #include "check_input.h"
67 #include "math_funcs.h"
68 #include "messages.h"
69 #include "special_interp.h"
70 
71 /*===========================================================================
72  === Point interpolation functions for atmospheric grids and fields
73  ===========================================================================*/
74 
76  const Index& atmosphere_dim,
77  const ArrayOfGridPos& gp_p,
78  const ArrayOfGridPos& gp_lat,
79  const ArrayOfGridPos& gp_lon) {
80  const Index n = gp_p.nelem();
81 
82  if (atmosphere_dim == 1) {
83  itw.resize(n, 2);
84  interpweights(itw, gp_p);
85  }
86 
87  else if (atmosphere_dim == 2) {
88  assert(gp_lat.nelem() == n);
89  itw.resize(n, 4);
90  interpweights(itw, gp_p, gp_lat);
91  }
92 
93  else if (atmosphere_dim == 3) {
94  assert(gp_lat.nelem() == n);
95  assert(gp_lon.nelem() == n);
96  itw.resize(n, 8);
97  interpweights(itw, gp_p, gp_lat, gp_lon);
98  }
99 }
100 
102  const Index& atmosphere_dim,
103  ConstTensor3View x_field,
104  const ArrayOfGridPos& gp_p,
105  const ArrayOfGridPos& gp_lat,
106  const ArrayOfGridPos& gp_lon,
107  ConstMatrixView itw) {
108  assert(x.nelem() == gp_p.nelem());
109 
110  if (atmosphere_dim == 1) {
111  assert(itw.ncols() == 2);
112  interp(x, itw, x_field(Range(joker), 0, 0), gp_p);
113  }
114 
115  else if (atmosphere_dim == 2) {
116  assert(itw.ncols() == 4);
117  interp(x, itw, x_field(Range(joker), Range(joker), 0), gp_p, gp_lat);
118  }
119 
120  else if (atmosphere_dim == 3) {
121  assert(itw.ncols() == 8);
122  interp(x, itw, x_field, gp_p, gp_lat, gp_lon);
123  }
124 }
125 
127  const Index& atmosphere_dim,
128  ConstTensor3View x_field,
129  const ArrayOfGridPos& gp_p,
130  const ArrayOfGridPos& gp_lat,
131  const ArrayOfGridPos& gp_lon) {
132  Matrix itw;
133 
134  interp_atmfield_gp2itw(itw, atmosphere_dim, gp_p, gp_lat, gp_lon);
135 
136  interp_atmfield_by_itw(x, atmosphere_dim, x_field, gp_p, gp_lat, gp_lon, itw);
137 }
138 
140  ConstTensor3View x_field,
141  const GridPos& gp_p,
142  const GridPos& gp_lat,
143  const GridPos& gp_lon) {
144  ArrayOfGridPos agp_p(1), agp_lat(0), agp_lon(0);
145 
146  gridpos_copy(agp_p[0], gp_p);
147 
148  if (atmosphere_dim > 1) {
149  agp_lat.resize(1);
150  gridpos_copy(agp_lat[0], gp_lat);
151  }
152 
153  if (atmosphere_dim > 2) {
154  agp_lon.resize(1);
155  gridpos_copy(agp_lon[0], gp_lon);
156  }
157 
158  Vector x(1);
159 
160  interp_atmfield_by_gp(x, atmosphere_dim, x_field, agp_p, agp_lat, agp_lon);
161 
162  return x[0];
163 }
164 
166  GridPos& gp_p_out,
167  GridPos& gp_lat_out,
168  GridPos& gp_lon_out,
169  const GridPos& gp_p_in,
170  const GridPos& gp_lat_in,
171  const GridPos& gp_lon_in,
172  const Index& atmosphere_dim,
173  const ArrayOfIndex& cloudbox_limits) {
174  // Shift grid positions to cloud box grids
175  if (atmosphere_dim == 1) {
176  gridpos_copy(gp_p_out, gp_p_in);
177  gp_p_out.idx -= cloudbox_limits[0];
179  assert(itw.nelem() == 2);
180  interpweights(itw, gp_p_out);
181  } else if (atmosphere_dim == 2) {
182  gridpos_copy(gp_p_out, gp_p_in);
183  gridpos_copy(gp_lat_out, gp_lat_in);
184  gp_p_out.idx -= cloudbox_limits[0];
185  gp_lat_out.idx -= cloudbox_limits[2];
188  assert(itw.nelem() == 4);
189  interpweights(itw, gp_p_out, gp_lat_out);
190  } else {
191  gridpos_copy(gp_p_out, gp_p_in);
192  gridpos_copy(gp_lat_out, gp_lat_in);
193  gridpos_copy(gp_lon_out, gp_lon_in);
194  gp_p_out.idx -= cloudbox_limits[0];
195  gp_lat_out.idx -= cloudbox_limits[2];
196  gp_lon_out.idx -= cloudbox_limits[4];
200  assert(itw.nelem() == 8);
201  interpweights(itw, gp_p_out, gp_lat_out, gp_lon_out);
202  }
203 }
204 
225  const Index& atmosphere_dim,
226  const ArrayOfGridPos& gp_lat,
227  const ArrayOfGridPos& gp_lon) {
228  if (atmosphere_dim == 1) {
229  itw.resize(1, 1);
230  itw = 1;
231  }
232 
233  else if (atmosphere_dim == 2) {
234  const Index n = gp_lat.nelem();
235  itw.resize(n, 2);
236  interpweights(itw, gp_lat);
237  }
238 
239  else if (atmosphere_dim == 3) {
240  const Index n = gp_lat.nelem();
241  assert(n == gp_lon.nelem());
242  itw.resize(n, 4);
243  interpweights(itw, gp_lat, gp_lon);
244  }
245 }
246 
248  const Index& atmosphere_dim,
249  ConstMatrixView x_surface,
250  const ArrayOfGridPos& gp_lat,
251  const ArrayOfGridPos& gp_lon,
252  ConstMatrixView itw) {
253  if (atmosphere_dim == 1) {
254  assert(itw.ncols() == 1);
255  x = x_surface(0, 0);
256  }
257 
258  else if (atmosphere_dim == 2) {
259  assert(x.nelem() == gp_lat.nelem());
260  assert(itw.ncols() == 2);
261  interp(x, itw, x_surface(Range(joker), 0), gp_lat);
262  }
263 
264  else if (atmosphere_dim == 3) {
265  assert(x.nelem() == gp_lat.nelem());
266  assert(itw.ncols() == 4);
267  interp(x, itw, x_surface, gp_lat, gp_lon);
268  }
269 }
270 
272  const Index& atmosphere_dim,
273  ConstMatrixView x_surface,
274  const ArrayOfGridPos& gp_lat,
275  const ArrayOfGridPos& gp_lon) {
276  Matrix itw;
277 
278  interp_atmsurface_gp2itw(itw, atmosphere_dim, gp_lat, gp_lon);
279 
280  interp_atmsurface_by_itw(x, atmosphere_dim, x_surface, gp_lat, gp_lon, itw);
281 }
282 
284  ConstMatrixView x_surface,
285  const GridPos& gp_lat,
286  const GridPos& gp_lon) {
287  ArrayOfGridPos agp_lat(0), agp_lon(0);
288 
289  if (atmosphere_dim > 1) {
290  agp_lat.resize(1);
291  gridpos_copy(agp_lat[0], gp_lat);
292  }
293 
294  if (atmosphere_dim > 2) {
295  agp_lon.resize(1);
296  gridpos_copy(agp_lon[0], gp_lon);
297  }
298 
299  Vector x(1);
300 
301  interp_atmsurface_by_gp(x, atmosphere_dim, x_surface, agp_lat, agp_lon);
302 
303  return x[0];
304 }
305 
306 /*===========================================================================
307  === Regridding
308  ===========================================================================*/
309 
311  const Index& atmosphere_dim,
312  ConstTensor3View field_old,
313  const ArrayOfGridPos& gp_p,
314  const ArrayOfGridPos& gp_lat,
315  const ArrayOfGridPos& gp_lon) {
316  const Index n1 = gp_p.nelem();
317 
318  if (atmosphere_dim == 1) {
319  field_new.resize(n1, 1, 1);
320  Matrix itw(n1, 2);
321  interpweights(itw, gp_p);
322  interp(field_new(joker, 0, 0), itw, field_old(joker, 0, 0), gp_p);
323  } else if (atmosphere_dim == 2) {
324  const Index n2 = gp_lat.nelem();
325  field_new.resize(n1, n2, 1);
326  Tensor3 itw(n1, n2, 4);
327  interpweights(itw, gp_p, gp_lat);
328  interp(field_new(joker, joker, 0),
329  itw,
330  field_old(joker, joker, 0),
331  gp_p,
332  gp_lat);
333  } else if (atmosphere_dim == 3) {
334  const Index n2 = gp_lat.nelem();
335  const Index n3 = gp_lon.nelem();
336  field_new.resize(n1, n2, n3);
337  Tensor4 itw(n1, n2, n3, 8);
338  interpweights(itw, gp_p, gp_lat, gp_lon);
339  interp(field_new, itw, field_old, gp_p, gp_lat, gp_lon);
340  }
341 }
342 
343 void regrid_atmsurf_by_gp(Matrix& field_new,
344  const Index& atmosphere_dim,
345  ConstMatrixView field_old,
346  const ArrayOfGridPos& gp_lat,
347  const ArrayOfGridPos& gp_lon) {
348  if (atmosphere_dim == 1) {
349  field_new = field_old;
350  } else if (atmosphere_dim == 2) {
351  const Index n1 = gp_lat.nelem();
352  field_new.resize(n1, 1);
353  Matrix itw(n1, 2);
354  interpweights(itw, gp_lat);
355  interp(field_new(joker, 0), itw, field_old(joker, 0), gp_lat);
356  } else if (atmosphere_dim == 3) {
357  const Index n1 = gp_lat.nelem();
358  const Index n2 = gp_lon.nelem();
359  field_new.resize(n1, n2);
360  Tensor3 itw(n1, n2, 4);
361  interpweights(itw, gp_lat, gp_lon);
362  interp(field_new, itw, field_old, gp_lat, gp_lon);
363  }
364 }
365 
367  ArrayOfGridPos& gp_lat,
368  ArrayOfGridPos& gp_lon,
369  const RetrievalQuantity& rq,
370  const Index& atmosphere_dim,
371  const Vector& p_grid,
372  const Vector& lat_grid,
373  const Vector& lon_grid) {
374  gp_p.resize(rq.Grids()[0].nelem());
375  p2gridpos(gp_p, p_grid, rq.Grids()[0], 0);
376  //
377  if (atmosphere_dim >= 2) {
378  gp_lat.resize(rq.Grids()[1].nelem());
379  gridpos(gp_lat, lat_grid, rq.Grids()[1], 0);
380  } else {
381  gp_lat.resize(0);
382  }
383  //
384  if (atmosphere_dim >= 3) {
385  gp_lon.resize(rq.Grids()[2].nelem());
386  gridpos(gp_lon, lon_grid, rq.Grids()[2], 0);
387  } else {
388  gp_lon.resize(0);
389  }
390 }
391 
393  ArrayOfGridPos& gp_lon,
394  const RetrievalQuantity& rq,
395  const Index& atmosphere_dim,
396  const Vector& lat_grid,
397  const Vector& lon_grid) {
398  if (atmosphere_dim >= 2) {
399  gp_lat.resize(rq.Grids()[0].nelem());
400  gridpos(gp_lat, lat_grid, rq.Grids()[0], 0);
401  } else {
402  gp_lat.resize(0);
403  }
404  //
405  if (atmosphere_dim >= 3) {
406  gp_lon.resize(rq.Grids()[1].nelem());
407  gridpos(gp_lon, lon_grid, rq.Grids()[1], 0);
408  } else {
409  gp_lon.resize(0);
410  }
411 }
412 
414  ArrayOfGridPos& gp_lat,
415  ArrayOfGridPos& gp_lon,
416  Index& n_p,
417  Index& n_lat,
418  Index& n_lon,
419  const ArrayOfVector& ret_grids,
420  const Index& atmosphere_dim,
421  const Vector& p_grid,
422  const Vector& lat_grid,
423  const Vector& lon_grid) {
424  // We want here an extrapolation to infinity ->
425  // extremly high extrapolation factor
426  const Numeric inf_proxy = 1.0e99;
427 
428  gp_p.resize(p_grid.nelem());
429  n_p = ret_grids[0].nelem();
430  if (n_p > 1) {
431  p2gridpos(gp_p, ret_grids[0], p_grid, inf_proxy);
433  } else {
434  gp4length1grid(gp_p);
435  }
436 
437  if (atmosphere_dim >= 2) {
438  gp_lat.resize(lat_grid.nelem());
439  n_lat = ret_grids[1].nelem();
440  if (n_lat > 1) {
441  gridpos(gp_lat, ret_grids[1], lat_grid, inf_proxy);
442  jacobian_type_extrapol(gp_lat);
443  } else {
444  gp4length1grid(gp_lat);
445  }
446  } else {
447  gp_lat.resize(0);
448  n_lat = 1;
449  }
450  //
451  if (atmosphere_dim >= 3) {
452  gp_lon.resize(lon_grid.nelem());
453  n_lon = ret_grids[2].nelem();
454  if (n_lon > 1) {
455  gridpos(gp_lon, ret_grids[2], lon_grid, inf_proxy);
456  jacobian_type_extrapol(gp_lon);
457  } else {
458  gp4length1grid(gp_lon);
459  }
460  } else {
461  gp_lon.resize(0);
462  n_lon = 1;
463  }
464 }
465 
467  ArrayOfGridPos& gp_lon,
468  Index& n_lat,
469  Index& n_lon,
470  const ArrayOfVector& ret_grids,
471  const Index& atmosphere_dim,
472  const Vector& lat_grid,
473  const Vector& lon_grid) {
474  // We want here an extrapolation to infinity ->
475  // extremly high extrapolation factor
476  const Numeric inf_proxy = 1.0e99;
477 
478  if (atmosphere_dim >= 2) {
479  gp_lat.resize(lat_grid.nelem());
480  n_lat = ret_grids[0].nelem();
481  if (n_lat > 1) {
482  gridpos(gp_lat, ret_grids[0], lat_grid, inf_proxy);
483  jacobian_type_extrapol(gp_lat);
484  } else {
485  gp4length1grid(gp_lat);
486  }
487  } else {
488  gp_lat.resize(0);
489  n_lat = 1;
490  }
491  //
492  if (atmosphere_dim >= 3) {
493  gp_lon.resize(lon_grid.nelem());
494  n_lon = ret_grids[1].nelem();
495  if (n_lon > 1) {
496  gridpos(gp_lon, ret_grids[1], lon_grid, inf_proxy);
497  jacobian_type_extrapol(gp_lon);
498  } else {
499  gp4length1grid(gp_lon);
500  }
501  } else {
502  gp_lon.resize(0);
503  n_lon = 1;
504  }
505 }
506 
508  const Index& atmosphere_dim,
509  ConstTensor3View field_old,
510  const ArrayOfGridPos& gp_p,
511  const ArrayOfGridPos& gp_lat,
512  const ArrayOfGridPos& gp_lon) {
513  const Index n1 = gp_p.nelem();
514 
515  const bool np_is1 = field_old.npages() == 1 ? true : false;
516  const bool nlat_is1 =
517  atmosphere_dim > 1 && field_old.nrows() == 1 ? true : false;
518  const bool nlon_is1 =
519  atmosphere_dim > 2 && field_old.ncols() == 1 ? true : false;
520 
521  // If no length 1, we can use standard function
522  if (!np_is1 && !nlat_is1 && !nlon_is1) {
524  field_new, atmosphere_dim, field_old, gp_p, gp_lat, gp_lon);
525  } else {
526  //--- 1D (1 possibilities left) -------------------------------------------
527  if (atmosphere_dim == 1) { // 1: No interpolation at all
528  field_new.resize(n1, 1, 1);
529  field_new(joker, 0, 0) = field_old(0, 0, 0);
530  }
531 
532  //--- 2D (3 possibilities left) -------------------------------------------
533  else if (atmosphere_dim == 2) {
534  const Index n2 = gp_lat.nelem();
535  field_new.resize(n1, n2, 1);
536  //
537  if (np_is1 && nlat_is1) // 1: No interpolation at all
538  {
539  // Here we need no interpolation at all
540  field_new(joker, joker, 0) = field_old(0, 0, 0);
541  } else if (np_is1) // 2: Latitude interpolation
542  {
543  Matrix itw(n2, 2);
544  interpweights(itw, gp_lat);
545  Vector tmp(n2);
546  interp(tmp, itw, field_old(0, joker, 0), gp_lat);
547  for (Index p = 0; p < n1; p++) {
548  assert(gp_p[p].fd[0] < 1e-6);
549  field_new(p, joker, 0) = tmp;
550  }
551  } else // 3: Pressure interpolation
552  {
553  Matrix itw(n1, 2);
554  interpweights(itw, gp_p);
555  Vector tmp(n1);
556  interp(tmp, itw, field_old(joker, 0, 0), gp_p);
557  for (Index lat = 0; lat < n2; lat++) {
558  assert(gp_lat[lat].fd[0] < 1e-6);
559  field_new(joker, lat, 0) = tmp;
560  }
561  }
562  }
563 
564  //--- 3D (7 possibilities left) -------------------------------------------
565  else if (atmosphere_dim == 3) {
566  const Index n2 = gp_lat.nelem();
567  const Index n3 = gp_lon.nelem();
568  field_new.resize(n1, n2, n3);
569  //
570  if (np_is1 && nlat_is1 && nlon_is1) // 1: No interpolation at all
571  {
572  field_new(joker, joker, joker) = field_old(0, 0, 0);
573  }
574 
575  else if (np_is1) // No pressure interpolation --------------
576  {
577  if (nlat_is1) // 2: Just longitude interpolation
578  {
579  Matrix itw(n3, 2);
580  interpweights(itw, gp_lon);
581  Vector tmp(n3);
582  interp(tmp, itw, field_old(0, 0, joker), gp_lon);
583  for (Index p = 0; p < n1; p++) {
584  assert(gp_p[p].fd[0] < 1e-6);
585  for (Index lat = 0; lat < n2; lat++) {
586  assert(gp_lat[lat].fd[0] < 1e-6);
587  field_new(p, lat, joker) = tmp;
588  }
589  }
590  } else if (nlon_is1) // 3: Just latitude interpolation
591  {
592  Matrix itw(n2, 2);
593  interpweights(itw, gp_lat);
594  Vector tmp(n2);
595  interp(tmp, itw, field_old(0, joker, 0), gp_lat);
596  for (Index p = 0; p < n1; p++) {
597  assert(gp_p[p].fd[0] < 1e-6);
598  for (Index lon = 0; lon < n3; lon++) {
599  assert(gp_lon[lon].fd[0] < 1e-6);
600  field_new(p, joker, lon) = tmp;
601  }
602  }
603  } else // 4: Both lat and lon interpolation
604  {
605  Tensor3 itw(n2, n3, 4);
606  interpweights(itw, gp_lat, gp_lon);
607  Matrix tmp(n2, n3);
608  interp(tmp, itw, field_old(0, joker, joker), gp_lat, gp_lon);
609  for (Index p = 0; p < n1; p++) {
610  assert(gp_p[p].fd[0] < 1e-6);
611  field_new(p, joker, joker) = tmp;
612  }
613  }
614  }
615 
616  else // Pressure interpolation --------------
617  {
618  if (nlat_is1 && nlon_is1) // 5: Just pressure interpolatiom
619  {
620  Matrix itw(n1, 2);
621  interpweights(itw, gp_p);
622  Vector tmp(n1);
623  interp(tmp, itw, field_old(joker, 0, 0), gp_p);
624  for (Index lat = 0; lat < n2; lat++) {
625  assert(gp_lat[lat].fd[0] < 1e-6);
626  for (Index lon = 0; lon < n3; lon++) {
627  assert(gp_lon[lon].fd[0] < 1e-6);
628  field_new(joker, lat, lon) = tmp;
629  }
630  }
631  } else if (nlat_is1) // 6: Both p and lon interpolation
632  {
633  Tensor3 itw(n1, n3, 4);
634  interpweights(itw, gp_p, gp_lon);
635  Matrix tmp(n1, n3);
636  interp(tmp, itw, field_old(joker, 0, joker), gp_p, gp_lon);
637  for (Index lat = 0; lat < n2; lat++) {
638  assert(gp_lat[lat].fd[0] < 1e-6);
639  field_new(joker, lat, joker) = tmp;
640  }
641  } else // 7: Both p and lat interpolation
642  {
643  Tensor3 itw(n1, n2, 4);
644  interpweights(itw, gp_p, gp_lat);
645  Matrix tmp(n1, n2);
646  interp(tmp, itw, field_old(joker, joker, 0), gp_p, gp_lat);
647  for (Index lon = 0; lon < n3; lon++) {
648  assert(gp_lon[lon].fd[0] < 1e-6);
649  field_new(joker, joker, lon) = tmp;
650  }
651  }
652  }
653  }
654  }
655 }
656 
657 /* So far just a temporary test */
659  const Index& atmosphere_dim,
660  ConstMatrixView field_old,
661  const ArrayOfGridPos& gp_lat,
662  const ArrayOfGridPos& gp_lon) {
663  // As 1D is so simple, let's do it here and not go to standard function
664  if (atmosphere_dim == 1) {
665  field_new = field_old;
666  } else {
667  const bool nlat_is1 = field_old.nrows() == 1 ? true : false;
668  const bool nlon_is1 =
669  atmosphere_dim > 2 && field_old.ncols() == 1 ? true : false;
670 
671  // If no length 1, we can use standard function
672  if (!nlat_is1 && !nlon_is1) {
674  field_new, atmosphere_dim, field_old, gp_lat, gp_lon);
675  } else {
676  if (atmosphere_dim == 2) { // 1: No interpolation at all
677  const Index n1 = gp_lat.nelem();
678  field_new.resize(n1, 1);
679  field_new(joker, 0) = field_old(0, 0);
680  } else {
681  const Index n1 = gp_lat.nelem();
682  const Index n2 = gp_lon.nelem();
683  field_new.resize(n1, n2);
684  //
685  if (nlat_is1 && nlon_is1) // 1: No interpolation at all
686  {
687  field_new(joker, joker) = field_old(0, 0);
688  } else if (nlon_is1) // 2: Just latitude interpolation
689  {
690  Matrix itw(n1, 2);
691  interpweights(itw, gp_lat);
692  Vector tmp(n1);
693  interp(tmp, itw, field_old(joker, 0), gp_lat);
694  for (Index lon = 0; lon < n2; lon++) {
695  assert(gp_lon[lon].fd[0] < 1e-6);
696  field_new(joker, lon) = tmp;
697  }
698  } else // 2: Just longitude interpolation
699  {
700  Matrix itw(n2, 2);
701  interpweights(itw, gp_lon);
702  Vector tmp(n2);
703  interp(tmp, itw, field_old(0, joker), gp_lon);
704  for (Index lat = 0; lat < n1; lat++) {
705  assert(gp_lat[lat].fd[0] < 1e-6);
706  field_new(lat, joker) = tmp;
707  }
708  }
709  }
710  }
711  }
712 }
713 
714 /*===========================================================================
715  === Conversion altitudes / pressure
716  ===========================================================================*/
717 
718 void itw2p(VectorView p_values,
720  const ArrayOfGridPos& gp,
721  ConstMatrixView itw) {
722  assert(itw.ncols() == 2);
723  assert(p_values.nelem() == itw.nrows());
724 
725  // Local variable to store log of the pressure grid:
726  Vector logpgrid(p_grid.nelem());
727 
728  transform(logpgrid, log, p_grid);
729 
730  interp(p_values, itw, logpgrid, gp);
731 
732  transform(p_values, exp, p_values);
733 }
734 
760  ConstVectorView old_pgrid,
761  ConstVectorView new_pgrid,
762  const Numeric& extpolfac) {
763  // Local variable to store log of the pressure grids
764  Vector logold(old_pgrid.nelem());
765  Vector lognew(new_pgrid.nelem());
766 
767  transform(logold, log, old_pgrid);
768  transform(lognew, log, new_pgrid);
769 
770  gridpos(gp, logold, lognew, extpolfac);
771 }
772 
774  ConstVectorView old_pgrid,
775  ConstVectorView new_pgrid,
776  const Index order,
777  const Numeric& extpolfac) {
778  // Local variable to store log of the pressure grids
779  Vector logold(old_pgrid.nelem());
780  Vector lognew(new_pgrid.nelem());
781 
782  transform(logold, log, old_pgrid);
783  transform(lognew, log, new_pgrid);
784 
785  gridpos_poly(gp, logold, lognew, order, extpolfac);
786 }
787 
789  GridPos& gp_lat,
790  GridPos& gp_lon,
791  const Index& atmosphere_dim,
798 
799  if (atmosphere_dim == 1) {
801  "Altitude interpolation", z_field(joker, 0, 0), rte_pos[0]);
802  gridpos(gp_p, z_field(joker, 0, 0), rte_pos[0]);
803  } else {
804  // Determine z at lat/lon (z_grid) by blue interpolation
805  const Index np = p_grid.nelem();
806  Vector z_grid(np);
807  ArrayOfGridPos agp_z, agp_lat(np);
808  //
809  gridpos_1to1(agp_z, p_grid);
810  //
811  chk_interpolation_grids("Latitude interpolation", lat_grid, rte_pos[1]);
812  gridpos(gp_lat, lat_grid, rte_pos[1]);
813 
814  if (atmosphere_dim == 2) {
815  for (Index i = 0; i < np; i++) {
816  agp_lat[i] = gp_lat;
817  }
818  Matrix itw(np, 4);
819  interpweights(itw, agp_z, agp_lat);
820  interp(z_grid, itw, z_field(joker, joker, 0), agp_z, agp_lat);
821  } else {
822  chk_interpolation_grids("Longitude interpolation", lon_grid, rte_pos[2]);
823  gridpos(gp_lon, lon_grid, rte_pos[2]);
824  ArrayOfGridPos agp_lon(np);
825  for (Index i = 0; i < np; i++) {
826  agp_lat[i] = gp_lat;
827  agp_lon[i] = gp_lon;
828  }
829  Matrix itw(np, 8);
830  interpweights(itw, agp_z, agp_lat, agp_lon);
831  interp(z_grid, itw, z_field, agp_z, agp_lat, agp_lon);
832  }
833 
834  // And use z_grid to get gp_p (gp_al and gp_lon determined above)
835  chk_interpolation_grids("Altitude interpolation", z_grid, rte_pos[0]);
836  gridpos(gp_p, z_grid, rte_pos[0]);
837  }
838 }
839 
841  GridPos& gp_lon,
842  const Index& atmosphere_dim,
847 
848  if (atmosphere_dim == 1) {
849  } else {
850  chk_interpolation_grids("Latitude interpolation", lat_grid, rte_pos[1]);
851  gridpos(gp_lat, lat_grid, rte_pos[1]);
852 
853  if (atmosphere_dim == 3) {
854  chk_interpolation_grids("Longitude interpolation", lon_grid, rte_pos[2]);
855  gridpos(gp_lon, lon_grid, rte_pos[2]);
856  }
857  }
858 }
859 
862 // FIXME only used in assertion
863 #ifndef NDEBUG
865 #else
867 #endif
869  const GridPos& gp_lat) {
870  const Index np = p_grid.nelem();
871 
872  assert(z.nelem() == np);
873  assert(z_field.nrows() == np);
874  assert(z_field.ncols() == lat_grid.nelem());
875 
876  Matrix z_matrix(np, 1);
877  ArrayOfGridPos gp_z(np), agp_lat(1);
878  Tensor3 itw(np, 1, 4);
879 
880  gridpos_copy(agp_lat[0], gp_lat);
881  gridpos(gp_z, p_grid, p_grid);
882  interpweights(itw, gp_z, agp_lat);
883  interp(z_matrix, itw, z_field, gp_z, agp_lat);
884 
885  z = z_matrix(Range(joker), 0);
886 }
887 
890 //FIXME only used in assertion
891 #ifndef NDEBUG
894 #else
897 #endif
899  const GridPos& gp_lat,
900  const GridPos& gp_lon) {
901  const Index np = p_grid.nelem();
902 
903  assert(z.nelem() == np);
904  assert(z_field.npages() == np);
905  assert(z_field.nrows() == lat_grid.nelem());
906  assert(z_field.ncols() == lon_grid.nelem());
907 
908  Tensor3 z_tensor(np, 1, 1);
909  ArrayOfGridPos agp_z(np), agp_lat(1), agp_lon(1);
910  Tensor4 itw(np, 1, 1, 8);
911 
912  gridpos_copy(agp_lat[0], gp_lat);
913  gridpos_copy(agp_lon[0], gp_lon);
914  gridpos(agp_z, p_grid, p_grid);
915  interpweights(itw, agp_z, agp_lat, agp_lon);
916 
917  interp(z_tensor, itw, z_field, agp_z, agp_lat, agp_lon);
918 
919  z = z_tensor(Range(joker), 0, 0);
920 }
921 
922 /*===========================================================================
923  === Interpolation of GriddedFields
924  ===========================================================================*/
925 
927  MatrixView n_imag,
928  const GriddedField3& complex_n,
929  const String& varname,
931  ConstVectorView t_grid) {
932  // Set expected order of grids
933  Index gfield_fID = 0;
934  Index gfield_tID = 1;
935  Index gfield_compID = 2;
936 
937  // Check of complex_n
938  //
939  complex_n.checksize_strict();
940  //
941  chk_griddedfield_gridname(complex_n, gfield_fID, "Frequency");
942  chk_griddedfield_gridname(complex_n, gfield_tID, "Temperature");
943  chk_griddedfield_gridname(complex_n, gfield_compID, "Complex");
944  //
945  if (complex_n.data.ncols() != 2) {
946  ostringstream os;
947  os << "The data in *" << varname
948  << "* must have exactly two pages. One page "
949  << "each\nfor the real and imaginary part of the complex refractive index.";
950  }
951 
952  // Frequency and temperature grid sizes
953  const Index nf_in = complex_n.data.npages();
954  const Index nt_in = complex_n.data.nrows();
955  const Index nf_out = f_grid.nelem();
956  const Index nt_out = t_grid.nelem();
957 
958  //Assert size of output variables
959  assert(n_real.nrows() == nf_out && n_real.ncols() == nt_out);
960  assert(n_imag.nrows() == nf_out && n_imag.ncols() == nt_out);
961 
962  const Vector& f_grid_in = complex_n.get_numeric_grid(gfield_fID);
963  const Vector& t_grid_in = complex_n.get_numeric_grid(gfield_tID);
964 
965  // Expand/interpolate in frequency dimension
966  //
967  Matrix nrf(nf_out, nt_in), nif(nf_out, nt_in);
968  //
969  if (nf_in == 1) {
970  for (Index i = 0; i < nf_out; i++) {
971  nrf(i, joker) = complex_n.data(0, joker, 0);
972  nif(i, joker) = complex_n.data(0, joker, 1);
973  }
974  } else {
975  chk_interpolation_grids("Frequency interpolation", f_grid_in, f_grid);
976  //
977  ArrayOfGridPos gp(nf_out);
978  Matrix itw(nf_out, 2);
979  gridpos(gp, f_grid_in, f_grid);
980  interpweights(itw, gp);
981  for (Index i = 0; i < nt_in; i++) {
982  interp(nrf(joker, i), itw, complex_n.data(joker, i, 0), gp);
983  interp(nif(joker, i), itw, complex_n.data(joker, i, 1), gp);
984  }
985  }
986 
987  // Expand/interpolate in temperature dimension
988  //
989  if (nt_in == 1) {
990  for (Index i = 0; i < nt_out; i++) {
991  n_real(joker, i) = nrf(joker, 0);
992  n_imag(joker, i) = nif(joker, 0);
993  }
994  } else {
995  chk_interpolation_grids("Temperature interpolation", t_grid_in, t_grid);
996  //
997  ArrayOfGridPos gp(nt_out);
998  Matrix itw(nt_out, 2);
999  gridpos(gp, t_grid_in, t_grid);
1000  interpweights(itw, gp);
1001  for (Index i = 0; i < nf_out; i++) {
1002  interp(n_real(i, joker), itw, nrf(i, joker), gp);
1003  interp(n_imag(i, joker), itw, nif(i, joker), gp);
1004  }
1005  }
1006 }
gridpos_1to1
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
Definition: interpolation.cc:449
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
transform
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1476
interp_atmsurface_by_gp
void interp_atmsurface_by_gp(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates a surface-type variable given the grid positions.
Definition: special_interp.cc:271
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
ARTS::Var::lat
Numeric lat(Workspace &ws) noexcept
Definition: autoarts.h:3933
auto_md.h
gridpos_copy
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Definition: interpolation.cc:473
RetrievalQuantity::Grids
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:255
complex_n_interp
void complex_n_interp(MatrixView n_real, MatrixView n_imag, const GriddedField3 &complex_n, const String &varname, ConstVectorView f_grid, ConstVectorView t_grid)
General function for interpolating data of complex n type.
Definition: special_interp.cc:926
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
itw2p
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Definition: special_interp.cc:718
p2gridpos
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Definition: special_interp.cc:759
regrid_atmfield_by_gp
void regrid_atmfield_by_gp(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric field, for precalculated grid positions.
Definition: special_interp.cc:310
jacobian_type_extrapol
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:885
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
regrid_atmsurf_by_gp
void regrid_atmsurf_by_gp(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric surface, for precalculated grid positions.
Definition: special_interp.cc:343
z_at_latlon
void z_at_latlon(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, const GridPos &gp_lat, const GridPos &gp_lon)
Returns the geomtrical altitudes of p_grid for one latitude and one longitude.
Definition: special_interp.cc:888
interp_atmfield_by_itw
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Definition: special_interp.cc:101
chk_interpolation_grids
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:989
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
regrid_atmsurf_by_gp_oem
void regrid_atmsurf_by_gp_oem(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of surface field OEM-type.
Definition: special_interp.cc:658
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
z_at_lat_2d
void z_at_lat_2d(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstMatrixView z_field, const GridPos &gp_lat)
Returns the geomtrical altitudes of p_grid for one latitude.
Definition: special_interp.cc:860
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
ARTS::Var::lon
Numeric lon(Workspace &ws) noexcept
Definition: autoarts.h:4064
GriddedField3
Definition: gridded_fields.h:284
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
interp_atmsurface_gp2itw
void interp_atmsurface_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of a surface-type variable.
Definition: special_interp.cc:224
GriddedField3::checksize_strict
void checksize_strict() const final
Strict consistency check.
Definition: gridded_fields.h:307
get_gp_atmgrids_to_rq
void get_gp_atmgrids_to_rq(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Definition: special_interp.cc:366
Array
This can be used to make arrays out of anything.
Definition: array.h:108
GriddedField3::data
Tensor3 data
Definition: gridded_fields.h:337
GriddedField::get_numeric_grid
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
Definition: gridded_fields.cc:87
messages.h
Declarations having to do with the four output streams.
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
chk_griddedfield_gridname
void chk_griddedfield_gridname(const GriddedField &gf, const Index gridindex, const String &gridname)
Check name of grid in GriddedField.
Definition: check_input.cc:1891
gridpos_upperend_check
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
Definition: interpolation.cc:589
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
my_basic_string< char >
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
math_funcs.h
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
chk_rte_pos
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
Definition: check_input.cc:1770
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
Range
The range class.
Definition: matpackI.h:160
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
interp_atmfield_gp2itw
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
Definition: special_interp.cc:75
gp4length1grid
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
Definition: interpolation.cc:636
interp_atmfield_by_gp
void interp_atmfield_by_gp(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates an atmospheric field given the grid positions.
Definition: special_interp.cc:126
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
special_interp.h
Header file for special_interp.cc.
rte_pos2gridpos
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
Definition: special_interp.cc:788
GridPos::idx
Index idx
Definition: interpolation.h:74
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
p2gridpos_poly
void p2gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Index order, const Numeric &extpolfac)
p2gridpos_poly
Definition: special_interp.cc:773
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
RetrievalQuantity
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:120
interp_atmsurface_by_itw
void interp_atmsurface_by_itw(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates a surface-type variable with pre-calculated weights by interp_atmsurface_gp2itw.
Definition: special_interp.cc:247
get_gp_rq_to_atmgrids
void get_gp_rq_to_atmgrids(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, Index &n_p, Index &n_lat, Index &n_lon, const ArrayOfVector &ret_grids, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Definition: special_interp.cc:413
check_input.h
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
Vector
The Vector class.
Definition: matpackI.h:860
get_gp_atmsurf_to_rq
void get_gp_atmsurf_to_rq(ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric surfaces to retrieval grids.
Definition: special_interp.cc:392
interp_cloudfield_gp2itw
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
Definition: special_interp.cc:165
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
regrid_atmfield_by_gp_oem
void regrid_atmfield_by_gp_oem(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of atmospheric field OEM-type.
Definition: special_interp.cc:507