as_tibble no longer working on raster brick

[afredston]

I'm re-running some old code that I haven't revisited since the update to R 4.0.X. Some calls to tabularaster::as_tibble() work fine, but I keep getting an error converting a raster brick into a df here that says:

Error: Assigned data `dimindex` must be compatible with existing data.
x Existing data has 6811182 rows.
x Assigned data has 20433546 rows.
ℹ Only vectors of size 1 are recycled.

I've tried updating packages (raster, tabularaster, sf, and rerddap), clearing the rerddap cache on my local machine, and testing different OISST raster bricks. The HadISST ones work fine, but none of the OISST ones do -- they just have slightly different dimension sizes in the error. Here's my session info:

> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rnaturalearth_0.1.0 lubridate_1.7.8     tabularaster_0.6.0  here_0.1           
 [5] forcats_0.5.0       stringr_1.4.0       dplyr_1.0.0         purrr_0.3.4        
 [9] readr_1.3.1         tidyr_1.0.3         tibble_3.1.0        ggplot2_3.3.2      
[13] tidyverse_1.3.0     sf_0.9-7            rerddap_0.7.4       raster_3.3-13      
[17] sp_1.4-4           

loaded via a namespace (and not attached):
 [1] httr_1.4.2            jsonlite_1.7.2        splines_4.0.2         modelr_0.1.8         
 [5] assertthat_0.2.1      triebeard_0.3.0       urltools_1.7.3        cellranger_1.1.0     
 [9] pillar_1.5.1          backports_1.1.10      lattice_0.20-41       glue_1.4.2           
[13] digest_0.6.27         polyclip_1.10-0       rvest_0.3.5           colorspace_1.4-1     
[17] Matrix_1.3-2          pkgconfig_2.0.3       httpcode_0.3.0        broom_0.5.6          
[21] haven_2.2.0           scales_1.1.1          tensor_1.5            spatstat.utils_1.17-0
[25] mgcv_1.8-31           generics_0.0.2        ellipsis_0.3.1        withr_2.4.1          
[29] cli_2.3.1             magrittr_2.0.1        crayon_1.4.1          readxl_1.3.1         
[33] deldir_0.1-29         fs_1.5.0              ncdf4_1.17            fansi_0.4.2          
[37] nlme_3.1-148          xml2_1.3.2            class_7.3-17          tools_4.0.2          
[41] data.table_1.12.8     hms_0.5.3             lifecycle_1.0.0       munsell_0.5.0        
[45] reprex_0.3.0          packrat_0.5.0         compiler_4.0.2        e1071_1.7-4          
[49] rlang_0.4.10          classInt_0.4-3        units_0.6-7           grid_4.0.2           
[53] rstudioapi_0.13       rappdirs_0.3.3        goftest_1.2-2         gtable_0.3.0         
[57] codetools_0.2-16      abind_1.4-5           DBI_1.1.0             curl_4.3             
[61] R6_2.5.0              rgdal_1.5-18          utf8_1.1.4            rprojroot_2.0.2      
[65] KernSmooth_2.23-17    hoardr_0.5.2          stringi_1.5.3         spatstat.data_1.4-3  
[69] spatstat_1.64-1       crul_0.9.0            Rcpp_1.0.6            vctrs_0.3.6          
[73] rpart_4.1-15          dbplyr_1.4.3          tidyselect_1.1.0     

[mdsumner]

I don't know how to reproduce this, I've tried a few test bricks including with daily oisst and can't get this error. Can you isolate it to a particular object?

[afredston]

I haven't but I will try to write a reprex, stand by!

[afredston]

I isolated the problem to the masking step. Even though (in reprex below) wc_oisst_brick1 and wc_oisst_crop1 are both raster bricks and basically the same size, as_tibble() works for the former but not the latter. (It also works on a hadISST cropped raster brick, oddly.)

