grid.c 40.6 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

5
  Copyright (C) 2003-2016 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

Uwe Schulzweida's avatar
Uwe Schulzweida committed
17
#if defined(HAVE_CONFIG_H)
18
19
20
#  include "config.h"
#endif

21
22
23
24
#if defined(_OPENMP)
#  include <omp.h>
#endif

Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
#include <stdio.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
#include <stdarg.h> /* va_list */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27

Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
#if defined(HAVE_LIBPROJ)
29
#  include "proj_api.h"
30
31
#endif

Ralf Mueller's avatar
Ralf Mueller committed
32
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
34
35
36
37
38
#include "cdo.h"
#include "cdo_int.h"
#include "error.h"
#include "grid.h"


39
static
40
void scale_vec(double scalefactor, long nvals, double *restrict values)
41
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
#if defined(_OPENMP)
43
#pragma omp parallel for default(none) shared(nvals, scalefactor, values)
44
#endif
45
  for ( long n = 0; n < nvals; ++n )
46
    {
47
      values[n] *= scalefactor;
48
49
50
51
    }
}


52
void grid_to_radian(const char *units, long nvals, double *restrict values, const char *description)
53
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
  if ( cmpstr(units, "degree") == 0 )
55
    {
56
      scale_vec(DEG2RAD, nvals, values);
57
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
  else if ( cmpstr(units, "radian") == 0 )
59
60
61
62
63
64
65
66
67
    {
      /* No conversion necessary */
    }
  else
    {
      cdoWarning("Unknown units [%s] supplied for %s; proceeding assuming radians!", units, description);
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
68

69
70
void grid_to_degree(const char *units, long nvals, double *restrict values, const char *description)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
  if ( cmpstr(units, "radian") == 0 )
72
73
74
    {
      for ( long n = 0; n < nvals; ++n ) values[n] *= RAD2DEG;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
  else if ( cmpstr(units, "degree") == 0 )
76
77
78
79
80
81
82
83
84
85
    {
      /* No conversion necessary */
    }
  else
    {
      cdoWarning("Unknown units [%s] supplied for %s; proceeding assuming degress!", units, description);
    }
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
88
89
90
91
92
93
94
95
96
int gridToZonal(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
  double  xval = 0;
  double *yvals;

  gridtype = gridInqType(gridID1);
  gridsize = gridInqYsize(gridID1);
  gridID2  = gridCreate(gridtype, gridsize);
	  
97
98
99
  if ( gridtype == GRID_LONLAT   ||
       gridtype == GRID_GAUSSIAN ||
       gridtype == GRID_GENERIC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
102
103
104
105
106
107
    {
      gridDefXsize(gridID2, 1);
      gridDefYsize(gridID2, gridsize);

      gridDefXvals(gridID2, &xval);

      if ( gridInqYvals(gridID1, NULL) )
	{
108
	  yvals = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
110
111
112

	  gridInqYvals(gridID1, yvals);
	  gridDefYvals(gridID2, yvals);

113
	  Free(yvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
	}
    }
116
117
  else
    {
118
      Error("Gridtype %s unsupported!", gridNamePtr(gridtype));
119
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120

121
  return gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
}


int gridToMeridional(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
  double *xvals;
  double  yval = 0;

  gridtype = gridInqType(gridID1);
  gridsize = gridInqXsize(gridID1);
  gridID2  = gridCreate(gridtype, gridsize);
	  
136
137
138
  if ( gridtype == GRID_LONLAT   ||
       gridtype == GRID_GAUSSIAN ||
       gridtype == GRID_GENERIC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
    {
140
141
      gridDefXsize(gridID2, gridsize);
      gridDefYsize(gridID2, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142

143
144
      if ( gridInqXvals(gridID1, NULL) )
	{
145
	  xvals = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146

147
148
	  gridInqXvals(gridID1, xvals);
	  gridDefXvals(gridID2, xvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149

150
	  Free(xvals);
151
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152

153
154
155
156
      gridDefYvals(gridID2, &yval);
    }
  else
    {
157
      Error("Gridtype %s unsupported!", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
159
    }

160
  return gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162
163
}


164
void grid_gen_corners(int n, const double* restrict vals, double* restrict corners)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
{
166
  int i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167

168
  if ( n == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
    {
170
171
      corners[0] = vals[0];
      corners[1] = vals[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
    }
173
174
175
176
  else
    {
      for ( i = 0; i < n-1; ++i )
	corners[i+1] = 0.5*(vals[i] + vals[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177

178
179
180
      corners[0] = 2*vals[0] - corners[1];
      corners[n] = 2*vals[n-1] - corners[n-1];
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
}


184
void grid_gen_bounds(int n, const double *restrict vals, double *restrict bounds)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
{
186
  for ( int i = 0; i < n-1; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
    {
188
189
      bounds[2*i+1]   = 0.5*(vals[i] + vals[i+1]);
      bounds[2*(i+1)] = 0.5*(vals[i] + vals[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
    }

192
193
  bounds[0]     = 2*vals[0] - bounds[1];
  bounds[2*n-1] = 2*vals[n-1] - bounds[2*(n-1)];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
196
}


197
void grid_check_lat_borders(int n, double *ybounds)
198
{
199
  if ( ybounds[0] > ybounds[n-1] )
200
    {
201
202
203
204
205
206
207
      if ( ybounds[0]   >  88 ) ybounds[0]   =  90;
      if ( ybounds[n-1] < -88 ) ybounds[n-1] = -90;
    }
  else
    {
      if ( ybounds[0]   < -88 ) ybounds[0]   = -90;
      if ( ybounds[n-1] >  88 ) ybounds[n-1] =  90;
208
209
210
211
    }
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
void grid_cell_center_to_bounds_X2D(const char* xunitstr, long xsize, long ysize, const double* restrict grid_center_lon, 
				    double* restrict grid_corner_lon, double dlon)
{
  long i, j, index;
  double minlon, maxlon;

  if ( ! (dlon > 0) ) dlon = 360./xsize;
  /*
  if ( xsize == 1 || (grid_center_lon[xsize-1]-grid_center_lon[0]+dlon) < 359 )
    cdoAbort("Cannot calculate Xbounds for %d vals with dlon = %g", xsize, dlon);
  */
  for ( i = 0; i < xsize; ++i )
    {
      minlon = grid_center_lon[i] - 0.5*dlon;
      maxlon = grid_center_lon[i] + 0.5*dlon;
      for ( j = 0; j < ysize; ++j )
	{
	  index = (j<<2)*xsize + (i<<2);
	  grid_corner_lon[index  ] = minlon;
	  grid_corner_lon[index+1] = maxlon;
	  grid_corner_lon[index+2] = maxlon;
	  grid_corner_lon[index+3] = minlon;
	}
    }
}

static
double genYmin(double y1, double y2)
{
  double ymin, dy;

  dy = y2 - y1;
  ymin = y1 - dy/2;

  if ( y1 < -85 && ymin < -87.5 ) ymin = -90;

  if ( cdoVerbose )
    cdoPrint("genYmin: y1 = %g  y2 = %g  dy = %g  ymin = %g", y1, y2, dy, ymin);

251
  return ymin;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
}

static
double genYmax(double y1, double y2)
{
  double ymax, dy;

  dy = y1 - y2;
  ymax = y1 + dy/2;

  if ( y1 > 85 && ymax > 87.5 ) ymax = 90;

  if ( cdoVerbose )
    cdoPrint("genYmax: y1 = %g  y2 = %g  dy = %g  ymax = %g", y1, y2, dy, ymax);

267
  return ymax;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
}



/*****************************************************************************/

void grid_cell_center_to_bounds_Y2D(const char* yunitstr, long xsize, long ysize, const double* restrict grid_center_lat, double* restrict grid_corner_lat)
{
  long i, j, index;
  double minlat, maxlat;
  double firstlat, lastlat;

  firstlat = grid_center_lat[0];
  lastlat  = grid_center_lat[xsize*ysize-1];

  // if ( ysize == 1 ) cdoAbort("Cannot calculate Ybounds for 1 value!");

  for ( j = 0; j < ysize; ++j )
    {
      if ( ysize == 1 )
	{
	  minlat = grid_center_lat[0] - 360./ysize;
	  maxlat = grid_center_lat[0] + 360./ysize;
	}
      else
	{
	  index = j*xsize;
	  if ( firstlat > lastlat )
	    {
	      if ( j == 0 )
		maxlat = genYmax(grid_center_lat[index], grid_center_lat[index+xsize]);
	      else
		maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);

	      if ( j == (ysize-1) )
		minlat = genYmin(grid_center_lat[index], grid_center_lat[index-xsize]);
	      else
		minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
	    }
	  else
	    {
	      if ( j == 0 )
		minlat = genYmin(grid_center_lat[index], grid_center_lat[index+xsize]);
	      else
		minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);

	      if ( j == (ysize-1) )
		maxlat = genYmax(grid_center_lat[index], grid_center_lat[index-xsize]);
	      else
		maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
	    }
	}

      for ( i = 0; i < xsize; ++i )
	{
	  index = (j<<2)*xsize + (i<<2);
	  grid_corner_lat[index  ] = minlat;
	  grid_corner_lat[index+1] = minlat;
	  grid_corner_lat[index+2] = maxlat;
	  grid_corner_lat[index+3] = maxlat;
	}
    }
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
334
335
void gridGenRotBounds(int gridID, int nx, int ny,
		      double *xbounds, double *ybounds, double *xbounds2D, double *ybounds2D)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
  long i, j, index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
338
  double minlon, maxlon;
  double minlat, maxlat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
  double xpole, ypole, angle;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
341
342

  xpole = gridInqXpole(gridID);
  ypole = gridInqYpole(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
  angle = gridInqAngle(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363

  for ( j = 0; j < ny; j++ )
    {
      if ( ybounds[0] > ybounds[1] )
	{
	  maxlat = ybounds[2*j];
	  minlat = ybounds[2*j+1];
	}
      else
	{
	  maxlat = ybounds[2*j+1];
	  minlat = ybounds[2*j];
	}

      for ( i = 0; i < nx; i++ )
	{
	  minlon = xbounds[2*i];
	  maxlon = xbounds[2*i+1];

	  index = j*4*nx + 4*i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
364
365
366
367
368
369
370
371
372
	  xbounds2D[index+0] = lamrot_to_lam(minlat, minlon, ypole, xpole, angle);
	  xbounds2D[index+1] = lamrot_to_lam(minlat, maxlon, ypole, xpole, angle);
	  xbounds2D[index+2] = lamrot_to_lam(maxlat, maxlon, ypole, xpole, angle);
	  xbounds2D[index+3] = lamrot_to_lam(maxlat, minlon, ypole, xpole, angle);

	  ybounds2D[index+0] = phirot_to_phi(minlat, minlon, ypole, angle);
	  ybounds2D[index+1] = phirot_to_phi(minlat, maxlon, ypole, angle);
	  ybounds2D[index+2] = phirot_to_phi(maxlat, maxlon, ypole, angle);
	  ybounds2D[index+3] = phirot_to_phi(maxlat, minlon, ypole, angle);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
374
375
376
	}
    }
}

377
378

void grid_gen_xbounds2D(int nx, int ny, const double *restrict xbounds, double *restrict xbounds2D)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
{
380
  int index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
381
382
  double minlon, maxlon;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
#if defined(_OPENMP)
384
#pragma omp parallel for default(none)        \
385
386
                          shared(nx, ny, xbounds, xbounds2D) \
                         private(minlon, maxlon, index)
387
#endif
388
  for ( int i = 0; i < nx; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389
    {
390
      minlon = xbounds[2*i  ];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
392
      maxlon = xbounds[2*i+1];

393
      for ( int j = 0; j < ny; ++j )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
395
	{
	  index = j*4*nx + 4*i;
396
	  xbounds2D[index  ] = minlon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
398
399
400
401
402
403
	  xbounds2D[index+1] = maxlon;
	  xbounds2D[index+2] = maxlon;
	  xbounds2D[index+3] = minlon;
	}
    }
}

404
405

void grid_gen_ybounds2D(int nx, int ny, const double *restrict ybounds, double *restrict ybounds2D)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
{
407
  int index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408
409
  double minlat, maxlat;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
#if defined(_OPENMP)
411
#pragma omp parallel for default(none)        \
412
413
                          shared(nx, ny, ybounds, ybounds2D) \
                         private(minlat, maxlat, index)
414
#endif
415
  for ( int j = 0; j < ny; ++j )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
417
418
    {
      if ( ybounds[0] > ybounds[1] )
	{
419
	  maxlat = ybounds[2*j  ];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
421
422
423
424
	  minlat = ybounds[2*j+1];
	}
      else
	{
	  maxlat = ybounds[2*j+1];
425
	  minlat = ybounds[2*j  ];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
427
	}

428
      for ( int i = 0; i < nx; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
430
	{
	  index = j*4*nx + 4*i;
431
	  ybounds2D[index  ] = minlat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
432
433
434
435
436
437
438
	  ybounds2D[index+1] = minlat;
	  ybounds2D[index+2] = maxlat;
	  ybounds2D[index+3] = maxlat;
	}
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453
char *gen_param(const char *fmt, ...)
{
  va_list args;
  char str[256];
  char *rstr;
  int len;

  va_start(args, fmt);

  len = vsprintf(str, fmt, args);

  va_end(args);

  len++;
454
  rstr = (char*) Malloc(len*sizeof(char));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
456
  memcpy(rstr, str, len*sizeof(char));

457
  return rstr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
458
459
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
461
462
463
464
465
466
467
static
void lcc_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
  double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
  double zlat, zlon;
  double xi, xj;
  int projflag, scanflag;
  long i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468
  proj_info_t proj;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469

Uwe Schulzweida's avatar
q    
Uwe Schulzweida committed
470
  gridInqLCC(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm, &projflag, &scanflag);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
472
473
474
475
  /*
    while ( originLon < 0 ) originLon += 360;
    while ( lonParY   < 0 ) lonParY   += 360;
  */
  if ( IS_NOT_EQUAL(xincm, yincm) )
476
    Warning("X and Y increment must be equal on Lambert Conformal grid (Xinc = %g, Yinc = %g)\n", 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
	    xincm, yincm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
  if ( IS_NOT_EQUAL(lat1, lat2) )
480
    Warning("Lat1 and Lat2 must be equal on Lambert Conformal grid (Lat1 = %g, Lat2 = %g)\n", 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
	    lat1, lat2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
483
484
  */
  map_set(PROJ_LC, originLat, originLon, xincm, lonParY, lat1, lat2, &proj);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
486
487
488
  for ( i = 0; i < gridsize; i++ )
    {
      xi = xvals[i];
      xj = yvals[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
490
      // status = W3FB12(xi, xj, originLat, originLon, xincm, lonParY, lat1, &zlat, &zlon);
      ijll_lc(xi, xj, proj, &zlat, &zlon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
492
493
494
495
496
497
498
      xvals[i] = zlon;
      yvals[i] = zlat;
    }
}

static
void sinusoidal_to_geo(int gridsize, double *xvals, double *yvals)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
499
#if defined(HAVE_LIBPROJ)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
501
502
503
  char *params[20];
  projUV data, res;
  long i;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
  int nbpar = 0;
505
506
  params[nbpar++] = gen_param("proj=sinu");
  params[nbpar++] = gen_param("ellps=WGS84");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
508
509
510
511

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%ld] = %s", i+1, params[i]);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
512
513
  projPJ proj = pj_init(nbpar, params);
  if ( !proj )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
514
515
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

516
517
  for ( i = 0; i < nbpar; ++i ) Free(params[i]);

518
  /* proj->over = 1; */		/* allow longitude > 180 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
520
521
522
523

  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524
      res = pj_inv(data, proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
526
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
Uwe Schulzweida's avatar
q    
Uwe Schulzweida committed
527
528
      if ( xvals[i] < -9000. || xvals[i] > 9000. ) xvals[i] = -9999.;
      if ( yvals[i] < -9000. || yvals[i] > 9000. ) yvals[i] = -9999.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
531

  pj_free(proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
533
534
535
536
537
538
539
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

static
void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
#if defined(HAVE_LIBPROJ)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
541
542
543
544
545
546
547
  char *params[20];
  projUV data, res;
  double a, lon_0, lat_0;
  long i;

  gridInqLaea(gridID, &a , &lon_0, &lat_0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
  int nbpar = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
  params[nbpar++] = gen_param("proj=laea");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
  if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
552
553
554
555
556
557
  params[nbpar++] = gen_param("lon_0=%g", lon_0);
  params[nbpar++] = gen_param("lat_0=%g", lat_0);

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%d] = %s", i+1, params[i]);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
559
  projPJ proj = pj_init(nbpar, &params[0]);
  if ( !proj )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
561
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

562
  for ( i = 0; i < nbpar; ++i ) Free(params[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563

564
  /* proj->over = 1; */		/* allow longitude > 180 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
565
566
567
568
569

  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
570
      res = pj_inv(data, proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
Uwe Schulzweida's avatar
q    
Uwe Schulzweida committed
573
574
      if ( xvals[i] < -9000. || xvals[i] > 9000. ) xvals[i] = -9999.;
      if ( yvals[i] < -9000. || yvals[i] > 9000. ) yvals[i] = -9999.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
577

  pj_free(proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
578
579
580
581
582
583
584
585
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

static
void lcc2_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
#if defined(HAVE_LIBPROJ)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
587
588
589
590
591
592
593
  char *params[20];
  projUV data, res;
  double a, lon_0, lat_0, lat_1, lat_2;
  long i;

  gridInqLcc2(gridID, &a , &lon_0, &lat_0, &lat_1, &lat_2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
594
  int nbpar = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
595
596
597
598
599
600
601
602
603
604
605
  params[nbpar++] = gen_param("proj=lcc");
  if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
  params[nbpar++] = gen_param("lon_0=%g", lon_0);
  params[nbpar++] = gen_param("lat_0=%g", lat_0);
  params[nbpar++] = gen_param("lat_1=%g", lat_1);
  params[nbpar++] = gen_param("lat_2=%g", lat_2);

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
606
607
  projPJ proj = pj_init(nbpar, &params[0]);
  if ( !proj )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
608
609
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

610
  for ( i = 0; i < nbpar; ++i ) Free(params[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611

612
  /* proj->over = 1; */		/* allow longitude > 180 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
614
615
616
617
  
  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
      res = pj_inv(data, proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
620
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
621
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
623

  pj_free(proj);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
624
625
626
627
628
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
/*
 * grib_get_reduced_row: code from GRIB_API 1.10.4
 *
 * Description:
 *   computes the number of points within the range lon_first->lon_last and the zero based indexes
 *   ilon_first,ilon_last of the first and last point given the number of points along a parallel (pl)
 *
 */
static
void grib_get_reduced_row(long pl,double lon_first,double lon_last,long* npoints,long* ilon_first, long* ilon_last )
{
  double range=0,dlon_first=0,dlon_last=0;
  long irange;

  range=lon_last-lon_first;
  if (range<0) {range+=360;lon_first-=360;}

  /* computing integer number of points and coordinates without using floating point resolution*/
  *npoints=(range*pl)/360.0+1;
  *ilon_first=(lon_first*pl)/360.0;
  *ilon_last=(lon_last*pl)/360.0;

  irange=*ilon_last-*ilon_first+1;

  if (irange != *npoints) {
    if (irange > *npoints) {
      /* checking if the first point is out of range*/
      dlon_first=((*ilon_first)*360.0)/pl;
      if (dlon_first < lon_first) {(*ilon_first)++;irange--;
      }

      /* checking if the last point is out of range*/
      dlon_last=((*ilon_last)*360.0)/pl;
      if (dlon_last > lon_last) {(*ilon_last)--;irange--;
      }
    } else {
      int ok=0;
      /* checking if the point before the first is in the range*/
      dlon_first=((*ilon_first-1)*360.0)/pl;
      if (dlon_first > lon_first) {(*ilon_first)--;irange++;ok=1;
      }

      /* checking if the point after the last is in the range*/
      dlon_last=((*ilon_last+1)*360.0)/pl;
      if (dlon_last < lon_last) {(*ilon_last)++;irange++;ok=1;
      }

      /* if neither of the two are triggered then npoints is too large */
      if (!ok) {(*npoints)--;
      }
    }

    //   assert(*npoints==irange);
  } else {
	  /* checking if the first point is out of range*/
	  dlon_first=((*ilon_first)*360.0)/pl;
	  if (dlon_first < lon_first) {
		  (*ilon_first)++;(*ilon_last)++;
	  }
  }

  if (*ilon_first<0) *ilon_first+=pl;
}
692

693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
static
int qu2reg_subarea(int gridsize, int np, double xfirst, double xlast, 
		   double *array, int *rowlon, int ny, double missval, int *iret, int lmiss, int lperio, int lveggy)
{
  int nx = 0;
  /* sub area (longitudes) */
  long ilon_firstx;
  long ilon_first, ilon_last;
  int i, j;
  long row_count;
  int rlon, nwork = 0;
  int np4 = np*4;
  int size = 0;
  int wlen;
  double *work, **pwork;
  int ii;
  double *parray;

  if ( np <= 0 ) cdoAbort("Number of values between pole and equator missing!");

  grib_get_reduced_row(np4, xfirst, xlast, &row_count, &ilon_firstx, &ilon_last);
  nx = row_count;
  // printf("nx %d  %ld %ld lon1 %g lon2 %g\n", nx, ilon_firstx, ilon_last, (ilon_firstx*360.)/np4, (ilon_last*360.)/np4);

  for ( j = 0; j < ny; ++j ) nwork += rowlon[j];

719
720
  pwork = (double **) Malloc(ny*sizeof(double *));
  work  = (double*) Malloc(ny*np4*sizeof(double));
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
  wlen = 0;
  pwork[0] = work;
  for ( j = 1; j < ny; ++j )
    {
      wlen += rowlon[j-1];
      pwork[j] = work + wlen;
    } 
  // printf(" ny, np4, nwork %d %d %d wlen %d\n", ny, np4, nwork, wlen);

  for ( j = 0; j < ny; ++j )
    {
      rlon = rowlon[j];
      for ( i = 0; i < rlon; ++i ) pwork[j][i] = missval;
    } 

  parray = array;
  for ( j = 0; j < ny; ++j )
    {
      rlon = rowlon[j];
      row_count = 0;
      grib_get_reduced_row(rlon, xfirst, xlast, &row_count, &ilon_first, &ilon_last);
      // printf("j %d xfirst %g xlast %g rowlon %d %ld %ld %ld %g %g\n", j, xfirst, xlast, rlon, row_count, ilon_first, ilon_last, (ilon_first*360.)/rlon, (ilon_last*360.)/rlon);

      for ( i = ilon_first; i < (ilon_first+row_count); ++i )
	{
	  ii = i;
	  if ( ii >= rlon ) ii -= rlon; 
	  pwork[j][ii] = *parray;
	  parray++;
	}
      size += row_count;
    }

  if ( gridsize != size ) cdoAbort("gridsize1 inconsistent!");

756
  (void) qu2reg3_double(work, rowlon, ny, np4, missval, iret, lmiss, lperio, lveggy);
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778

  wlen = 0;
  pwork[0] = work;
  for ( j = 1; j < ny; ++j )
    {
      wlen += np4;
      pwork[j] = work + wlen;
    } 

  // printf("nx, ilon_firstx %d %ld\n", nx, ilon_firstx);
  parray = array;
  for ( j = 0; j < ny; ++j )
    {
      for ( i = ilon_firstx; i < (ilon_firstx+nx); ++i )
	{
	  ii = i;
	  if ( ii >= np4 ) ii -= np4; 
	  *parray = pwork[j][ii];
	  parray++;
	}
    }

779
780
  Free(work);
  Free(pwork);
781

782
  return nx;
783
784
785
}


786
787
void field2regular(int gridID1, int gridID2, double missval, double *array, int nmiss)
{
788
  int nx = 0, ny, np;
789
790
791
  int gridtype;
  int lmiss, lperio, lveggy;
  int iret;
792
793
794
  int *rowlon;
  double xfirstandlast[2];
  double xfirst, xlast;
795
796
797

  gridtype = gridInqType(gridID1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
798
  if ( gridtype != GRID_GAUSSIAN_REDUCED ) Error("Not a reduced Gaussian grid!");
799

800
801
802
  lmiss = nmiss > 0;
  lperio = 1;
  lveggy = 0;
803

804
  ny = gridInqYsize(gridID1);
805
  np = gridInqNP(gridID1);
806

807
  rowlon = (int*) Malloc(ny*sizeof(int));
808
  gridInqRowlon(gridID1, rowlon);
809

810
811
812
813
814
  xfirstandlast[0] = 0.;
  xfirstandlast[1] = 0.;
  gridInqXvals(gridID1, xfirstandlast);
  xfirst = xfirstandlast[0];
  xlast  = xfirstandlast[1];
815

816
  if ( fabs(xfirst) > 0 || (np > 0 && fabs(xlast - (360.0-90.0/np)) > 90.0/np) )
817
    {
818
      nx = qu2reg_subarea(gridInqSize(gridID1), np, xfirst, xlast, array, rowlon, ny, missval, &iret, lmiss, lperio, lveggy);
819
820
821
    }
  else
    {
822
      nx = 2*ny;
823
      (void) qu2reg3_double(array, rowlon, ny, nx, missval, &iret, lmiss, lperio, lveggy);
824
    }
825

826
  if ( gridInqSize(gridID2) != nx*ny ) Error("Gridsize differ!");
827

828
  Free(rowlon);
829
830
831
832
833
834
835
}


int gridToRegular(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
836
  int nx = 0, ny, np;
837
838
  long i;
  double *xvals = NULL, *yvals = NULL;
839
840
  double xfirstandlast[2];
  double xfirst, xlast;
841
842
843

  gridtype = gridInqType(gridID1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
844
  if ( gridtype != GRID_GAUSSIAN_REDUCED ) Error("Not a reduced Gaussian grid!");
845
846

  ny = gridInqYsize(gridID1);
847
848
  np = gridInqNP(gridID1);

849
  yvals = (double*) Malloc(ny*sizeof(double));
850
851
852
853
854
855
856
857
  gridInqYvals(gridID1, yvals);

  xfirstandlast[0] = 0.;
  xfirstandlast[1] = 0.;
  gridInqXvals(gridID1, xfirstandlast);
  xfirst = xfirstandlast[0];
  xlast  = xfirstandlast[1];

858
  if ( fabs(xfirst) > 0 || (np > 0 && fabs(xlast - (360.0-90.0/np)) > 90.0/np) )
859
860
    {
      /* sub area (longitudes) */
861
      long ilon_first, ilon_last;
862
863
864
865
866
867
868
869
870
      long row_count;
      int np4 = np*4;
      int *rowlon = NULL;

      if ( np <= 0 ) cdoAbort("Number of values between pole and equator missing!");

      grib_get_reduced_row(np4, xfirst, xlast, &row_count, &ilon_first, &ilon_last);

      nx = row_count;
871
      xvals = (double*) Malloc(nx*sizeof(double));
872
873
874
875
876
877
      for ( i = 0; i < nx; ++i )
	{
	  xvals[i] = ((ilon_first+i)*360.)/np4;
	  if ( xfirst > xlast ) xvals[i] -= 360.;
	}

878
      Free(rowlon);
879
880
881
882
    }
  else
    {
      nx = 2*ny;
883
      xvals = (double*) Malloc(nx*sizeof(double));
884
885
886
      for ( i = 0; i < nx; ++i ) xvals[i] = i * 360./nx;
    }

887
888
889
890
891
892
893
894
895
  gridsize = nx*ny;

  gridID2  = gridCreate(GRID_GAUSSIAN, gridsize);
	  
  gridDefXsize(gridID2, nx);
  gridDefYsize(gridID2, ny);
  
  gridDefXvals(gridID2, xvals);
  gridDefYvals(gridID2, yvals);
896
  gridDefNP(gridID2, np);
897

898
899
  Free(xvals);
  Free(yvals);
900

901
  return gridID2;
902
903
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
904
905
906
907
908
static
void gridCopyMask(int gridID1, int gridID2, long gridsize)
{
  if ( gridInqMask(gridID1, NULL) )
    {
909
      int *mask = (int*) Malloc(gridsize*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
910
911
      gridInqMask(gridID1, mask);
      gridDefMask(gridID2, mask);
912
      Free(mask);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
913
914
915
    }
}

916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
static
int check_range(long n, double *vals, double valid_min, double valid_max)
{
  int status = 0;
  long i;

  for ( i = 0; i < n; ++i )
    {
      if ( vals[i] < valid_min || vals[i] > valid_max )
	{
	  status = 1;
	  break;
	}
    }

931
  return status;
932
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
933

934

935
int gridToCurvilinear(int gridID1, int lbounds)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
936
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
937
  long index;
938
939
940
  int gridtype = gridInqType(gridID1);
  size_t gridsize = (size_t) gridInqSize(gridID1);
  int gridID2  = gridCreate(GRID_CURVILINEAR, (int) gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
941
  gridDefPrec(gridID2, DATATYPE_FLT32);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
943
944
945
946
947
	  
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
    case GRID_LCC:
948
    case GRID_LCC2:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
949
    case GRID_LAEA:
950
    case GRID_SINUSOIDAL:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
951
      {
952
	double xscale = 1, yscale = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953
954
	double *xvals = NULL, *yvals = NULL;
	double *xbounds = NULL, *ybounds = NULL;
955
	char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956

957
958
	int nx = gridInqXsize(gridID1);
	int ny = gridInqYsize(gridID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
959

960
961
962
	gridInqXunits(gridID1, xunits);
	gridInqYunits(gridID1, yunits);

963
964
965
966
967
968
969
	if ( gridtype == GRID_LAEA )
	  {
	    int lvalid_xunits = FALSE;
	    int lvalid_yunits = FALSE;
	    int len;
	    len = (int) strlen(xunits);
	    if ( len == 1 && memcmp(xunits, "m",  1) == 0 ) lvalid_xunits = TRUE;
970
	    if ( len == 2 && memcmp(xunits, "km", 2) == 0 ) lvalid_xunits = TRUE;
971
972
	    len = (int) strlen(yunits);
	    if ( len == 1 && memcmp(yunits, "m",  1) == 0 ) lvalid_yunits = TRUE;
973
	    if ( len == 2 && memcmp(yunits, "km", 2) == 0 ) lvalid_yunits = TRUE;
974
975

	    if ( lvalid_xunits == FALSE )
976
	      cdoWarning("Possibly wrong result! Invalid x-coordinate units: \"%s\" (expected \"m\" or \"km\")", xunits);
977
	    if ( lvalid_yunits == FALSE )
978
	      cdoWarning("Possibly wrong result! Invalid y-coordinate units: \"%s\" (expected \"m\" or \"km\")", yunits);
979
980
	  }

981
982
	if ( memcmp(xunits, "km", 2) == 0 ) xscale = 1000;
	if ( memcmp(yunits, "km", 2) == 0 ) yscale = 1000;
983

Uwe Schulzweida's avatar
Uwe Schulzweida committed
984
985
986
	gridDefXsize(gridID2, nx);
	gridDefYsize(gridID2, ny);

987
988
	double *xvals2D = (double*) Malloc(gridsize*sizeof(double));
	double *yvals2D = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
989
990
991
992


	if ( gridtype == GRID_LCC )
	  {
993
994
	    for ( int j = 0; j < ny; j++ )
	      for ( int i = 0; i < nx; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
996
997
		  xvals2D[j*nx+i] = i+1;
		  yvals2D[j*nx+i] = j+1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
998
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
999
1000

	    lcc_to_geo(gridID1, gridsize, xvals2D, yvals2D);
For faster browsing, not all history is shown. View entire blame