search

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
558 stars 94 forks source link

write_mdim() error "Error: NA value in which_dims: logic error" , but replacing values works? #701

Open jvandens opened 1 month ago

jvandens commented 1 month ago

I have an issue that might be similar/related to #699

Consider this stars object output from dput():

m = 

structure(list(`F Coli` = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(5L, z = 10L, 
times = 8L)), Entero = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(5L, z = 10L, 
times = 8L))), dimensions = structure(list(geometry = structure(list(
    from = 1L, to = 5L, offset = NA_real_, delta = NA_real_, 
    refsys = structure(list(input = "NAD83 / UTM zone 18N", wkt = "PROJCRS[\"NAD83 / UTM zone 18N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    ID[\"EPSG\",26918]]"), class = "crs"), 
    point = FALSE, values = structure(list(structure(list(structure(c(498756.910095215, 
    510957.220092773, 522344.630126953, 512631.440124512, 498756.910095215, 
    4285731.00012207, 4304042.50012207, 4298020.00012207, 4279347.50012207, 
    4285731.00012207), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
    "sfg")), structure(list(structure(c(512631.440124512, 522344.630126953, 
    533818.380126953, 525469.630126953, 512631.440124512, 4279347.50012207, 
    4298020.00012207, 4292652.50012207, 4273846.50012207, 4279347.50012207
    ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
        structure(list(structure(c(525469.630126953, 533818.380126953, 
        545675.130126953, 538335.190124512, 525469.630126953, 
        4273846.50012207, 4292652.50012207, 4287713.50012207, 
        4268440.50012207, 4273846.50012207), .Dim = c(5L, 2L))), class = c("XY", 
        "POLYGON", "sfg")), structure(list(structure(c(538335.190124512, 
        545675.130126953, 558535.25012207, 552430.810119629, 
        538335.190124512, 4268440.50012207, 4287713.50012207, 
        4283077.50012207, 4262513.00012207, 4268440.50012207), .Dim = c(5L, 
        2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
            structure(c(552430.810119629, 558535.25012207, 572934.75012207, 
            568734.440124512, 552430.810119629, 4262513.00012207, 
            4283077.50012207, 4279116.00012207, 4256515.00012207, 
            4262513.00012207), .Dim = c(5L, 2L))), class = c("XY", 
        "POLYGON", "sfg"))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = 498756.910095215, 
    ymin = 4256515.00012207, xmax = 572934.75012207, ymax = 4304042.50012207
    ), class = "bbox"), crs = structure(list(input = "NAD83 / UTM zone 18N", 
        wkt = "PROJCRS[\"NAD83 / UTM zone 18N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    ID[\"EPSG\",26918]]"), class = "crs"), n_empty = 0L)), class = "dimension"), 
    z = structure(list(from = 1L, to = 10L, offset = 1, delta = 1, 
        refsys = NA_character_, point = FALSE, values = NULL), class = "dimension"), 
    times = structure(list(from = 1L, to = 8L, offset = NA_real_, 
        delta = NA_real_, refsys = NA_character_, point = FALSE, 
        values = structure(list(start = c(0, 0.020833333954215, 
        0.0625, 0.10416666418314, 0.145833343267441, 0.1875, 
        0.22916667163372, 0.270833313465118), end = c(0.020833333954215, 
        0.0625, 0.10416666418314, 0.145833343267441, 0.1875, 
        0.22916667163372, 0.270833313465118, 0.3125)), class = "intervals")), class = "dimension")), raster = structure(list(
    affine = c(0, 0), dimensions = c(NA_character_, NA_character_
    ), curvilinear = FALSE, blocksizes = NULL), class = "stars_raster"), class = "dimensions"), class = "stars")

When I try to save it with write_mdim() to netcdf I get the following error:

write_mdim(m, "m.nc")  # Fails
# Error: NA value in which_dims: logic error

But, if I simply replace the values of the attributes with themselves, it saves:

m2=m
m2$`F Coli` = m$`F Coli` |> as.vector()
m2$Entero = m$Entero |> as.vector()

# write
write_mdim(m2, "m2.nc")  # Succeeds
m_test <- read_mdim("m2.nc")
m_test

# but note that is False, but seems like they are the same?
identical(m2, m_test) # False?

Any ideas why this happening or what could be "wrong" with the attribute values?