commit 9334b1ae918dadaef0d73bd137bac2265bb0a495
parent 543e2b5a84cf11da087ba47f92c9945cfa1edee0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 26 Apr 2018 18:37:09 +0200
Make test_nc a program that convers LES into a specific fileformat
Diffstat:
| M | src/test_nc.c | | | 628 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------- |
1 file changed, 557 insertions(+), 71 deletions(-)
diff --git a/src/test_nc.c b/src/test_nc.c
@@ -13,14 +13,147 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 2 /* getopt support */
+
+#include <rsys/rsys.h>
+#include <rsys/math.h>
#include <rsys/mem_allocator.h>
+
+#include <alloca.h>
#include <netcdf.h>
+#include <string.h>
+#include <unistd.h>
+
+/*******************************************************************************
+ * Command arguments
+ ******************************************************************************/
+struct args {
+ const char* output;
+ const char* input;
+ int force_overwrite;
+ int check;
+ int no_output;
+ int quit; /* Quit the application */
+};
+#define ARGS_DEFAULT__ {NULL,NULL,0,0,0,0}
+static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
-#define NCID_INVALID -1
+static void
+print_help(const char* cmd)
+{
+ ASSERT(cmd);
-/*
- * Simple program used to play with the NetCDF API and the LES fileformat
- */
+ printf(
+"Usage: %s -i INPUT [OPTIONS]\n"
+"Convert the LES data stored into INPUT from NetCDF to the htLES fileformat.\n\n",
+ cmd);
+ printf(
+" -c advanced check of the validity of the submitted LES file\n"
+" with respect to the converter pre-requisites. Note that\n"
+" this option can increase significantly the conversion\n"
+" time.\n");
+ printf(
+" -f overwrite the OUTPUT file if it already exists.\n");
+ printf(
+" -h display this help and exit.\n");
+ printf(
+" -i INPUT path of the LES file to convert.\n");
+ printf(
+" -o OUTPUT write results to OUTPUT. If not defined, write results to\n"
+" standard output.\n");
+ printf(
+" -q do not write results to OUTPUT.\n");
+ printf("\n");
+ printf(
+"%s (C) 2018 |Meso|Star>. This is free software released under the GNU GPL\n"
+"license, version 3 or later. You are free to change or redistribute it under\n"
+"certain conditions <http://gnu.org/licenses/gpl.html>.\n", cmd);
+}
+
+static void
+args_release(struct args* args)
+{
+ ASSERT(args);
+ *args = ARGS_DEFAULT;
+}
+
+static res_T
+args_init(struct args* args, const int argc, char** argv)
+{
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(args && argc && argv);
+
+ while((opt = getopt(argc, argv, "chi:o:q")) != -1) {
+ switch(opt) {
+ case 'c': args->check = 1; break;
+ case 'f': args->force_overwrite = 1; break;
+ case 'h':
+ print_help(argv[0]);
+ args_release(args);
+ args->quit = 1;
+ goto exit;
+ case 'i': args->input = optarg; break;
+ case 'o': args->output = optarg; break;
+ case 'q': args->no_output = 1; break;
+ default: RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argumet '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ if(!args->input) {
+ fprintf(stderr, "Missing input file.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ args_release(args);
+ goto exit;
+}
+
+/*******************************************************************************
+ * Grid header
+ ******************************************************************************/
+struct grid {
+ int32_t nx, ny, nz, ntimes; /* Definition */
+ int8_t is_z_irregular;
+ double lower[3]; /* Lower position of the grid */
+ double vxsz_x; /* Size of the voxel is X */
+ double vxsz_y; /* Size of the voxel in Y */
+ double* vxsz_z; /* Size of the voxel in Z */
+};
+
+#define GRID_NULL__ {0,0,0,0,0,{0,0,0},0,0,NULL}
+static const struct grid GRID_NULL = GRID_NULL__;
+
+static void
+grid_release(struct grid* grid)
+{
+ CHK(grid);
+ if(grid->vxsz_z) free(grid->vxsz_z);
+}
+
+/*******************************************************************************
+ * NetCDF helper functions and macros
+ ******************************************************************************/
+#define INVALID_ID -1
+
+#define NC(Func) { \
+ const int err__ = nc_ ## Func; \
+ if(err__ != NC_NOERR) { \
+ fprintf(stderr, "error:%i:%s\n", __LINE__, ncerr_to_str(err__)); \
+ abort(); \
+ } \
+ } (void)0
static INLINE const char*
ncerr_to_str(const int ncerr)
@@ -85,88 +218,441 @@ sizeof_nctype(const nc_type type)
case NC_INT64:
case NC_UINT64:
sz = 8; break;
- default: FATAL("Unreachable code.\n"); break;
+ default: FATAL("Unreachable cde.\n"); break;
}
return sz;
}
-#define NC(Func) { \
- const int err__ = nc_ ## Func; \
- if(err__ != NC_NOERR) { \
- fprintf(stderr, "error:%s:%i:%s:%s\n", \
- __FILE__, __LINE__, FUNC_NAME, ncerr_to_str(err__)); \
- abort(); \
- } \
- } (void)0
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_definition
+ (int nc, int32_t* nx, int32_t* ny, int32_t* nz, int32_t* ntimes)
+{
+ size_t len[4];
+ int dimids[4];
+ int id;
+ int ndims;
+ int err = NC_NOERR;
+ res_T res = RES_OK;
+ ASSERT(nx && ny && nz && ntimes);
-int
-main(int argc, char** argv)
+ err = nc_inq_varid(nc, "RCT", &id);
+ if(err != NC_NOERR) {
+ fprintf(stderr, "Could not inquire the 'RCT' variable -- %s\n",
+ ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_varndims(nc, id, &ndims));
+ if(ndims != 4) {
+ fprintf(stderr, "The dimension of the 'RCT' variable must be 4.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_vardimid(nc, id, dimids));
+ NC(inq_dimlen(nc, dimids[0], len+0));
+ NC(inq_dimlen(nc, dimids[1], len+1));
+ NC(inq_dimlen(nc, dimids[2], len+2));
+ NC(inq_dimlen(nc, dimids[3], len+3));
+
+ *nx = (int32_t)len[3];
+ *ny = (int32_t)len[2];
+ *nz = (int32_t)len[1];
+ *ntimes = (int32_t)len[0];
+
+ if(*nx < 1 || *ny < 1 || *nz < 1 || *ntimes < 1) {
+ fprintf(stderr,
+ "The spatial and time definitions cannot be null.\n"
+ " #x = %i; #y = %i; #z = %i; #time = %i\n",
+ *nx, *ny, *nz, *ntimes);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_regular_dimension
+ (const int nc,
+ const char* var_name,
+ double* voxel_size,
+ double* lower_pos,
+ const int check)
{
- int ncid = NCID_INVALID;
- int id = NCID_INVALID;
+ char* mem = NULL;
+ size_t len, start, i;
+ double vxsz;
nc_type type;
- int ndims;
- int dimid[BUFSIZ];
- size_t startid[BUFSIZ];
- int natts;
- size_t len[BUFSIZ];
- char str[BUFSIZ];
- int err;
- void* mem = NULL;
-
- const char* att_name[] = {
- "W_E_direction",
- "S_N_direction",
- "VLEV",
- "RVT",
- "RCT"
- };
- const size_t natt_names = sizeof(att_name)/sizeof(att_name[0]);
- size_t iatt;
- size_t idim;
- size_t ival;
-
- if(argc < 2) {
- fprintf(stderr, "Usage: %s FILENAME\n", argv[0]);
- goto error;
- }
-
- NC(open(argv[1], NC_WRITE|NC_MMAP, &ncid));
-
- FOR_EACH(iatt, 0, natt_names) {
- NC(inq_varid(ncid, att_name[iatt], &id));
- NC(inq_var(ncid, id, str, &type, &ndims, dimid, &natts));
- printf("%s: type: %s; #dims: %d; #atts: %d\n",
- str, nctype_to_str(type), ndims, natts);
-
- FOR_EACH(idim, 0, ndims) {
- NC(inq_dim(ncid, dimid[idim], str, len));
- printf("\t%s: len: %lu\n", str, (unsigned long)len[0]);
+ int varid, dimid, ndims;
+ int err = NC_NOERR;
+ res_T res = RES_OK;
+ ASSERT(nc && var_name && voxel_size && lower_pos);
+
+ err = nc_inq_varid(nc, var_name, &varid);
+ if(err != NC_NOERR) {
+ fprintf(stderr, "Could not inquire the '%s' variable -- %s\n",
+ var_name, ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ NC(inq_varndims(nc, varid, &ndims));
+ if(ndims != 1) {
+ fprintf(stderr, "The dimension of the '%s' variable must be 1 .\n", var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_vardimid(nc, varid, &dimid));
+ NC(inq_dimlen(nc, dimid, &len));
+ ASSERT(len > 0);
+
+ NC(inq_vartype(nc, varid, &type));
+ if(type != NC_DOUBLE && type != NC_FLOAT) {
+ fprintf(stderr, "The type of the '%s' variable must be float or double.\n",
+ var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define AT(Mem, Id) (type == NC_FLOAT ? ((float*)Mem)[Id] : ((double*)Mem)[Id])
+ if(!check) {
+ char* raw = alloca(8);
+ start = 0;
+ len = 1;
+ NC(get_vara(nc, varid, &start, &len, raw));
+ vxsz = AT(raw, 0) * 2; /* Assume that the grid starts at 0 */
+ } else {
+ mem = mem_alloc(len*sizeof_nctype(type));
+ if(!mem) {
+ fprintf(stderr, "Could not allocate memory for the '%s' variable.\n", var_name);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ start = 0;
+ NC(get_vara(nc, varid, &start, &len, mem));
+ vxsz = AT(mem, 0) * 2; /* Assume that the grid starts from 0 */
+ }
+
+ if(vxsz <= 0) {
+ fprintf(stderr, "The '%s' variable can't have negative or null values.\n",
+ var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(check) {
+ /* Check that the dimension is regular */
+ FOR_EACH(i, 1, len) {
+ if(!eq_eps(AT(mem, i) - AT(mem, i-1), vxsz, 1.e-6)) {
+ fprintf(stderr, "The voxel size of the '%s' variable must be regular.\n",
+ var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ }
+ #undef AT
+
+ *voxel_size = vxsz;
+ *lower_pos = 0;
+exit:
+ if(mem) mem_rm(mem);
+ return res;
+error:
+ goto exit;
+}
+
+/* Return 1 if Z is irregular */
+static int
+setup_Z_dimension
+ (const int nc,
+ double** voxel_size,
+ double* lower_pos,
+ int8_t* is_z_irregular,
+ const int check)
+{
+ enum { Z, X, Y, NDIMS };
+ char* mem = NULL;
+ char* mem2 = NULL;
+ double* pvxsz = NULL;
+ size_t len[NDIMS], start[NDIMS], dim_len[NDIMS], iz;
+ size_t i;
+ double vxsz;
+ nc_type type;
+ int varid, dimids[NDIMS], irregular, ndims;
+ int err = NC_NOERR;
+ res_T res = RES_OK;
+ ASSERT(nc && voxel_size && lower_pos && is_z_irregular);
+
+ err = nc_inq_varid(nc, "VLEV", &varid);
+ if(err != NC_NOERR) {
+ fprintf(stderr, "Could not inquire the 'VLEV' variable -- %s\n",
+ ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_varndims(nc, varid, &ndims));
+ if(ndims != NDIMS) {
+ fprintf(stderr, "The dimension of the 'VLEV' variable must be %i .\n",
+ NDIMS);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_vardimid(nc, varid, dimids));
+ FOR_EACH(i, 0, NDIMS) NC(inq_dimlen(nc, dimids[i], len+i));
+
+ NC(inq_vartype(nc, varid, &type));
+ if(type != NC_DOUBLE && type != NC_FLOAT) {
+ fprintf(stderr, "The type of 'VLEV' variable must be float or double.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ mem = mem_alloc(len[Z]*sizeof_nctype(type));
+ if(!mem) {
+ fprintf(stderr, "Could not allocate memory for the 'VLEV' variable.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ start[X] = 0; dim_len[X] = 1;
+ start[Y] = 0; dim_len[Y] = 1;
+ start[Z] = 0; dim_len[Z] = len[Z];
+
+ NC(get_vara(nc, varid, start, dim_len, mem));
+
+ #define AT(Mem, Id) (type == NC_FLOAT ? ((float*)Mem)[Id] : ((double*)Mem)[Id])
+
+ /* Assume that the Z dimension starts from 0 */
+ vxsz = AT(mem, 0) * 2.0;
+ if(vxsz <= 0) {
+ fprintf(stderr, "The 'VLEV' variable can't have negative or null values.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check if the Z dimension is regular */
+ FOR_EACH(iz, 1, len[Z]) {
+ if(!eq_eps(AT(mem, iz) - AT(mem, iz-1), vxsz, 1.e-6)) break;
+ }
+ irregular = (iz != len[Z]);
+
+ /* Setup the size of the voxel */
+ pvxsz = mem_alloc(sizeof(double) * (irregular ? len[Z] : 1));
+ if(!pvxsz) {
+ fprintf(stderr,
+ "Could not allocate memory for the voxel size along the Z dimension.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ pvxsz[0] = vxsz;
+
+ if(irregular) {
+ FOR_EACH(iz, 1, len[Z]) {
+ pvxsz[iz] = (AT(mem, iz) - (AT(mem, iz-1) + pvxsz[iz-1]*0.5)) * 2.0;
+ if(pvxsz[iz] <= 0) {
+ fprintf(stderr,
+ "The 'VLEV' variable can't have negative or null values.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ }
+
+ /* Check that the others columns have the same voxel size */
+ if(check) {
+ size_t ix, iy;
+
+ mem2 = mem_alloc(len[Z]*sizeof_nctype(type));
+ if(!mem2) {
+ fprintf(stderr,
+ "Could not allocate memory to check the 'VLEV' variable.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(iy, 0, len[Y]) {
+ start[Y] = iy;
+ FOR_EACH(ix, 1, len[X]) {
+ start[X] = ix;
+ NC(get_vara(nc, varid, start, dim_len, mem2));
+ FOR_EACH(iz, 0, len[Z]) {
+ if(!eq_eps(AT(mem2, iz), AT(mem, iz), 1.e-6)) {
+ fprintf(stderr,
+ "The Z columns of the 'VLEV' variable must be the same.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ }
}
}
+ #undef AT
+
+ *lower_pos = 0;
+ *voxel_size = pvxsz;
+ *is_z_irregular = (int8_t)irregular;
+
+exit:
+ if(mem2) mem_rm(mem2);
+ if(mem) mem_rm(mem);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+write_data(int nc, const char* var_name, FILE* stream)
+{
+ enum { TIME, Z, X, Y, NDIMS };
+ char* mem = NULL;
+ size_t start[NDIMS], dim_len[NDIMS], len[NDIMS], i, grid_len;
+ nc_type type;
+ int varid, ndims, dimids[NDIMS];
+ int err = NC_NOERR;
+ res_T res = RES_OK;
+ CHK(var_name);
+
+ err = nc_inq_varid(nc, var_name, &varid);
+ if(err != NC_NOERR) {
+ fprintf(stderr, "Could not inquire the '%s' variable -- %s\n",
+ var_name, ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
- NC(inq_varid(ncid, "RVT", &id));
- NC(inq_var(ncid, id, NULL, &type, &ndims, dimid, NULL));
- CHK(type == NC_FLOAT);
- CHK(ndims == 4);
- NC(inq_dim(ncid, dimid[0], NULL, len+0));
- NC(inq_dim(ncid, dimid[1], NULL, len+1));
- NC(inq_dim(ncid, dimid[2], NULL, len+2));
- NC(inq_dim(ncid, dimid[3], NULL, len+3));
- startid[0] = 0;
- startid[1] = 0;
- startid[2] = 0;
- startid[3] = 0;
+ NC(inq_varndims(nc, varid, &ndims));
+ if(ndims != NDIMS) {
+ fprintf(stderr, "The dimension of the '%s' variable must be %i .\n",
+ var_name, NDIMS);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ NC(inq_vardimid(nc, varid, dimids));
+ FOR_EACH(i, 0, NDIMS) NC(inq_dimlen(nc, dimids[i], len+i));
+
+ NC(inq_vartype(nc, varid, &type));
+ if(type != NC_DOUBLE && type != NC_FLOAT) {
+ fprintf(stderr,
+ "The type of the '%s' variable cannot be %s. Expecting floating point data.\n",
+ var_name, nctype_to_str(type));
+ res = RES_BAD_ARG;
+ goto error;
+ }
- CHK(mem = mem_alloc(len[0]*len[1]*len[2]*len[3]*sizeof_nctype(type)));
- NC(get_vara(ncid, id, startid, len, mem));
- FOR_EACH(ival, 0, len[0]*len[1]*len[2]*len[3]) {
- printf("%g\n", ((float*)mem)[ival]);
+ grid_len = len[X]*len[Y]*len[Z];
+ mem = mem_alloc(grid_len*sizeof_nctype(type));
+ if(!mem) {
+ fprintf(stderr, "Could not allocate memory for the '%s' variable.\n",
+ var_name);
+ res = RES_MEM_ERR;
+ goto error;
}
+ start[X] = 0; dim_len[X] = len[X];
+ start[Y] = 0; dim_len[Y] = len[Y];
+ start[Z] = 0; dim_len[Z] = len[Z];
+
+ FOR_EACH(i, 0, len[TIME]) {
+ start[TIME] = i;
+ dim_len[TIME] = 1;
+ NC(get_vara(nc, varid, start, dim_len, mem));
+ if(fwrite(mem, sizeof_nctype(type), grid_len, stream) != grid_len) {
+ fprintf(stderr, "Error writing data of the '%s' variable -- '%s'.\n",
+ var_name, strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+ }
exit:
if(mem) mem_rm(mem);
- if(ncid != NCID_INVALID) NC(close(ncid));
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Program
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ FILE* stream = stdout;
+ struct grid grid = GRID_NULL;
+ struct args args = ARGS_DEFAULT;
+ int err = 0;
+ int err_nc = NC_NOERR;
+ int nc = INVALID_ID;
+ res_T res = RES_OK;
+
+ res = args_init(&args, argc, argv);
+ if(res != RES_OK) goto error;
+ if(args.quit) goto exit;
+
+ if(args.output) {
+ if(access(args.output, F_OK) == 0) {
+ fprintf(stderr, "The OUTPUT file '%s' already exists.\n", args.output);
+ res = RES_IO_ERR;
+ goto error;
+ }
+ stream = fopen(args.output, "w");
+ if(!stream) {
+ fprintf(stderr, "Could not create the OUTPUT file '%s'.\n", args.output);
+ res = RES_IO_ERR;
+ goto error;
+ }
+ }
+
+ err_nc = nc_open(args.input, NC_WRITE|NC_MMAP, &nc);
+ if(err_nc != NC_NOERR) {
+ fprintf(stderr, "Error opening file `%s' -- %s.\n",
+ args.input, ncerr_to_str(err_nc));
+ goto exit;
+ }
+
+ #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0
+ CALL(setup_definition(nc, &grid.nx, &grid.ny, &grid.nz, &grid.ntimes));
+ CALL(setup_regular_dimension
+ (nc, "W_E_direction", &grid.vxsz_x, grid.lower+0, args.check));
+ CALL(setup_regular_dimension
+ (nc, "S_N_direction", &grid.vxsz_y, grid.lower+1, args.check));
+ CALL(setup_Z_dimension
+ (nc, &grid.vxsz_z, grid.lower+2, &grid.is_z_irregular, args.check));
+ #undef CALL
+
+ if(!args.no_output) {
+ CHK(fwrite(&grid.is_z_irregular, sizeof(int8_t), 1, stream) == 1);
+ CHK(fwrite(&grid.nx, sizeof(int32_t), 1, stream) == 1);
+ CHK(fwrite(&grid.ny, sizeof(int32_t), 1, stream) == 1);
+ CHK(fwrite(&grid.nz, sizeof(int32_t), 1, stream) == 1);
+ CHK(fwrite(&grid.ntimes, sizeof(int32_t), 1, stream) == 1);
+ CHK(fwrite(grid.lower, sizeof(double), 3, stream) == 3);
+ CHK(fwrite(&grid.vxsz_x, sizeof(double), 1, stream) == 1);
+ CHK(fwrite(&grid.vxsz_y, sizeof(double), 1, stream) == 1);
+ if(!grid.is_z_irregular) {
+ CHK(fwrite(grid.vxsz_z, sizeof(double), 1, stream) == 1);
+ } else {
+ CHK(fwrite(grid.vxsz_z, sizeof(double), (size_t)grid.nz, stream)
+ == (size_t)grid.nz);
+ }
+ write_data(nc, "RVT", stream);
+ write_data(nc, "RCT", stream);
+ }
+
+exit:
+ grid_release(&grid);
+ if(nc != INVALID_ID) NC(close(nc));
return err;
error:
err = -1;
