Merged branch develop into master

DESCRIPTION

@@ -3,6 +3,7 @@ Type: Package
 Title: Pre- and post-processing of MODFLOW files
 Version: 0.4.0
 Date: 2016-06-27
+Authors@R: person("Bart", "Rogiers", role = c("aut", "cre"), email = "brogiers@sckcen.be")
 Author: Bart Rogiers <brogiers@sckcen.be>
 Maintainer: Bart Rogiers <brogiers@sckcen.be>
 Description: The RMODFLOW package provides a set of tools for groundwater flow
@@ -22,8 +23,10 @@ Imports:
     sp,
     rgdal,
     animation,
-    tools
+    tools,
+    lubridate
 URL: https://github.com/rogiersbart/RMODFLOW
+BugReports: https://github.com/rogiersbart/RMODFLOW/issues
 RoxygenNote: 5.0.1
 Suggests: knitr,
     rmarkdown

NAMESPACE

@@ -24,8 +24,11 @@ export(cell_info)
 export(convert_bud_to_darcy)
 export(convert_dis_to_saturated_dis)
 export(convert_grid_to_xyz)
+export(convert_hob_to_time_series)
 export(convert_huf_to_dis)
 export(convert_huf_to_grid)
+export(convert_huf_to_mask)
+export(convert_huf_to_nlay)
 export(convert_ibound_to_neighbours)
 export(convert_id_to_ijk)
 export(convert_ijk_to_id)

R/cell_coordinates.dis.R

@@ -20,12 +20,14 @@ cell_coordinates.dis <- function(dis,
   cell_coordinates$y[,,] <- rev(cumsum(rev(dis$delc))-rev(dis$delc)/2)
   cell_coordinates$x[,,] <- rep(c(cumsum(dis$delr)-dis$delr/2),each=dis$nrow)
   if(include_faces) {
+    dis$delr <- array(rep(dis$delr,each=dis$nrow),dim=c(dis$nrow,dis$ncol,dis$nlay))
+    dis$delc <- array(rep(dis$delc,dis$ncol),dim=c(dis$nrow,dis$ncol,dis$nlay))
     cell_coordinates$lower <- dis$botm
     cell_coordinates$upper <- 2 * cell_coordinates$z - dis$botm
     cell_coordinates$left <- cell_coordinates$x - dis$delr/2
     cell_coordinates$right <- cell_coordinates$x + dis$delr/2 
-    cell_coordinates$front <- cell_coordinates$y - rev(dis$delc)/2
-    cell_coordinates$back <- cell_coordinates$y + rev(dis$delc)/2
+    cell_coordinates$front <- cell_coordinates$y - dis$delc/2
+    cell_coordinates$back <- cell_coordinates$y + dis$delc/2
   }
   return(cell_coordinates)
 }

R/cell_coordinates.huf.R

@@ -1,13 +1,16 @@
 #' Get cell coordinates from a huf object
 #' 
 #' @param huf huf object
+#' @param dis dis object, corresponding to the huf object
+#' @param include_faces logical; should face coordinates be included?
 #' @return 3d array with cell coordinates
 #'
 #' @rdname cell_coordinates
 #' @method cell_coordinates huf
 #' @export
-cell_coordinates.huf <- function(huf, dis = NULL)
-{
+cell_coordinates.huf <- function(huf,
+                                 dis = NULL,
+                                 include_faces = FALSE) {
   cell_coordinates <- NULL
   cell_coordinates$z <- huf$top - huf$thck/2
   class(cell_coordinates$z) <- 'rmodflow_3d_array'
@@ -17,5 +20,15 @@ cell_coordinates.huf <- function(huf, dis = NULL)
     cell_coordinates$y[,,] <- rev(cumsum(rev(dis$delc))-rev(dis$delc)/2)
     cell_coordinates$x[,,] <- rep(c(cumsum(dis$delr)-dis$delr/2),each=dis$nrow)
   }
+  if(include_faces) {
+    dis$delr <- array(rep(dis$delr,each=dis$nrow),dim=c(dis$nrow,dis$ncol,huf$nhuf))
+    dis$delc <- array(rep(dis$delc,dis$ncol),dim=c(dis$nrow,dis$ncol,huf$nhuf))
+    cell_coordinates$lower <- cell_coordinates$z - huf$thck/2
+    cell_coordinates$upper <- cell_coordinates$z + huf$thck/2
+    cell_coordinates$left <- cell_coordinates$x - dis$delr/2
+    cell_coordinates$right <- cell_coordinates$x + dis$delr/2 
+    cell_coordinates$front <- cell_coordinates$y - dis$delc/2
+    cell_coordinates$back <- cell_coordinates$y + dis$delc/2
+  }
   return(cell_coordinates)
 }

