NetCDF library support integrated in configuration.

Projection specific metadata generation implemented in r3.out.netcdf.
Added r3.out.netcdf to 3D raster Makefile.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@51634 15284696-431f-4ddb-bdfa-cd5b030d7da7

configure.in

@@ -282,6 +282,12 @@ AC_ARG_WITH(wxwidgets,
                           e.g. '--with-wxwidgets=/usr/local/bin/wx-config',
                           default: no)],, with_wxwidgets="no")
 
+AC_ARG_WITH(netcdf,
+[  --with-netcdf[=path/nc-config]
+                          enable NetCDF support (nc-config with path,
+                          e.g. '--with-nc=/usr/local/bin/nc-config',
+                          default: no)],, with_netcdf="no")
+
 # With includes and libs options
 
 AC_ARG_WITH(geos,
@@ -826,12 +832,60 @@ AC_SUBST(LIBLAS_CFLAGS)
 AC_SUBST(LIBLAS_INC)
 AC_SUBST(USE_LIBLAS)
 
+# NetCDF option
+
+AC_MSG_CHECKING(whether to use NetCDF)
+
+NETCDF_LIBS=
+NETCDF_CFLAGS=
+USE_NETCDF=
+
+if test "`basename xx/$with_netcdf`" = "nc-config" ; then
+  NETCDF_CONFIG="$with_netcdf"
+fi
+
+if test "$with_netcdf" = "no" ; then
+  AC_MSG_RESULT(no)
+else
+  AC_MSG_RESULT(yes)
+  AC_PATH_PROG(NETCDF_CONFIG, nc-config, no)
+
+  if test "$NETCDF_CONFIG" = "no" ; then
+    AC_MSG_ERROR([*** couldn't find nc-config])
+  fi
+
+  if test "$NETCDF_CONFIG" != "" ; then
+    NETCDF_LIBS=`"$NETCDF_CONFIG" --libs`
+    NETCDF_CFLAGS=`"$NETCDF_CONFIG" --cflags`
+    USE_NETCDF=1
+  fi
+
+  NETCDF=
+  ac_save_libs="$LIBS"
+  ac_save_cflags="$CFLAGS"
+  LIBS="$LIBS $NETCDF_LIBS"
+  CFLAGS="$CFLAGS $NETCDF_CFLAGS"
+  AC_TRY_LINK([#include <netcdf.h>],[nc_create("foo", NC_CLOBBER, NULL);],,[
+  AC_TRY_LINK([#include <netcdf.h>],[nc_create("foo", NC_CLOBBER, NULL);],NETCDF_LIBS="$NETCDF_LIBS",[
+  AC_MSG_ERROR([*** Unable to locate NetCDF library.])
+  ])
+  ])
+  LIBS=${ac_save_libs}
+  CFLAGS=${ac_save_cflags}
+
+  AC_DEFINE(HAVE_NETCDF)
+fi
+
+AC_SUBST(NETCDF_LIBS)
+AC_SUBST(NETCDF_CFLAGS)
+AC_SUBST(USE_NETCDF)
+
 # GEOS option
 
 AC_MSG_CHECKING(whether to use GEOS)
 
 GEOS_LIBS=
-GGEOS_CFLAGS=
+GEOS_CFLAGS=
 USE_GEOS=
 
 # FIXME: "quote" $with_geos ?
@@ -1975,6 +2029,7 @@ LOC_MSG_USE(FFMPEG support,USE_FFMPEG)
 LOC_MSG_USE(FFTW support,USE_FFTW)
 LOC_MSG_USE(FreeType support,USE_FREETYPE)
 LOC_MSG_USE(GDAL support,USE_GDAL)
+LOC_MSG_USE(NETCDF support,USE_NETCDF)
 LOC_MSG_USE(GEOS support,USE_GEOS)
 LOC_MSG_USE(LAPACK support,USE_LAPACK)
 LOC_MSG_USE(Large File support (LFS), USE_LARGEFILES)

include/Make/Platform.make.in

@@ -183,6 +183,11 @@ GDALCFLAGS          = @GDAL_CFLAGS@
 USE_GDAL            = @USE_GDAL@
 USE_OGR             = @USE_OGR@
 
+#NetCDF
+NETCDFLIBS           = @NETCDF_LIBS@
+NETCDFCFLAGS        = @NETCDF_CFLAGS@    
+USE_NETCDF          = @USE_NETCDF@
+
 #LAS LiDAR through libLAS
 LASLIBS             = @LIBLAS_LIBS@
 LASCFLAGS           = @LIBLAS_CFLAGS@

raster3d/Makefile

@@ -10,6 +10,7 @@ SUBDIRS = \
 	r3.mkdspf \
 	r3.null \
 	r3.out.ascii \
+	r3.out.netcdf \
 	r3.out.v5d \
 	r3.out.vtk \
 	r3.retile \

raster3d/r3.out.netcdf/Makefile

@@ -2,13 +2,14 @@ MODULE_TOPDIR = ../..
 
 PGM=r3.out.netcdf
 
-LIBES = $(RASTER3DLIB) $(GISLIB) $(GPROJLIB)
-DEPENDENCIES = $(RASTER3DDEP) $(GISDEP) $(GPROJDEP)
+LIBES = $(RASTER3DLIB) $(GPROJLIB) $(GISLIB) $(GDALLIBS) $(NETCDFLIBS) $(PROJLIB) $(DATETIMELIB)
+DEPENDENCIES = $(RASTER3DDEP) $(GISDEP) $(GPROJDEP) $(DATETIMEDEP)
 
 EXTRA_INC = $(PROJINC)
-
-EXTRA_LIBS = -lnetcdf
+EXTRA_CFLAGS = $(NETCDF_CFLAGS)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 
+ifneq ($(strip $(NETCDFLIBS)),)
 default: cmd
+endif

raster3d/r3.out.netcdf/main.c

@@ -14,16 +14,18 @@
  *   	    	for details.
  * 
  * TODO: Add time zone support to time variable
- * TODO: Implement support for coordinate reference system defined here:
+ * TODO: Implement better support for CF coordinate reference system defined here:
  *       http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#coordinate-system
  *       https://cf-pcmdi.llnl.gov/trac/wiki/Cf2CrsWkt
  *       http://trac.osgeo.org/gdal/wiki/NetCDF_ProjectionTestingStatus
  *
  *****************************************************************************/
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <netcdf.h>
+#include <grass/datetime.h>
 #include <grass/gis.h>
 #include <grass/raster3d.h>
 #include <grass/glocale.h>
@@ -45,6 +47,7 @@
 #define Y_STANDARD_NAME "projection_y_coordinate"
 #define Z_NAME "z"
 #define Z_LONG_NAME "z coordinate of projection"
+#define Z_STANDARD_NAME "projection_z_coordinate"
 #define UNITS "units"
 #define DEGREES_EAST "degrees_east"
 #define DEGREES_NORTH "degrees_north"
@@ -54,7 +57,8 @@
 #define ERR(e) {fatalError(nc_strerror(e));}
 
 /* structs */
-typedef struct {
+typedef struct
+{
     struct Option *input, *output;
     struct Flag *mask;
 } paramType;
@@ -70,9 +74,9 @@ paramType param;
 static void fatalError(const char *errorMsg)
 {
     if (map != NULL) {
-        /* should unopen map here! */
-        if (!Rast3d_close(map))
-            fatalError(_("Unable to close 3D raster map"));
+	/* 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("%s", errorMsg);
 }
@@ -83,22 +87,25 @@ static void fatalError(const char *errorMsg)
  */
 static void setParams()
 {
-    param.input = G_define_standard_option(G_OPT_R3_INPUTS);
-    
+    param.input = G_define_standard_option(G_OPT_R3_INPUT);
+
     param.output = G_define_standard_option(G_OPT_F_OUTPUT);
     param.output->key = "output";
     param.output->description = _("Name for netcdf output file");
 
     param.mask = G_define_flag();
     param.mask->key = 'm';
-    param.mask->description = _("Use 3D raster mask (if exists) with input map");
+    param.mask->description =
+	_("Use 3D raster mask (if exists) with input map");
 }
 
-static void write_netcdf_header(int ncid, RASTER3D_Region *region, int *varid)
+static void write_netcdf_header(int ncid, RASTER3D_Region * region,
+				int *varid)
 {
     int retval, typeIntern, time;
     size_t i;
     int is_time = 0;
+    int is_absolute_time = 0;
     float c;
     int lat_dimid = 0, lon_dimid = 0, time_dimid = 0;
     int lat_varid = 0, lon_varid = 0, time_varid = 0;
@@ -106,179 +113,353 @@ static void write_netcdf_header(int ncid, RASTER3D_Region *region, int *varid)
     struct Cell_head window;
 
     /* global attributes */
-    if ((retval = nc_put_att_text(ncid, NC_GLOBAL, "Conventions", strlen(CF_SUPPORT), CF_SUPPORT))) ERR(retval);
-    if ((retval = nc_put_att_text(ncid, NC_GLOBAL, "history", strlen(HISTORY_TEXT), HISTORY_TEXT))) ERR(retval);
+    if ((retval =
+	 nc_put_att_text(ncid, NC_GLOBAL, "Conventions", strlen(CF_SUPPORT),
+			 CF_SUPPORT)))
+	ERR(retval);
+    if ((retval =
+	 nc_put_att_text(ncid, NC_GLOBAL, "history", strlen(HISTORY_TEXT),
+			 HISTORY_TEXT)))
+	ERR(retval);
 
     G_get_window(&window);
-     
+
     /* Projection parameter */
-    if(window.proj != PROJECTION_XY) {
-        struct Key_Value *pkv, *ukv;
-        struct pj_info pjinfo;
-        char *proj4, *proj4mod;
-        const char *unfact;
+    if (window.proj != PROJECTION_XY) {
+	struct Key_Value *pkv, *ukv;
+	struct pj_info pjinfo;
+	char *proj4, *proj4mod;
+	const char *unfact;
 	int crs_dimid = 0, crs_varid;
 
-        if ((retval = nc_def_var(ncid, "crs", NC_CHAR, 0, &crs_dimid, &crs_varid))) ERR(retval);
+	if ((retval =
+	     nc_def_var(ncid, "crs", NC_CHAR, 0, &crs_dimid, &crs_varid)))
+	    ERR(retval);
+
+	pkv = G_get_projinfo();
+	ukv = G_get_projunits();
 
-     	pkv = G_get_projinfo();
-        ukv = G_get_projunits();
- 
-        pj_get_kv(&pjinfo, pkv, ukv);
-        proj4 = pj_get_def(pjinfo.pj, 0);
-        pj_free(pjinfo.pj);
+	pj_get_kv(&pjinfo, pkv, ukv);
+	proj4 = pj_get_def(pjinfo.pj, 0);
+	pj_free(pjinfo.pj);
 
 #ifdef HAVE_OGR
 	/* We support the CF suggestion crs_wkt and the gdal spatil_ref attribute */
-        if ((retval = nc_put_att_text(ncid, crs_varid, "crs_wkt", strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)), GPJ_grass_to_wkt(pkv, ukv, 0, 0)))) ERR(retval);
-        if ((retval = nc_put_att_text(ncid, crs_varid, "spatial_ref", strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)), GPJ_grass_to_wkt(pkv, ukv, 0, 0)))) ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, crs_varid, "crs_wkt",
+			     strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)),
+			     GPJ_grass_to_wkt(pkv, ukv, 0, 0))))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, crs_varid, "spatial_ref",
+			     strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)),
+			     GPJ_grass_to_wkt(pkv, ukv, 0, 0))))
+	    ERR(retval);
 #endif
