YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
read_fesom_grid.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#ifdef HAVE_CONFIG_H
6#include "config.h"
7#endif
8
9#include <stdlib.h>
10#include <stdio.h>
11
12#include "read_fesom_grid.h"
13#include "geometry.h"
14#include "utils_common.h"
15#include "io_utils.h"
16
17#ifdef YAC_NETCDF_ENABLED
18
19#include <netcdf.h>
20
21static void get_fesom_vertices ( int ncid,
22 double **vertex_lon,
23 double **vertex_lat,
24 int * nbr_cells,
25 int * nbr_vertices_per_cell);
26
27static void get_fesom_cell_center ( int ncid,
28 double **cell_lon,
29 double **cell_lat,
30 size_t *nbr_cells );
31
32static void remove_duplicated_vertices(double * temp_vertex_lon,
33 double * temp_vertex_lat,
34 int * temp_nbr_vertices,
35 int * old_to_new_id);
36
37void yac_read_fesom_grid_information(const char * filename, int * nbr_vertices,
38 int * nbr_cells, int ** num_vertices_per_cell,
39 int ** cell_to_vertex, double ** x_vertices,
40 double ** y_vertices, double ** x_cells,
41 double ** y_cells) {
42
43 int ncid;
44
45 /* Open file */
46
47 yac_nc_open(filename, NC_NOWRITE, &ncid);
48
49 /* Get vertex longitudes and latitudes of cells and
50 relations between vertices and cells*/
51
52 {
53 double * temp_vertex_lon;
54 double * temp_vertex_lat;
55 int nbr_vertices_per_cell;
56 int temp_nbr_vertices;
57
58 get_fesom_vertices (ncid, &temp_vertex_lon, &temp_vertex_lat,
59 nbr_cells, &nbr_vertices_per_cell);
60
61 temp_nbr_vertices = *nbr_cells * nbr_vertices_per_cell;
62
63 int * old_to_new_id = xmalloc(temp_nbr_vertices * sizeof(*old_to_new_id));
64
65 remove_duplicated_vertices(temp_vertex_lon, temp_vertex_lat,
66 &temp_nbr_vertices, old_to_new_id);
67
68 *nbr_vertices = temp_nbr_vertices;
69
70 *x_vertices = xrealloc(temp_vertex_lon,
71 temp_nbr_vertices * sizeof(**x_vertices));
72 *y_vertices = xrealloc(temp_vertex_lat,
73 temp_nbr_vertices * sizeof(**y_vertices));
74
75 *num_vertices_per_cell = xmalloc(*nbr_cells * sizeof(**num_vertices_per_cell));
76 for (int i = 0; i < *nbr_cells; ++i)
77 (*num_vertices_per_cell)[i] = nbr_vertices_per_cell;
78
79 *cell_to_vertex = xmalloc(*nbr_cells * nbr_vertices_per_cell *
80 sizeof(**cell_to_vertex));
81
82 // Unfortunately the data is only available in Fortran order
83 for (int i = 0, k = 0; i < *nbr_cells; ++i)
84 for (int j = 0; j < nbr_vertices_per_cell; ++j, ++k)
85 (*cell_to_vertex)[k] = old_to_new_id[k];
86
87 int new_size_cell_to_vertex = 0;
88
89 // remove duplicated corners in each cell
90 for (int i = 0; i < *nbr_cells; ++i) {
91 (*cell_to_vertex)[new_size_cell_to_vertex++] =
92 (*cell_to_vertex)[i * nbr_vertices_per_cell];
93 for (int j = 1; j < nbr_vertices_per_cell; ++j)
94 if ((*cell_to_vertex)[i * nbr_vertices_per_cell + j] !=
95 (*cell_to_vertex)[i * nbr_vertices_per_cell + j - 1])
96 (*cell_to_vertex)[new_size_cell_to_vertex++] =
97 (*cell_to_vertex)[i * nbr_vertices_per_cell + j];
98 else
99 (*num_vertices_per_cell)[i]--;
100 }
101
102 *cell_to_vertex = xrealloc(*cell_to_vertex, new_size_cell_to_vertex *
103 sizeof(**cell_to_vertex));
104
105 free(old_to_new_id);
106
107 }
108
109 {
110 /* Get longitudes and latitudes of cell center points */
111
112 size_t temp_nbr_cells;
113
114 get_fesom_cell_center (ncid, x_cells, y_cells, &temp_nbr_cells);
115 }
116 YAC_HANDLE_ERROR(nc_close(ncid));
117}
118
120 char const * filename) {
121
122 int nbr_vertices;
123 int nbr_cells;
124 int * num_vertices_per_cell = NULL;
125 int * cell_to_vertex = NULL;
126
127 double * x_vertices = NULL;
128 double * y_vertices = NULL;
129 double * x_cells = NULL;
130 double * y_cells = NULL;
131
132 yac_read_fesom_grid_information(filename, &nbr_vertices, &nbr_cells,
134 &x_vertices, &y_vertices,
135 &x_cells, &y_cells);
136
137 free(x_cells);
138 free(y_cells);
139
140 struct yac_basic_grid_data grid_data =
142 (size_t)nbr_vertices, (size_t)nbr_cells, num_vertices_per_cell,
143 x_vertices, y_vertices, cell_to_vertex);
145 free(x_vertices);
146 free(y_vertices);
147 free(cell_to_vertex);
148
149 return grid_data;
150}
151
153 char const * filename, char const * gridname) {
154
155 return
157 gridname,
159}
160
161/* ---------------------------------------------------------------- */
162
163static void get_fesom_vertices ( int ncid,
164 double **vertex_lon,
165 double **vertex_lat,
166 int * nbr_cells,
167 int * nbr_vertices_per_cell) {
168
169 int glat_id; // Various NetCDF IDs
170 int glon_id;
171
172 int glon_ndims; // number of dimensions in file
173 int glat_ndims;
174
175 int glon_dimids[NC_MAX_VAR_DIMS]; // dimension NetCDF IDs
176 int glat_dimids[NC_MAX_VAR_DIMS];
177
178 nc_type glon_type; // variable type
179 nc_type glat_type;
180
181 int nbr_atts; // number of attributes
182
183 size_t nbr_vertices;
184
185 yac_nc_inq_varid (ncid, "lon_vertices", &glon_id);
186 yac_nc_inq_varid (ncid, "lat_vertices", &glat_id);
187
189 nc_inq_var (ncid, glon_id, 0, &glon_type, &glon_ndims, glon_dimids, &nbr_atts));
190
192 nc_inq_var (ncid, glat_id, 0, &glat_type, &glat_ndims, glat_dimids, &nbr_atts));
193
194 /* Allocate memory to read in coordinates */
196 glon_type == glat_type,
197 "ERROR(get_fesom_vertices): lon and lat datatypes do not match");
199 glon_ndims == 2 && glat_ndims == 2,
200 "ERROR(get_fesom_vertices): unsupported number of dimensions for lon or lat");
201
202 /* get size of arrays */
203
204 size_t dimlen_lat[2]; // dimension size for latitude
205 size_t dimlen_lon[2]; // dimension size for longitude
206
207 for (int i = 0; i < 2; ++i) {
208 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, glon_dimids[i], dimlen_lon + i));
209 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, glat_dimids[i], dimlen_lat + i));
211 dimlen_lon[i] == dimlen_lat[i],
212 "ERROR(get_fesom_vertices): mismatching dimension size for lon and lat");
213 }
214
215 nbr_vertices = dimlen_lon[0] * dimlen_lon[1];
216 *nbr_cells = dimlen_lon[0];
217 *nbr_vertices_per_cell = dimlen_lon[1];
218
219 /* read coordinates and convert radians into degrees */
220
221
222 *vertex_lon = (double * ) xmalloc ( nbr_vertices * sizeof ( double ) );
223 *vertex_lat = (double * ) xmalloc ( nbr_vertices * sizeof ( double ) );
224
225 YAC_HANDLE_ERROR(nc_get_var_double (ncid, glon_id, (*vertex_lon)));
226 YAC_HANDLE_ERROR(nc_get_var_double (ncid, glat_id, (*vertex_lat)));
227
228 for (size_t i = 0; i < nbr_vertices; ++i) {
229 (*vertex_lon)[i] *= YAC_RAD;
230 (*vertex_lat)[i] *= YAC_RAD;
231 }
232}
233
234/* ---------------------------------------------------------------- */
235
236static void get_fesom_cell_center ( int ncid,
237 double **cell_lon,
238 double **cell_lat,
239 size_t *nbr_cells ) {
240
241 int glat_id; // Various NetCDF IDs
242 int glon_id;
243
244 int glon_ndims; // number of dimensions in file
245 int glat_ndims;
246
247 int glon_dimids[NC_MAX_VAR_DIMS]; // dimension NetCDF IDs
248 int glat_dimids[NC_MAX_VAR_DIMS];
249
250 nc_type glon_type; // variable type
251 nc_type glat_type;
252
253 int nbr_atts; // number of attributes
254
255 yac_nc_inq_varid (ncid, "lon", &glon_id);
256 yac_nc_inq_varid (ncid, "lat", &glat_id);
257
259 nc_inq_var (ncid, glon_id, 0, &glon_type, &glon_ndims, glon_dimids, &nbr_atts));
261 nc_inq_var (ncid, glat_id, 0, &glat_type, &glat_ndims, glat_dimids, &nbr_atts));
262
263 /* Allocate memory to read in coordinates */
265 glon_type == glat_type,
266 "ERROR(get_fesom_cell_center): lon and lat datatypes do not match");
268 glon_ndims == 1 && glat_ndims == 1,
269 "ERROR(get_fesom_cell_center): unsupported number of dimensions for lon or lat");
270
271 /* get size of arrays */
272
273 size_t dimlen_lat; // dimension size for latitude
274 size_t dimlen_lon; // dimension size for longitude
275
276 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, glon_dimids[0], &dimlen_lon));
277 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, glat_dimids[0], &dimlen_lat));
279 dimlen_lon == dimlen_lat,
280 "ERROR(get_fesom_cell_center): mismatching dimension size for lon and lat");
281
282 *nbr_cells = dimlen_lon;
283
284 /* read coordinates and convert radians into degrees */
285
286 *cell_lon = (double * ) xmalloc ( *nbr_cells * sizeof ( **cell_lon ) );
287 *cell_lat = (double * ) xmalloc ( *nbr_cells * sizeof ( **cell_lat ) );
288
289 YAC_HANDLE_ERROR(nc_get_var_double (ncid, glon_id, *cell_lon));
290 YAC_HANDLE_ERROR(nc_get_var_double (ncid, glat_id, *cell_lat));
291
292 for (size_t i = 0; i < *nbr_cells; ++i) {
293 (*cell_lon)[i] *= YAC_RAD;
294 (*cell_lat)[i] *= YAC_RAD;
295 }
296}
297
298/* ---------------------------------------------------------------- */
299
301
302 struct {
303 double lon, lat;
304 } p;
305 int i;
306};
307
308static int compare_point_with_index(const void * a,const void * b) {
309
310 const struct point_with_index * a_ = (struct point_with_index*)a;
311 const struct point_with_index * b_ = (struct point_with_index*)b;
312
313 int lon_diff = fabs(a_->p.lon - b_->p.lon) > 1e-9;
314 int lat_diff = fabs(a_->p.lat - b_->p.lat) > 1e-9;
315
316 if (lon_diff) {
317
318 if (a_->p.lon > b_->p.lon) return -1;
319 else return 1;
320
321 } else if (lat_diff) {
322
323 if (a_->p.lat > b_->p.lat) return -1;
324 else return 1;
325 } else
326 return 0;
327}
328
329static void remove_duplicated_vertices(double * temp_vertex_lon,
330 double * temp_vertex_lat,
331 int * temp_nbr_vertices,
332 int * old_to_new_id) {
333
334 struct point_with_index * sort_array =
335 xmalloc(*temp_nbr_vertices * sizeof(*sort_array));
336
337 for (int i = 0; i < *temp_nbr_vertices; ++i) {
338
339 double curr_lon, curr_lat;
340
341 curr_lon = ((double*)temp_vertex_lon)[i];
342 curr_lat = ((double*)temp_vertex_lat)[i];
343
344 while (curr_lon < 0.0) curr_lon += 360.0;
345 while (curr_lon >= 360) curr_lon -= 360.0;
346
347 sort_array[i].p.lon = curr_lon;
348 sort_array[i].p.lat = curr_lat;
349
350 sort_array[i].i = i;
351 }
352
353 yac_mergesort(sort_array, *temp_nbr_vertices, sizeof(*sort_array),
355
356 old_to_new_id[sort_array[0].i] = 1;
357
358 int last_unique_idx = sort_array[0].i;
359
360 for (int i = 1; i < *temp_nbr_vertices; ++i) {
361
362 if (compare_point_with_index(sort_array + i - 1, sort_array + i)) {
363
364 old_to_new_id[sort_array[i].i] = 1;
365 last_unique_idx = sort_array[i].i;
366
367 } else {
368
369 old_to_new_id[sort_array[i].i] = -last_unique_idx;
370 }
371 }
372
373 free(sort_array);
374
375 size_t new_nbr_vertices = 0;
376
377 for (int i = 0; i < *temp_nbr_vertices; ++i) {
378
379 if (old_to_new_id[i] == 1) {
380
381 temp_vertex_lon[new_nbr_vertices] = temp_vertex_lon[i];
382 temp_vertex_lat[new_nbr_vertices] = temp_vertex_lat[i];
383
384 old_to_new_id[i] = new_nbr_vertices;
385
386 new_nbr_vertices++;
387 }
388 }
389
390 for (int i = 0; i < *temp_nbr_vertices; ++i)
391 if (old_to_new_id[i] <= 0)
392 old_to_new_id[i] = old_to_new_id[-old_to_new_id[i]];
393
394 *temp_nbr_vertices = new_nbr_vertices;
395}
396
397#else
398
400 char const * filename) {
401
402 UNUSED(filename);
403 die(
404 "ERROR(yac_read_fesom_basic_grid_data): "
405 "YAC is built without the NetCDF support");
406
407 return
409 (size_t[]){0,0}, (int[]){0,0}, NULL, NULL);
410}
411
413 char const * filename, char const * gridname) {
414
415 UNUSED(filename);
416 UNUSED(gridname);
417 die(
418 "ERROR(yac_read_fesom_basic_grid): "
419 "YAC is built without the NetCDF support");
420
421 return NULL;
422}
423
424void yac_read_fesom_grid_information(const char * filename, int * nbr_vertices,
425 int * nbr_cells, int ** num_vertices_per_cell,
426 int ** cell_to_vertex, double ** x_vertices,
427 double ** y_vertices, double ** x_cells,
428 double ** y_cells) {
429
430 UNUSED(filename);
431 UNUSED(nbr_vertices);
432 UNUSED(nbr_cells);
433 UNUSED(num_vertices_per_cell);
434 UNUSED(cell_to_vertex);
435 UNUSED(x_vertices);
436 UNUSED(y_vertices);
437 UNUSED(x_cells);
438 UNUSED(y_cells);
439 die(
440 "ERROR(yac_read_fesom_grid_information): "
441 "YAC is built without the NetCDF support");
442}
443
444
445
446#endif // YAC_NETCDF_ENABLED
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:45
struct yac_basic_grid_data yac_generate_basic_grid_data_reg_2d(size_t nbr_vertices[2], int cyclic[2], double *lon_vertices, double *lat_vertices)
Definition grid_reg2d.c:65
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
#define UNUSED(x)
Definition core.h:71
#define YAC_RAD
Definition geometry.h:30
void yac_nc_inq_varid(int ncid, char const *name, int *varidp)
Definition io_utils.c:405
void yac_nc_open(const char *path, int omode, int *ncidp)
Definition io_utils.c:344
#define YAC_HANDLE_ERROR(exp)
Definition io_utils.h:30
void yac_mergesort(void *base, size_t num, size_t size, int(*compar)(const void *, const void *))
Definition mergesort.c:134
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xmalloc(size)
Definition ppm_xfuncs.h:66
static int compare_point_with_index(const void *a, const void *b)
static void remove_duplicated_vertices(double *temp_vertex_lon, double *temp_vertex_lat, int *temp_nbr_vertices, int *old_to_new_id)
void yac_read_fesom_grid_information(const char *filename, int *nbr_vertices, int *nbr_cells, int **num_vertices_per_cell, int **cell_to_vertex, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells)
static void get_fesom_cell_center(int ncid, double **cell_lon, double **cell_lat, size_t *nbr_cells)
struct yac_basic_grid * yac_read_fesom_basic_grid(char const *filename, char const *gridname)
static void get_fesom_vertices(int ncid, double **vertex_lon, double **vertex_lat, int *nbr_cells, int *nbr_vertices_per_cell)
struct yac_basic_grid_data yac_read_fesom_basic_grid_data(char const *filename)
struct point_with_index::@64 p
#define die(msg)
Definition yac_assert.h:12
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16