R/convert-hob-to-time-series.R

@@ -0,0 +1,14 @@
+#' Convert a hob object to a time series data frame
+#' 
+#' @param hob hob object
+#' @param dis dis object
+#' @param prj prj object
+#' @return time series data frame containing name, time and head columns
+#' @export
+convert_hob_to_time_series <- function(hob,
+                                       dis,
+                                       prj) {
+  toffset <- lubridate::days(ifelse(hob$irefsp==1,0,cumsum(dis$perlen)[hob$irefsp-1]) + hob$toffset * hob$tomulth)
+  time_series <- data.frame(name = hob$obsloc, time = prj$starttime + toffset, head = hob$hobs)
+  return(time_series)
+}

R/convert-huf-to-mask.R

@@ -0,0 +1,14 @@
+#' Convert a huf to a mask object
+#' 
+#' @param huf huf object
+#' @param dis dis object, corresponding to the huf object
+#' @param bas bas object, corresponding to the huf object
+#' @return mask rmodflow_3d_array
+#' @export
+convert_huf_to_mask <- function(huf, dis, bas) {
+  mask <- RMODFLOW::convert_huf_to_nlay(huf = huf, dis = dis, bas = bas)
+  mask[which(mask==0)] <- NA
+  mask <- mask/mask
+  mask[which(huf$thck==0)] <- NA
+  return(mask)
+}

R/convert-huf-to-nlay.R

@@ -0,0 +1,19 @@
+#' Convert a huf object to an rmodflow_3d_array with the number of numerical layers per hydrogeological unit
+#' 
+#' @param huf huf object
+#' @param dis dis object, corresponding to the huf object
+#' @param bas bas object, corresponding to the huf object
+#' @return nlay rmodflow_3d_array
+#' @export
+convert_huf_to_nlay <- function(huf, dis, bas) {
+  nlay <- huf$top * 0
+  huf_coordinates <- RMODFLOW::cell_coordinates(huf, dis = dis, include_faces = TRUE)
+  dis_coordinates <- RMODFLOW::cell_coordinates(dis, include_faces = TRUE)
+  ibound <- abs(bas$ibound)
+  for(i in 1:huf$nhuf) {
+    for(j in 1:dis$nlay) {
+      nlay[,,i] <- nlay[,,i] + (!(dis_coordinates$upper[,,j] < huf_coordinates$lower[,,i] | dis_coordinates$lower[,,j] > huf_coordinates$upper[,,i])) * ibound[,,j]
+    }
+  }
+  return(nlay)
+}

R/convert_bud_to_darcy.R