-        /* Code from g.proj:
-         * GRASS-style PROJ.4 strings don't include a unit factor as this is
-         * handled separately in GRASS - must include it here though */
-        unfact = G_find_key_value("meters", ukv);
-        if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
-            G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
-        else
-            proj4mod = G_store(proj4);
-        pj_dalloc(proj4); 
-
-        if ((retval = nc_put_att_text(ncid, crs_varid, "crs_proj4", strlen(proj4mod), proj4mod))) ERR(retval);
-
-        if(pkv)
-            G_free_key_value(pkv);
-        if(ukv)
-            G_free_key_value(ukv);
+	/* Code from g.proj:
+	 * GRASS-style PROJ.4 strings don't include a unit factor as this is
+	 * handled separately in GRASS - must include it here though */
+	unfact = G_find_key_value("meters", ukv);
+	if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
+	    G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
+	else
+	    proj4mod = G_store(proj4);
+	pj_dalloc(proj4);
+
+	if ((retval =
+	     nc_put_att_text(ncid, crs_varid, "crs_proj4", strlen(proj4mod),
+			     proj4mod)))
+	    ERR(retval);
+
+	if (pkv)
+	    G_free_key_value(pkv);
+	if (ukv)
+	    G_free_key_value(ukv);
     }
 
     typeIntern = Rast3d_tile_type_map(map);
