[GRASS-SVN] r52062 - in grass/trunk/raster3d/r3.out.netcdf: . test_suite
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 13 15:14:08 PDT 2012
Author: huhabla
Date: 2012-06-13 15:14:07 -0700 (Wed, 13 Jun 2012)
New Revision: 52062
Modified:
grass/trunk/raster3d/r3.out.netcdf/main.c
grass/trunk/raster3d/r3.out.netcdf/test_suite/test.r3.out.netcdf.sh
Log:
Missing value and FillValue support added. Projection information is now
optional, since some software get confused by it (ParaView).
Modified: grass/trunk/raster3d/r3.out.netcdf/main.c
===================================================================
--- grass/trunk/raster3d/r3.out.netcdf/main.c 2012-06-13 16:29:53 UTC (rev 52061)
+++ grass/trunk/raster3d/r3.out.netcdf/main.c 2012-06-13 22:14:07 UTC (rev 52062)
@@ -59,11 +59,12 @@
/* structs */
typedef struct
{
- struct Option *input, *output;
- struct Flag *mask;
+ struct Option *input, *output, *null;
+ struct Flag *mask, *proj;
} paramType;
RASTER3D_Map *map;
+
paramType param;
/*---------------------------------------------------------------------------*/
@@ -76,7 +77,9 @@
if (map != NULL) {
/* should unopen map here! */
if (!Rast3d_close(map))
- G_fatal_error(_("Unable to close 3D raster map while catching error: %s"), errorMsg);
+ G_fatal_error(_
+ ("Unable to close 3D raster map while catching error: %s"),
+ errorMsg);
}
G_fatal_error("%s", errorMsg);
}
@@ -93,6 +96,20 @@
param.output->key = "output";
param.output->description = _("Name for netcdf output file");
+ param.null = G_define_option();
+ param.null->key = "null";
+ param.null->type = TYPE_DOUBLE;
+ param.null->required = NO;
+ param.null->multiple = NO;
+ param.null->description =
+ _
+ ("The value to be used for null values, default is the NetCDF standard");
+
+ param.proj = G_define_flag();
+ param.proj->key = 'p';
+ param.proj->description =
+ _("Export projection information as wkt and proj4 parameter");
+
param.mask = G_define_flag();
param.mask->key = 'm';
param.mask->description =
@@ -100,7 +117,7 @@
}
static void write_netcdf_header(int ncid, RASTER3D_Region * region,
- int *varid)
+ int *varid, char write_proj, char *null)
{
int retval, typeIntern, time;
size_t i;
@@ -111,6 +128,7 @@
int lat_varid = 0, lon_varid = 0, time_varid = 0;
int dimids[NDIMS];
struct Cell_head window;
+ double min, max;
/* global attributes */
if ((retval =
@@ -125,7 +143,7 @@
G_get_window(&window);
/* Projection parameter */
- if (window.proj != PROJECTION_XY) {
+ if (window.proj != PROJECTION_XY && write_proj) {
struct Key_Value *pkv, *ukv;
struct pj_info pjinfo;
char *proj4, *proj4mod;
@@ -164,6 +182,7 @@
G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
else
proj4mod = G_store(proj4);
+
pj_dalloc(proj4);
if ((retval =
@@ -180,7 +199,7 @@
typeIntern = Rast3d_tile_type_map(map);
if (window.proj == PROJECTION_LL) {
- /* X-Axis */
+ /* X-Axis */
if ((retval = nc_def_dim(ncid, LON_NAME, region->cols, &lon_dimid)))
ERR(retval);
@@ -192,30 +211,28 @@
nc_put_att_text(ncid, lon_varid, UNITS, strlen(DEGREES_EAST),
DEGREES_EAST)))
ERR(retval);
- if ((retval =
- nc_put_att_text(ncid, lon_varid, LONG_NAME,
- strlen(LON_LONG_NAME), LON_LONG_NAME)))
+ if ((retval = nc_put_att_text(ncid, lon_varid, LONG_NAME,
+ strlen(LON_LONG_NAME), LON_LONG_NAME)))
ERR(retval);
if ((retval =
nc_put_att_text(ncid, lon_varid, STANDARD_NAME, strlen(LON_NAME),
LON_NAME)))
ERR(retval);
-
- /* Y-Axis */
+
+ /* Y-Axis */
if ((retval = nc_def_dim(ncid, LAT_NAME, region->rows, &lat_dimid)))
ERR(retval);
-
+
if ((retval =
nc_def_var(ncid, LAT_NAME, NC_FLOAT, 1, &lat_dimid, &lat_varid)))
ERR(retval);
-
+
if ((retval =
nc_put_att_text(ncid, lat_varid, UNITS, strlen(DEGREES_NORTH),
DEGREES_NORTH)))
ERR(retval);
- if ((retval =
- nc_put_att_text(ncid, lat_varid, LONG_NAME,
- strlen(LAT_LONG_NAME), LAT_LONG_NAME)))
+ if ((retval = nc_put_att_text(ncid, lat_varid, LONG_NAME,
+ strlen(LAT_LONG_NAME), LAT_LONG_NAME)))
ERR(retval);
if ((retval =
nc_put_att_text(ncid, lat_varid, STANDARD_NAME, strlen(LAT_NAME),
@@ -223,7 +240,7 @@
ERR(retval);
}
else {
- /* X-Axis */
+ /* X-Axis */
if ((retval = nc_def_dim(ncid, X_NAME, region->cols, &lon_dimid)))
ERR(retval);
@@ -231,37 +248,35 @@
nc_def_var(ncid, X_NAME, NC_FLOAT, 1, &lon_dimid, &lon_varid)))
ERR(retval);
- if ((retval =
- nc_put_att_text(ncid, lon_varid, UNITS, strlen("meter"),
- "meter")))
+ if ((retval = nc_put_att_text(ncid, lon_varid, UNITS, strlen("meter"),
+ "meter")))
ERR(retval);
if ((retval =
nc_put_att_text(ncid, lon_varid, LONG_NAME, strlen(X_LONG_NAME),
X_LONG_NAME)))
ERR(retval);
- if ((retval =
- nc_put_att_text(ncid, lon_varid, STANDARD_NAME,
- strlen(X_STANDARD_NAME), X_STANDARD_NAME)))
+ if ((retval = nc_put_att_text(ncid, lon_varid, STANDARD_NAME,
+ strlen(X_STANDARD_NAME),
+ X_STANDARD_NAME)))
ERR(retval);
- /* Y-Axis */
+ /* Y-Axis */
if ((retval = nc_def_dim(ncid, Y_NAME, region->rows, &lat_dimid)))
ERR(retval);
-
+
if ((retval =
nc_def_var(ncid, Y_NAME, NC_FLOAT, 1, &lat_dimid, &lat_varid)))
ERR(retval);
-
- if ((retval =
- nc_put_att_text(ncid, lat_varid, UNITS, strlen("meter"),
- "meter")))
+
+ if ((retval = nc_put_att_text(ncid, lat_varid, UNITS, strlen("meter"),
+ "meter")))
ERR(retval);
if ((retval =
nc_put_att_text(ncid, lat_varid, LONG_NAME, strlen(Y_LONG_NAME),
Y_LONG_NAME)))
ERR(retval);
- if ((retval =
- nc_put_att_text(ncid, lat_varid, STANDARD_NAME,
- strlen(Y_STANDARD_NAME), Y_STANDARD_NAME)))
+ if ((retval = nc_put_att_text(ncid, lat_varid, STANDARD_NAME,
+ strlen(Y_STANDARD_NAME),
+ Y_STANDARD_NAME)))
ERR(retval);
}
@@ -297,7 +312,7 @@
G_read_raster3d_timestamp(map->fileName, map->mapset, &ts);
- /* Days since datum in ISO norm*/
+ /* Days since datum in ISO norm */
if (datetime_is_absolute(&ts.dt[0])) {
is_absolute_time = 1;
G_snprintf(time_unit, 1024,
@@ -310,7 +325,6 @@
G_snprintf(time_unit, 1024, "%s",
Rast3d_get_vertical_unit(map));
}
-
}
else {
G_snprintf(time_unit, 1024, "%s since %s",
@@ -342,7 +356,8 @@
else {
if ((retval =
nc_def_dim(ncid, Z_NAME, region->depths, &time_dimid)))
- ERR(retval);;
+ ERR(retval);
+ ;
if ((retval =
nc_def_var(ncid, Z_NAME, NC_FLOAT, 1, &time_dimid,
&time_varid)))
@@ -382,7 +397,7 @@
strlen(Z_STANDARD_NAME), Z_STANDARD_NAME)))
ERR(retval);
}
-
+
/* z - axis orientation */
if ((retval =
nc_put_att_text(ncid, time_varid, "positive", strlen("up"), "up")))
@@ -406,19 +421,70 @@
dimids[1] = lat_dimid;
dimids[2] = lon_dimid;
+
+ Rast3d_range_load(map);
+ Rast3d_range_min_max(map, &min, &max);
+
if (typeIntern == FCELL_TYPE) {
if ((retval =
nc_def_var(ncid, param.input->answer, NC_FLOAT, NDIMS, dimids,
varid)))
ERR(retval);
+ /* Set the range values */
+ float fmin = min;
+
+ float fmax = max;
+
+ if ((retval =
+ nc_put_att_float(ncid, *varid, "valid_min", NC_FLOAT, 1, &fmin)))
+ ERR(retval);
+ if ((retval =
+ nc_put_att_float(ncid, *varid, "valid_max", NC_FLOAT, 1, &fmax)))
+ ERR(retval);
+
+ if (null) {
+ float null_val = (float)atof(null);
+
+ if ((retval =
+ nc_put_att_float(ncid, *varid, "missing_value", NC_FLOAT, 1,
+ &null_val)))
+ ERR(retval);
+ if ((retval =
+ nc_put_att_float(ncid, *varid, "_FillValue", NC_FLOAT, 1,
+ &null_val)))
+ ERR(retval);
+ }
}
else {
if ((retval =
nc_def_var(ncid, param.input->answer, NC_DOUBLE, NDIMS, dimids,
varid)))
ERR(retval);
+ /* Set the range values */
+ if ((retval =
+ nc_put_att_double(ncid, *varid, "valid_min", NC_DOUBLE, 1,
+ &min)))
+ ERR(retval);
+ if ((retval =
+ nc_put_att_double(ncid, *varid, "valid_max", NC_DOUBLE, 1,
+ &max)))
+ ERR(retval);
+
+ if (null) {
+ double null_val = (double)atof(null);
+
+ if ((retval =
+ nc_put_att_double(ncid, *varid, "missing_value", NC_DOUBLE,
+ 1, &null_val)))
+ ERR(retval);
+ if ((retval =
+ nc_put_att_double(ncid, *varid, "_FillValue", NC_DOUBLE, 1,
+ &null_val)))
+ ERR(retval);
+ }
}
- if (window.proj != PROJECTION_XY) {
+
+ if (window.proj != PROJECTION_XY && write_proj) {
if ((retval =
nc_put_att_text(ncid, *varid, "grid_mapping", strlen("crs"),
"crs")))
@@ -446,7 +512,7 @@
for (i = 0; i < region->depths; i++) {
if (is_time) {
- c = region->bottom + i * region->tb_res;
+ c = i * region->tb_res;
time = (int)c;
nc_put_var1_int(ncid, time_varid, &i, &time);
}
@@ -508,11 +574,9 @@
int main(int argc, char *argv[])
{
RASTER3D_Region region;
-
int ncid;
int retval;
int varid;
-
int changemask = 0;
struct GModule *module;
@@ -556,7 +620,8 @@
if ((retval = nc_create(param.output->answer, NC_CLOBBER, &ncid)))
ERR(retval);
- write_netcdf_header(ncid, ®ion, &varid);
+ write_netcdf_header(ncid, ®ion, &varid, param.proj->answer,
+ param.null->answer);
/*if requested set the Mask on */
if (param.mask->answer) {
@@ -571,7 +636,7 @@
write_netcdf_data(ncid, ®ion, varid);
- /*We set the Mask off, if it was off bevor */
+ /*We set the Mask off, if it was off before */
if (param.mask->answer) {
if (Rast3d_mask_file_exists())
if (Rast3d_mask_is_on(map) && changemask)
Modified: grass/trunk/raster3d/r3.out.netcdf/test_suite/test.r3.out.netcdf.sh
===================================================================
--- grass/trunk/raster3d/r3.out.netcdf/test_suite/test.r3.out.netcdf.sh 2012-06-13 16:29:53 UTC (rev 52061)
+++ grass/trunk/raster3d/r3.out.netcdf/test_suite/test.r3.out.netcdf.sh 2012-06-13 22:14:07 UTC (rev 52062)
@@ -6,14 +6,14 @@
# @preprocess step of this test. We generate
# voxel data with r3.mapcalc. The region setting
# should work for UTM and LL test locations
-g.region s=-90 n=90 w=-180 e=180 b=0 t=50 res=10 res3=10 -p3
+g.region s=-90 n=90 w=-180 e=180 b=0 t=5 res=10 res3=10 tbres=1 -p3
# We create several (float, double, null value) voxel map
# with value = col + row + depth.
r3.mapcalc --o expr="volume_float = float(col() + row() + depth())"
r3.mapcalc --o expr="volume_double = double(col() + row() + depth())"
r3.mapcalc --o expr="volume_time_double = double(col() + row() + depth())"
r3.mapcalc --o expr="volume_time_float = float(col() + row() + depth())"
-r3.timestamp map=volume_time_double date='1 Jan 2001/5 Jan 2001'
+r3.timestamp map=volume_time_double date='1 Jan 1900/5 Jan 1900'
r3.support map=volume_time_double vunit="days"
r3.timestamp map=volume_time_float date='5 seconds/10 seconds'
r3.support map=volume_time_float vunit="seconds"
@@ -21,13 +21,13 @@
r3.out.netcdf --o input=volume_float output=test_float.nc
#r3.info volume_float
#ncdump -h test_float.nc
-r3.out.netcdf --o input=volume_double output=test_double.nc
+r3.out.netcdf --o null=-100 input=volume_double output=test_double.nc
#r3.info volume_double
#ncdump -h test_double.nc
-r3.out.netcdf --o input=volume_time_double output=test_time_double.nc
+r3.out.netcdf --o -p input=volume_time_double output=test_time_double.nc
#r3.info volume_time_double
#ncdump -h test_time_double.nc
-r3.out.netcdf --o input=volume_time_float output=test_time_float.nc
+r3.out.netcdf --o -p null=-1000 input=volume_time_float output=test_time_float.nc
#r3.info volume_time_float
#ncdump -h test_time_float.nc
More information about the grass-commit
mailing list