@@ -10,13 +10,16 @@ convert_bud_to_darcy <- function(bud, dis, hed = NULL) {
   delc <- create_rmodflow_array(rep(dis$delc,dis$ncol),dim=dim(bud$flow_right_face))
   delr <- create_rmodflow_array(rep(dis$delr,each=dis$nrow),dim=dim(bud$flow_right_face))
   darcy <- list()
+  ## temp fix for Cas' files
+  ## TODO: fix this when either of the components is missing
+    if(!"flow_front_face" %in% names(bud)) bud$flow_front_face <- bud$flow_right_face * 0
   darcy$right <- bud$flow_right_face
   darcy$front <- -bud$flow_front_face
   darcy$lower <- -bud$flow_lower_face/delc/delr
   darcy$left <- darcy$back <- darcy$upper <- darcy$right * 0
-  darcy$left[,c(2:dis$ncol),,] <- darcy$right[,c(1:(dis$ncol-1)),,]
-  darcy$back[c(2:dis$nrow),,,] <- darcy$front[c(1:(dis$nrow-1)),,,]
-  darcy$upper[,,c(2:dis$nlay),] <- darcy$lower[,,c(1:(dis$nlay-1)),]
+  if(dis$ncol > 1) darcy$left[,c(2:dis$ncol),,] <- darcy$right[,c(1:(dis$ncol-1)),,] else darcy$left <- darcy$right * 0
+  if(dis$nrow > 1) darcy$back[c(2:dis$nrow),,,] <- darcy$front[c(1:(dis$nrow-1)),,,] else darcy$back <- darcy$front * 0
+  if(dis$nlay > 1) darcy$upper[,,c(2:dis$nlay),] <- darcy$lower[,,c(1:(dis$nlay-1)),] else darcy$upper <- darcy$lower * 0
   darcy$right <- darcy$right/delc/thck
   darcy$left <- darcy$left/delc/thck
   darcy$front <- darcy$front/delr/thck

R/convert_xyz_to_grid.R

@@ -2,14 +2,18 @@
 #' 
 #' @param x real world x coordinate
 #' @param y real world y coordinate
-#' @param prj prj object
 #' @param z real world z coordinate; optional
+#' @param prj prj object
 #' @param dis dis object; optional
+#' @param output character; containing 'xyz','ijk' and/or 'off' for the return of x, y, z, i, j, k, roff, coff and loff modflow coordinates
 #' @details
-#' If dis is not provided, only x, y and z coordinates are returned. If z is not provided, no third dimension coordinates are returned.
-#' @return data frame with x, y, z, i, j, k, roff, coff and loff modflow coordinates
+#' If dis is not provided, only x, y and z coordinates are returned. If z is not provided, no third dimension coordinates are returned. For the x, y and z modflow coordinates, the origin is placed at the lower left corner of the grid.
+#' @return data frame with modflow coordinates
 #' @export
-convert_xyz_to_grid <- function(x,y,prj=NULL,z=NULL,dis=NULL) {
+convert_xyz_to_grid <- function(x,y,prj=NULL,z=NULL,dis=NULL,output='xyz') {
+  output_xyz <- grepl('xyz',output)
+  output_ijk <- grepl('ijk',output)
+  output_off <- grepl('off',output)
   if(!is.null(prj)) {
     x <- x-prj$origin[1]
     y <- y-prj$origin[2]
@@ -22,15 +26,14 @@ convert_xyz_to_grid <- function(x,y,prj=NULL,z=NULL,dis=NULL) {
   }
   dat <- data.frame(x=x,y=y)
   if(!is.null(z)) dat$z <- z
-  if(!is.null(dis)) {
+  if(output_ijk | output_off) {
+    if(is.null(dis)) error('Please provide dis argument ...')    
     if(ncol(dat)==3) {
       dis$thck <- dis$tops <- dis$botm
       dis$thck[,,1] <- dis$top - dis$botm[,,1]
       dis$tops[,,1] <- dis$top
-      for(k in 2:dis$nlay) {
-        dis$thck[,,k] <- dis$botm[,,(k-1)] - dis$botm[,,k]
-        dis$tops[,,k] <- dis$botm[,,(k-1)]
-      }
+      dis$thck[,,2:dis$nlay] <- dis$botm[,,(2:dis$nlay-1)] - dis$botm[,,2:dis$nlay]
+      dis$tops[,,2:dis$nlay] <- dis$botm[,,(2:dis$nlay-1)]
     }
     for(i in 1:nrow(dat)) {
       dat$i[i] <- which(cumsum(dis$delc) > sum(dis$delc)-dat$y[i])[1]
@@ -42,6 +45,13 @@ convert_xyz_to_grid <- function(x,y,prj=NULL,z=NULL,dis=NULL) {
         dat$loff[i] <- -(dat$z[i]-(dis$tops[dat$i[i],dat$j[i],dat$k[i]]+dis$botm[dat$i[i],dat$j[i],dat$k[i]])/2)/dis$thck[dat$i[i],dat$j[i],dat$k[i]]
       }
     }
+    if(output_xyz & output_ijk & output_off) return(dat)
+    if(output_xyz & output_ijk & !output_off) return(dat[,c('x','y','z','i','j','k')])
+    if(output_xyz & !output_ijk & output_off) return(dat[,c('x','y','z','roff','coff','loff')])
+    if(!output_xyz & output_ijk & output_off) return(dat[,c('i','j','k','roff','coff','loff')])
+    if(!output_xyz & !output_ijk & output_off) return(dat[,c('roff','coff','loff')])
+    if(!output_xyz & output_ijk & !output_off) return(dat[,c('i','j','k')])
+  } else {
+    return(dat)
   }
-  return(dat)
 }

R/plot.rmodflow_2d_array.R

@@ -47,16 +47,22 @@ plot.rmodflow_2d_array <- function(array,
                           height=NULL,
                           title = NULL) {
   if(plot3d) {
-    x <- (cumsum(dis$delr)-dis$delr/2)
-    y <- sum(dis$delc) - (cumsum(dis$delc)-dis$delc/2)
+    xyz <- cell_coordinates(dis)
+    x <- xyz$x[,,1]
+    y <- xyz$y[,,1]
+    if(!is.null(prj)) {
+      xyz <- convert_grid_to_xyz(x=c(x),y=c(y),prj=prj)
+      x[,] <- xyz$x
+      y[,] <- xyz$y
+    }
     z <- t(height)*height_exaggeration
     if(!add) rgl::open3d()
     colorlut <- colorRampPalette(colour_palette(nlevels))(25) # height color lookup table
     col <- colorlut[ round(approx(seq(zlim[1],zlim[2],length=25+1),seq(0.5,25+0.5,length=25+1),xout=c(t(array)),rule=2)$y) ] # assign colors to heights for each point
     alpha <- rep(1,length(col))
     alpha[which(c(t(mask))==0)] <- 0
-    if(type=='fill') rgl::surface3d(x,y,z,color=col,alpha=alpha,back='lines',smooth=FALSE) 
-    if(type=='grid') rgl::surface3d(x,y,z,front='lines',alpha=alpha,back='lines',smooth=FALSE) 
+    if(type=='fill') rgl::surface3d(t(x),t(y),z,color=col,alpha=alpha,back='lines',smooth=FALSE) 
+    if(type=='grid') rgl::surface3d(t(x),t(y),z,front='lines',alpha=alpha,back='lines',smooth=FALSE) 
   } else {
     xy <- expand.grid(cumsum(dis$delr)-dis$delr/2,sum(dis$delc)-(cumsum(dis$delc)-dis$delc/2))
     names(xy) <- c('x','y')
@@ -116,6 +122,14 @@ plot.rmodflow_2d_array <- function(array,
                  coord_equal() + ggtitle(title))
       }
     } else if(type=='contour') {
+      if(!is.null(prj)) {
+        new_xy <- convert_grid_to_xyz(x=xy$x,y=xy$y,prj=prj)
+        xy$x <- new_xy$x
+        xy$y <- new_xy$y
+      }
+      if(!is.null(crs)) {
+        xy <- RMODFLOW:::convert_coordinates(xy,from=CRS(prj$projection),to=crs)
+      }
       xy$z <- c(t(array*mask^2))
       xyBackup <- xy
       xy <- na.omit(xy)

R/read_bud.R

@@ -8,7 +8,9 @@
 #' @importFrom readr read_lines
 #' @export
 read_bud <- function(file = {cat('Please select bud file ...\n'); file.choose()},
-                     binary = TRUE) {
+                     binary = TRUE,
+                     double_precision = FALSE) {
+  if(double_precision) nbytes <- 8 else nbytes <- 4
   if(binary) {
     con <- file(file,open='rb')
     bud <- list()
@@ -38,12 +40,12 @@ read_bud <- function(file = {cat('Please select bud file ...\n'); file.choose()}
       bud[[name]][[kper]][[kstp]]$nlay <- readBin(con,what='integer',n=1)
 
       if(bud[[name]][[kper]][[kstp]]$nlay > 0) {
-        bud[[name]][[kper]][[kstp]]$data <- aperm(array(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow*bud[[name]][[kper]][[kstp]]$nlay,size=4),dim=c(bud[[name]][[kper]][[kstp]]$ncol,bud[[name]][[kper]][[kstp]]$nrow,bud[[name]][[kper]][[kstp]]$nlay)),c(2,1,3))
+        bud[[name]][[kper]][[kstp]]$data <- aperm(array(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow*bud[[name]][[kper]][[kstp]]$nlay,size = nbytes),dim=c(bud[[name]][[kper]][[kstp]]$ncol,bud[[name]][[kper]][[kstp]]$nrow,bud[[name]][[kper]][[kstp]]$nlay)),c(2,1,3))
       } else {
         bud[[name]][[kper]][[kstp]]$itype <- readBin(con,what='integer',n=1)
-        bud[[name]][[kper]][[kstp]]$delt <- readBin(con,what='numeric',n=1,size=4)
-        bud[[name]][[kper]][[kstp]]$pertim <- readBin(con,what='numeric',n=1,size=4)
-        bud[[name]][[kper]][[kstp]]$totim <- readBin(con,what='numeric',n=1,size=4)
+        bud[[name]][[kper]][[kstp]]$delt <- readBin(con,what='numeric',n=1,size = nbytes)
+        bud[[name]][[kper]][[kstp]]$pertim <- readBin(con,what='numeric',n=1,size = nbytes)
+        bud[[name]][[kper]][[kstp]]$totim <- readBin(con,what='numeric',n=1,size = nbytes)
         if(bud[[name]][[kper]][[kstp]]$itype==5) {
           bud[[name]][[kper]][[kstp]]$nval <- readBin(con,what='integer',n=1)
         } else {
@@ -64,49 +66,55 @@ read_bud <- function(file = {cat('Please select bud file ...\n'); file.choose()}
             names(bud[[name]][[kper]][[kstp]]$data)[2] <- 'value'
             # add reading ctmps here!
             for(nr in 1:bud[[name]][[kper]][[kstp]]$nlist) {
-              bud[[name]][[kper]][[kstp]]$data[nr,] <- c(readBin(con,what='integer',n=1),readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$nval,size=4))
+              bud[[name]][[kper]][[kstp]]$data[nr,] <- c(readBin(con,what='integer',n=1),readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$nval,size = nbytes))
             }
           }
         }
         if(bud[[name]][[kper]][[kstp]]$itype %in% c(0,1)) {
-          bud[[name]][[kper]][[kstp]]$data <- aperm(array(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow*abs(bud[[name]][[kper]][[kstp]]$nlay),size=4),dim=c(bud[[name]][[kper]][[kstp]]$ncol,bud[[name]][[kper]][[kstp]]$nrow,abs(bud[[name]][[kper]][[kstp]]$nlay))),c(2,1,3))
+          bud[[name]][[kper]][[kstp]]$data <- aperm(array(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow*abs(bud[[name]][[kper]][[kstp]]$nlay),size = nbytes),dim=c(bud[[name]][[kper]][[kstp]]$ncol,bud[[name]][[kper]][[kstp]]$nrow,abs(bud[[name]][[kper]][[kstp]]$nlay))),c(2,1,3))
           class(bud[[name]][[kper]][[kstp]]$data) <- 'rmodflow_3d_array'
         }
         if(bud[[name]][[kper]][[kstp]]$itype ==3) {
           bud[[name]][[kper]][[kstp]]$layer <- matrix(readBin(con,what='integer',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow),ncol=bud[[name]][[kper]][[kstp]]$ncol,nrow=bud[[name]][[kper]][[kstp]]$nrow,byrow=TRUE)
           class(bud[[name]][[kper]][[kstp]]$layer) <- 'rmodflow_2d_array'
-          bud[[name]][[kper]][[kstp]]$data <- matrix(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow,size=4),ncol=bud[[name]][[kper]][[kstp]]$ncol,nrow=bud[[name]][[kper]][[kstp]]$nrow,byrow=TRUE)
+          bud[[name]][[kper]][[kstp]]$data <- matrix(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow,size = nbytes),ncol=bud[[name]][[kper]][[kstp]]$ncol,nrow=bud[[name]][[kper]][[kstp]]$nrow,byrow=TRUE)
           class(bud[[name]][[kper]][[kstp]]$data) <- 'rmodflow_2d_array'
         }
         if(bud[[name]][[kper]][[kstp]]$itype ==4) {
-          bud[[name]][[kper]][[kstp]]$data <- matrix(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow,size=4),ncol=bud[[name]][[kper]][[kstp]]$ncol,nrow=bud[[name]][[kper]][[kstp]]$nrow,byrow=TRUE)
+          bud[[name]][[kper]][[kstp]]$data <- matrix(readBin(con,what='numeric',n=bud[[name]][[kper]][[kstp]]$ncol*bud[[name]][[kper]][[kstp]]$nrow,size = nbytes),ncol=bud[[name]][[kper]][[kstp]]$ncol,nrow=bud[[name]][[kper]][[kstp]]$nrow,byrow=TRUE)
           class(bud[[name]][[kper]][[kstp]]$data) <- 'rmodflow_2d_array'
         }
       }
 
       # set data as the main list item, and include all parameters as attributes