-    
+
     if (window.proj == PROJECTION_LL) {
-        if ((retval = nc_def_dim(ncid, LON_NAME, region->cols, &lon_dimid))) ERR(retval);
-        if ((retval = nc_def_dim(ncid, LAT_NAME, region->rows, &lat_dimid))) ERR(retval);
+        /* X-Axis */
+	if ((retval = nc_def_dim(ncid, LON_NAME, region->cols, &lon_dimid)))
+	    ERR(retval);
+
+	if ((retval =
+	     nc_def_var(ncid, LON_NAME, NC_FLOAT, 1, &lon_dimid, &lon_varid)))
+	    ERR(retval);
+
+	if ((retval =
+	     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)))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, lon_varid, STANDARD_NAME, strlen(LON_NAME),
+			     LON_NAME)))
+	    ERR(retval);
+        
+        /* Y-Axis */
+	if ((retval = nc_def_dim(ncid, LAT_NAME, region->rows, &lat_dimid)))
+	    ERR(retval);
         
-        if ((retval = nc_def_var(ncid, LON_NAME, NC_FLOAT, 1, &lon_dimid, &lon_varid))) 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, 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))) ERR(retval);
-        if ((retval = nc_put_att_text(ncid, lon_varid, STANDARD_NAME, strlen(LON_NAME), LON_NAME))) 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))) ERR(retval);
-        if ((retval = nc_put_att_text(ncid, lat_varid, STANDARD_NAME, strlen(LAT_NAME), LAT_NAME))) ERR(retval);
-    } else {
-        if ((retval = nc_def_dim(ncid, X_NAME, region->cols, &lon_dimid))) ERR(retval);
-        if ((retval = nc_def_dim(ncid, Y_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_def_var(ncid, X_NAME, NC_FLOAT, 1, &lon_dimid, &lon_varid))) 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, 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))) ERR(retval);
-        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))) 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)))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, lat_varid, STANDARD_NAME, strlen(LAT_NAME),
+			     LAT_NAME)))
+	    ERR(retval);
     }