library(rerddap)
library(raster)
library(sf)
library(tidyverse)
library(tabularaster)

select <- dplyr::select

# get OISST
wc_latrange <- c(34.4, 50)
wc_lonrange <- c(-126, -117)

oisstID <- "ncdcOisst2Agg_LonPM180"
oisst_time1 <- c("1982-01-01","1989-12-31") 
oisst_fields <- "sst"

wc_oisst_grid1 <- griddap(oisstID, time=oisst_time1, latitude = wc_latrange, longitude = wc_lonrange, fields=oisst_fields)

wc_oisst_nc_file1 <- wc_oisst_grid1$summary$filename
wc_oisst_brick1 <- brick(wc_oisst_nc_file1)

# get mask
bathy <- "etopo180" # https://coastwatch.pfeg.noaa.gov/erddap/griddap/etopo180.html 

wc.bathy.grid <- griddap(bathy, latitude=wc_latrange, longitude=wc_lonrange)

wc.topo.file <- wc.bathy.grid$summary$filename

wc.bathy.raster <- raster(wc.topo.file)

wc.bathy.crop <- wc.bathy.raster
wc.bathy.crop[wc.bathy.crop > 0] <- NA
wc.bathy.crop[wc.bathy.crop < -600] <- NA

wc.bathy.mask <- wc.bathy.crop %>% 
  as(., "SpatialPolygonsDataFrame") %>% 
  st_as_sf()

# mask OISST and convert to df
wc_oisst_crop1 <- mask(wc_oisst_brick1, as_Spatial(wc.bathy.mask))

# this produces an error, but not if you replace wc_oisst_crop1 with wc_oisst_brick1
wc_oisst_df1 <- tabularaster::as_tibble(wc_oisst_crop1, cell=FALSE, dim=TRUE, values=TRUE, xy=TRUE) %>% 
  filter(!is.na(cellvalue))%>%
  rename("sst" = cellvalue,
         "time" = dimindex)%>%
  mutate(time = as_date(time),
         year = year(time),
         month=month(time))

[mdsumner]

thanks for the reprex I'll try to look more closely soon

[mdsumner]

sorry to have lost this ... the problem is the dates on the raster are bogus

 str(getZ(wc_oisst_crop1))
 chr [1:8766] "1982-01-01 12" "00" "00" "1982-01-02 12" "00" "00" ...

but, I should be robust to that so will put in a fix. Thanks!

[mdsumner]

Actually it's a bug in rerddap, nlayers should match length of getZ

 str(getZ(wc_oisst_crop1))
 chr [1:8766] "1982-01-01 12" "00" "00" "1982-01-02 12" "00" "00" ...
 nlayers(wc_oisst_crop1)
[1] 2922

so a workaround is

z <- as.Date(getZ(wc_oisst_crop1))
 z <- z[!is.na(z)]
range(z)
#[1] "1982-01-01" "1989-12-31"
range(diff(z))
#Time differences in days
#[1] 1 1
wc_oisst_crop_fixed <- setZ(wc_oisst_crop1, z)
as_tibble(wc_oisst_crop_fixed, cell=FALSE, dim=TRUE, values=TRUE, xy=TRUE)
# A tibble: 6,919,296 × 4
   cellvalue dimindex       x     y
       <dbl> <date>     <dbl> <dbl>
 1        NA 1982-01-01 -126.  50.1
 2        NA 1982-01-01 -126.  50.1
 3        NA 1982-01-01 -125.  50.1
 4        NA 1982-01-01 -125.  50.1
 5        NA 1982-01-01 -125.  50.1
 6        NA 1982-01-01 -125.  50.1
 7        NA 1982-01-01 -124.  50.1
 8        NA 1982-01-01 -124.  50.1
 9        NA 1982-01-01 -124.  50.1
10        NA 1982-01-01 -124.  50.1
# … with 6,919,286 more rows
# ℹ Use `print(n = ...)` to see more rows