-        for(i in 1:(length(bud[[name]][[kper]][[kstp]])-1)) {
-          attr(bud[[name]][[kper]][[kstp]]$data,names(bud[[name]][[kper]][[kstp]])[i]) <- bud[[name]][[kper]][[kstp]][[i]]
+        if("data" %in% names(bud[[name]][[kper]][[kstp]])) {
+          for(i in 1:(length(bud[[name]][[kper]][[kstp]])-1)) {
+            attr(bud[[name]][[kper]][[kstp]]$data,names(bud[[name]][[kper]][[kstp]])[i]) <- bud[[name]][[kper]][[kstp]][[i]]
+          }
+          bud[[name]][[kper]][[kstp]] <- bud[[name]][[kper]][[kstp]]$data
         }
-        bud[[name]][[kper]][[kstp]] <- bud[[name]][[kper]][[kstp]]$data
-
+
       kstp <- readBin(con,what='integer',n=1)
       kper <- readBin(con,what='integer',n=1)
       desc <- readChar(con,nchars=16)
     }
 
     # create rmodflow_4d_array for list items with itype 0 or 1
       for (i in 1:length(bud)) {
-        if (attr(bud[[i]][[1]][[1]],'itype') %in% c(0,1)) {
-          bud_item <- bud[[i]]
-          bud[[i]] <- unlist(bud_item)
-          bud[[i]] <- create_rmodflow_array(bud[[i]], dim = c(attr(bud_item[[1]][[1]],'nrow'), attr(bud_item[[1]][[1]],'ncol'), abs(attr(bud_item[[1]][[1]],'nlay')), length(bud[[i]])/prod(attr(bud_item[[1]][[1]],'nrow'), attr(bud_item[[1]][[1]],'ncol'), abs(attr(bud_item[[1]][[1]],'nlay')))))
-          ats <- attributes(bud_item[[length(bud_item)]][[length(bud_item[[length(bud_item)]])]])
-          ats <- ats[-which(names(ats) == 'dim')]
-          for(at in 1:length(ats)) {
-            attr(bud[[i]], names(ats)[at]) <- ats[at][[1]]
-          } 
-          class(bud[[i]]) <- 'rmodflow_4d_array'
+        # problems when nlist is 0 and hence no data is present -> check for itype in names
+        # also for storage the first timestep seems to be empty -> not changed into 4d array at the moment
+        if (! "itype" %in% names(bud[[i]][[1]][[1]]) & !is.null(bud[[i]][[1]][[1]])) {
+          if (attr(bud[[i]][[1]][[1]],'itype') %in% c(0,1)) {
+            bud_item <- bud[[i]]
+            bud[[i]] <- unlist(bud_item)
+            bud[[i]] <- create_rmodflow_array(bud[[i]], dim = c(attr(bud_item[[1]][[1]],'nrow'), attr(bud_item[[1]][[1]],'ncol'), abs(attr(bud_item[[1]][[1]],'nlay')), length(bud[[i]])/prod(attr(bud_item[[1]][[1]],'nrow'), attr(bud_item[[1]][[1]],'ncol'), abs(attr(bud_item[[1]][[1]],'nlay')))))
+            ats <- attributes(bud_item[[length(bud_item)]][[length(bud_item[[length(bud_item)]])]])
+            ats <- ats[-which(names(ats) == 'dim')]
+            for(at in 1:length(ats)) {
+              attr(bud[[i]], names(ats)[at]) <- ats[at][[1]]
+            } 
+            class(bud[[i]]) <- 'rmodflow_4d_array'
+          }
         }
       }