+    else {
+        /* X-Axis */
+	if ((retval = nc_def_dim(ncid, X_NAME, region->cols, &lon_dimid)))
+	    ERR(retval);
+
+	if ((retval =
+	     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")))
+	    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)))
+	    ERR(retval);
+        /* 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")))
+	    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)))
+	    ERR(retval);
+    }
+
+    /* We set the vertical axis and its unit. Units can be spatial
+     * or temporal. Temporal can be absolute or relative. 
+     */
+
     is_time = 0;
-    
-    if(Rast3d_get_vertical_unit(map) && G_strncasecmp(Rast3d_get_vertical_unit(map), "units", 5) != 0) {
-        /* Time handling */
-        if(G_is_units_type_temporal(Rast3d_get_vertical_unit2(map))) {
-            is_time = 1;
-            char long_name[1024];
-            char time_unit[1024];
-            
-            G_snprintf(long_name, 1024, "Time in %s", Rast3d_get_vertical_unit(map));
-            
-            if ((retval = nc_def_dim(ncid, TIME_NAME, region->depths, &time_dimid))) ERR(retval);
-            if ((retval = nc_def_var(ncid, TIME_NAME, NC_INT, 1, &time_dimid, &time_varid))) ERR(retval);
-            
-            /* Temporal unit */
-            if(G_has_raster3d_timestamp(map->fileName, map->mapset)) {
-                struct TimeStamp ts;
-                G_read_raster3d_timestamp(map->fileName, map->mapset, &ts);
-                G_snprintf(time_unit, 1024, "%s since %d-%02d-%02d %02d:%02d:%02.0f", 
-                           Rast3d_get_vertical_unit(map), ts.dt[0].year, ts.dt[0].month, ts.dt[0].day,
-                           ts.dt[0].hour,ts.dt[0].minute,ts.dt[0].second);
-            } else {
-                G_snprintf(time_unit, 1024, "%s since %s", Rast3d_get_vertical_unit(map), "1900-01-01 00:00:00");
-            }
-            
-            if ((retval = nc_put_att_text(ncid, time_varid, UNITS, strlen(time_unit), time_unit))) ERR(retval);
-            if ((retval = nc_put_att_text(ncid, time_varid, LONG_NAME, strlen(long_name), long_name))) ERR(retval)
-            if ((retval = nc_put_att_text(ncid, time_varid, "calendar", strlen("gregorian"), "gregorian"))) ERR(retval);
-            
-        } else {
-            if ((retval = nc_def_dim(ncid, Z_NAME, region->depths, &time_dimid))) ERR(retval);;
-            if ((retval = nc_def_var(ncid, Z_NAME, NC_FLOAT, 1, &time_dimid, &time_varid))) ERR(retval);
-            if ((retval = nc_put_att_text(ncid, time_varid, UNITS, strlen(Rast3d_get_vertical_unit(map)), Rast3d_get_vertical_unit(map)))) ERR(retval);
-            if ((retval = nc_put_att_text(ncid, time_varid, LONG_NAME, strlen(Z_LONG_NAME), Z_LONG_NAME))) ERR(retval);
-        }
-    }else {
-        if ((retval = nc_def_dim(ncid, Z_NAME, region->depths, &time_dimid))) ERR(retval);
-        if ((retval = nc_def_var(ncid, Z_NAME, NC_FLOAT, 1, &time_dimid, &time_varid))) ERR(retval);
-        if ((retval = nc_put_att_text(ncid, time_varid, UNITS, strlen("meter"), "meter"))) ERR(retval);
-        if ((retval = nc_put_att_text(ncid, time_varid, LONG_NAME, strlen(Z_LONG_NAME), Z_LONG_NAME))) ERR(retval);
+    is_absolute_time = 0;
+
+    if (Rast3d_get_vertical_unit(map) &&
+	G_strncasecmp(Rast3d_get_vertical_unit(map), "units", 5) != 0) {
+	/* Time handling */
+	if (G_is_units_type_temporal(Rast3d_get_vertical_unit2(map))) {
+	    is_time = 1;
+	    char long_name[1024];
+	    char time_unit[1024];
+
+	    G_snprintf(long_name, 1024, "Time in %s",
+		       Rast3d_get_vertical_unit(map));
+
+	    if ((retval =
+		 nc_def_dim(ncid, TIME_NAME, region->depths, &time_dimid)))
+		ERR(retval);
+	    if ((retval =
+		 nc_def_var(ncid, TIME_NAME, NC_INT, 1, &time_dimid,
+			    &time_varid)))
+		ERR(retval);
+
+	    /* Temporal unit */
+	    if (G_has_raster3d_timestamp(map->fileName, map->mapset)) {
+		struct TimeStamp ts;
+
+		G_read_raster3d_timestamp(map->fileName, map->mapset, &ts);
+
+                /* Days since datum in ISO norm*/
+		if (datetime_is_absolute(&ts.dt[0])) {
+		    is_absolute_time = 1;
+		    G_snprintf(time_unit, 1024,
+			       "%s since %d-%02d-%02d %02d:%02d:%02.0f",
+			       Rast3d_get_vertical_unit(map), ts.dt[0].year,
+			       ts.dt[0].month, ts.dt[0].day, ts.dt[0].hour,
+			       ts.dt[0].minute, ts.dt[0].second);
+		}
+		else {
+		    G_snprintf(time_unit, 1024, "%s",
+			       Rast3d_get_vertical_unit(map));
+		}
+
+	    }
+	    else {
+		G_snprintf(time_unit, 1024, "%s since %s",
+			   Rast3d_get_vertical_unit(map),
+			   "1900-01-01 00:00:00");
+	    }
+
+	    if ((retval =
+		 nc_put_att_text(ncid, time_varid, UNITS, strlen(time_unit),
+				 time_unit)))
+		ERR(retval);
+	    if ((retval =
+		 nc_put_att_text(ncid, time_varid, LONG_NAME,
+				 strlen(long_name), long_name)))
+		ERR(retval)
+		    if (is_absolute_time) {
+		    if ((retval =
+			 nc_put_att_text(ncid, time_varid, "calendar",
+					 strlen("gregorian"), "gregorian")))
+			ERR(retval);
+		}
+		else {
+		    if ((retval =
+			 nc_put_att_text(ncid, time_varid, "calendar",
+					 strlen("none"), "none")))
+			ERR(retval);
+		}
+	}
+	else {
+	    if ((retval =
+		 nc_def_dim(ncid, Z_NAME, region->depths, &time_dimid)))
+		ERR(retval);;
+	    if ((retval =
+		 nc_def_var(ncid, Z_NAME, NC_FLOAT, 1, &time_dimid,
+			    &time_varid)))
+		ERR(retval);
+	    if ((retval =
+		 nc_put_att_text(ncid, time_varid, UNITS,
+				 strlen(Rast3d_get_vertical_unit(map)),
+				 Rast3d_get_vertical_unit(map))))
+		ERR(retval);
+	    if ((retval =
+		 nc_put_att_text(ncid, time_varid, LONG_NAME,
+				 strlen(Z_LONG_NAME), Z_LONG_NAME)))
+		ERR(retval);
+	    if ((retval =
+		 nc_put_att_text(ncid, time_varid, STANDARD_NAME,
+				 strlen(Z_STANDARD_NAME), Z_STANDARD_NAME)))
+		ERR(retval);
+	}
+    }
+    else {
+	/* Default is z unit meter */
+	if ((retval = nc_def_dim(ncid, Z_NAME, region->depths, &time_dimid)))
+	    ERR(retval);
+	if ((retval =
+	     nc_def_var(ncid, Z_NAME, NC_FLOAT, 1, &time_dimid, &time_varid)))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, time_varid, UNITS, strlen("meter"),
+			     "meter")))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, time_varid, LONG_NAME, strlen(Z_LONG_NAME),
+			     Z_LONG_NAME)))
+	    ERR(retval);
+	if ((retval =
+	     nc_put_att_text(ncid, time_varid, STANDARD_NAME,
+			     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"))) ERR(retval);
     
+    /* z - axis orientation */
+    if ((retval =
+	 nc_put_att_text(ncid, time_varid, "positive", strlen("up"), "up")))
+	ERR(retval);
+
     /* Axis identifier attributes */
-    if ((retval = nc_put_att_text(ncid, lon_varid, "axis", 1, "X"))) ERR(retval);
-    if ((retval = nc_put_att_text(ncid, lat_varid, "axis", 1, "Y"))) ERR(retval);
-    if(is_time) {
-        if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "T"))) ERR(retval);
-    } else {
-        if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "Z"))) ERR(retval);
+    if ((retval = nc_put_att_text(ncid, lon_varid, "axis", 1, "X")))
+	ERR(retval);
+    if ((retval = nc_put_att_text(ncid, lat_varid, "axis", 1, "Y")))
+	ERR(retval);
+    if (is_time) {
+	if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "T")))
+	    ERR(retval);
+    }
+    else {
+	if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "Z")))
+	    ERR(retval);
     }
 
     dimids[0] = time_dimid;
     dimids[1] = lat_dimid;
     dimids[2] = lon_dimid;
-    
-    if(typeIntern == FCELL_TYPE) {
-        if ((retval = nc_def_var(ncid, param.input->answer, NC_FLOAT, NDIMS, dimids, varid))) ERR(retval); 
-    } else {
-        if ((retval = nc_def_var(ncid, param.input->answer, NC_DOUBLE, NDIMS, dimids, varid))) ERR(retval); 
+
+    if (typeIntern == FCELL_TYPE) {
+	if ((retval =
+	     nc_def_var(ncid, param.input->answer, NC_FLOAT, NDIMS, dimids,
+			varid)))
+	    ERR(retval);
     }
-    if(window.proj != PROJECTION_XY) {
-        if ((retval = nc_put_att_text(ncid, *varid, "grid_mapping", strlen("crs"), "crs"))) ERR(retval);
+    else {
+	if ((retval =
+	     nc_def_var(ncid, param.input->answer, NC_DOUBLE, NDIMS, dimids,
+			varid)))
+	    ERR(retval);
     }
-    
-   /* End define mode. */
-   if ((retval = nc_enddef(ncid))) ERR(retval);  
-    
+    if (window.proj != PROJECTION_XY) {
+	if ((retval =
+	     nc_put_att_text(ncid, *varid, "grid_mapping", strlen("crs"),
+			     "crs")))
+	    ERR(retval);
+    }
+
+    /* End define mode. */
+    if ((retval = nc_enddef(ncid)))
+	ERR(retval);
+
     /* 
      * Build coordinates, we need to use the cell center in case of spatial dimensions 
      * */
 
-    for(i = 0; i < region->cols; i++) {
-        c = region->west + i*region->ew_res + 0.5*region->ew_res;
-        nc_put_var1_float(ncid, lon_varid, &i, &c);     
+    for (i = 0; i < region->cols; i++) {
+	c = region->west + i * region->ew_res + 0.5 * region->ew_res;
+	nc_put_var1_float(ncid, lon_varid, &i, &c);
     }
-    
-    for(i = 0; i < region->rows; i++) {
-        c = region->south + i*region->ns_res + 0.5*region->ns_res;
-        nc_put_var1_float(ncid, lat_varid, &i, &c);     
+
+    for (i = 0; i < region->rows; i++) {
+	//c = region->south + i * region->ns_res + 0.5 * region->ns_res;
+	c = region->north - i * region->ns_res - 0.5 * region->ns_res;
+	nc_put_var1_float(ncid, lat_varid, &i, &c);
     }
-    
-    for(i = 0; i < region->depths; i++) {
-        if(is_time) {
-            c = region->bottom + i*region->tb_res;
-            time = (int)c;
-            nc_put_var1_int(ncid, time_varid, &i, &time);  
-	} else {
-            c = region->bottom + i*region->tb_res + 0.5*region->tb_res;
-            nc_put_var1_float(ncid, time_varid, &i, &c);     
+
+    for (i = 0; i < region->depths; i++) {
+	if (is_time) {
+	    c = region->bottom + i * region->tb_res;
+	    time = (int)c;
+	    nc_put_var1_int(ncid, time_varid, &i, &time);
+	}
+	else {
+	    c = region->bottom + i * region->tb_res + 0.5 * region->tb_res;
+	    nc_put_var1_float(ncid, time_varid, &i, &c);
 	}
     }
 }
 
 /*---------------------------------------------------------------------------*/
 
-static void write_netcdf_data(int ncid, RASTER3D_Region *region, int varid)
+static void write_netcdf_data(int ncid, RASTER3D_Region * region, int varid)
 {
     DCELL dvalue;
     FCELL fvalue;
@@ -293,29 +474,30 @@ static void write_netcdf_data(int ncid, RASTER3D_Region *region, int varid)
     typeIntern = Rast3d_tile_type_map(map);
 
     for (z = 0; z < depths; z++) {
-        G_percent(z, depths, 1);
-        for (y = rows - 1; y >= 0; y--) { /* rows count from south to north */
-            for (x = 0; x < cols; x++) {
-                coords[0] = z;
-                coords[1] = y;
-                coords[2] = x;
-                 if (typeIntern == FCELL_TYPE) {
-                    
-                    Rast3d_get_value(map, x, y, z, &fvalue, FCELL_TYPE);
-                    
-                    if (!Rast3d_is_null_value_num(&fvalue, FCELL_TYPE)) {
-                        nc_put_var1_float(ncid, varid, coords, &fvalue);
-                    }
-                } else {
-                    
-                    Rast3d_get_value(map, x, y, z, &dvalue, DCELL_TYPE);
-                    
-                    if (!Rast3d_is_null_value_num(&dvalue, DCELL_TYPE)) {
-                        nc_put_var1_double(ncid, varid, coords, &dvalue);
-                    }
-                }
-            }
-        }
+	G_percent(z, depths, 1);
+	for (y = 0; y < rows; y++) {
+	    for (x = 0; x < cols; x++) {
+		coords[0] = z;
+		coords[1] = y;
+		coords[2] = x;
+		if (typeIntern == FCELL_TYPE) {
+
+		    Rast3d_get_value(map, x, y, z, &fvalue, FCELL_TYPE);
+
+		    if (!Rast3d_is_null_value_num(&fvalue, FCELL_TYPE)) {
+			nc_put_var1_float(ncid, varid, coords, &fvalue);
+		    }
+		}
+		else {
+
+		    Rast3d_get_value(map, x, y, z, &dvalue, DCELL_TYPE);
+
+		    if (!Rast3d_is_null_value_num(&dvalue, DCELL_TYPE)) {
+			nc_put_var1_double(ncid, varid, coords, &dvalue);
+		    }
+		}
+	    }
+	}
     }
     G_percent(1, 1, 1);
     G_percent_reset();
@@ -326,11 +508,11 @@ static void write_netcdf_data(int ncid, RASTER3D_Region *region, int varid)
 int main(int argc, char *argv[])
 {
     RASTER3D_Region region;
-    
+
     int ncid;
     int retval;
     int varid;
-    
+
     int changemask = 0;
     struct GModule *module;
 
@@ -340,18 +522,18 @@ int main(int argc, char *argv[])
     G_add_keyword(_("raster3d"));
     G_add_keyword(_("netcdf"));
     G_add_keyword(_("export"));
-    module->description =
-        _("Export a 3D raster map as netcdf file.");
+    module->description = _("Export a 3D raster map as netcdf file.");
 
     /* Get parameters from user */
     setParams();
 
     /* Have GRASS get inputs */
     if (G_parser(argc, argv))
-        exit(EXIT_FAILURE);
+	exit(EXIT_FAILURE);
 
     if (NULL == G_find_raster3d(param.input->answer, ""))
-        Rast3d_fatal_error(_("3D raster map <%s> not found"), param.input->answer);
+	Rast3d_fatal_error(_("3D raster map <%s> not found"),
+			   param.input->answer);
 
     /* Initiate the default settings */
     Rast3d_init_defaults();
@@ -360,43 +542,49 @@ int main(int argc, char *argv[])
     Rast3d_get_window(&region);
 
     /* Open the map and use XY cache mode */
-    map = Rast3d_open_cell_old(param.input->answer, G_find_raster3d(param.input->answer, ""), &region,
-                          RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
+    map =
+	Rast3d_open_cell_old(param.input->answer,
+			     G_find_raster3d(param.input->answer, ""),
+			     &region, RASTER3D_TILE_SAME_AS_FILE,
+			     RASTER3D_USE_CACHE_DEFAULT);
 
     if (map == NULL)
-        G_fatal_error(_("Error opening 3d raster map <%s>"), param.input->answer);
+	G_fatal_error(_("Error opening 3d raster map <%s>"),
+		      param.input->answer);
 
     /* Create netcdf file */
-    if ((retval = nc_create(param.output->answer, NC_CLOBBER, &ncid))) ERR(retval);
-    
+    if ((retval = nc_create(param.output->answer, NC_CLOBBER, &ncid)))
+	ERR(retval);
+
     write_netcdf_header(ncid, &region, &varid);
 
     /*if requested set the Mask on */
     if (param.mask->answer) {
-        if (Rast3d_mask_file_exists()) {
-            changemask = 0;
-            if (Rast3d_mask_is_off(map)) {
-                Rast3d_mask_on(map);
-                changemask = 1;
-            }
-        }
+	if (Rast3d_mask_file_exists()) {
+	    changemask = 0;
+	    if (Rast3d_mask_is_off(map)) {
+		Rast3d_mask_on(map);
+		changemask = 1;
+	    }
+	}
     }
-    
+
     write_netcdf_data(ncid, &region, varid);
 
     /*We set the Mask off, if it was off bevor */
     if (param.mask->answer) {
-        if (Rast3d_mask_file_exists())
-            if (Rast3d_mask_is_on(map) && changemask)
-                Rast3d_mask_off(map);
+	if (Rast3d_mask_file_exists())
+	    if (Rast3d_mask_is_on(map) && changemask)
+		Rast3d_mask_off(map);
     }
 
     /* Close files and exit */
     if (!Rast3d_close(map))
-        fatalError(_("Unable to close 3D raster map"));
-    
+	fatalError(_("Unable to close 3D raster map"));
+
     /* Close the netcdf file */
-    if ((retval = nc_close(ncid))) ERR(retval);
-    
+    if ((retval = nc_close(ncid)))
+	ERR(retval);
+
     return 0;
 }

raster3d/r3.out.netcdf/test_suite/test.r3.out.netcdf.sh

@@ -6,19 +6,28 @@
 # @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=-40 n=40 w=-20 e=100 b=-25 t=25 res=10 res3=10 -p3
+g.region s=-90 n=90 w=-180 e=180 b=0 t=50 res=10 res3=10 -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.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"
 # @test
 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.info volume_double
+#ncdump -h test_double.nc
 r3.out.netcdf --o 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.info volume_time_float
+#ncdump -h test_time_float.nc
 
-ncdump -h test_float.nc
-ncdump -h test_double.nc
-ncdump -h test_time_double.nc