# Package sample data
test_that("sample data", {
  # ERA5land: Y-axis runs from top to bottom
  fn <- system.file("extdata", "ERA5land_Rwanda_20160101.nc", package = "ncdfCF")
  ds <- open_ncdf(fn)
  expect_equal(names(ds), c("t2m", "pev", "tp"))
  t2m <- ds[["t2m"]]
  expect_equal(names(t2m$axes), c("longitude", "latitude", "time"))
  lat <- t2m$axes[["latitude"]]
  lat_crds <- lat$coordinates
  expect_length(lat_crds, 21)
  expect_equal(lat_crds[1], -1.0)
  expect_equal(lat_crds[21], -3.0)

  raw <- t2m$data()$raw()
  dn <- dimnames(raw)
  expect_equal(names(dn), c("longitude", "latitude", "time"))
  expect_equal(as.numeric(dn[["latitude"]][1]), lat_crds[1])
  expect_equal(as.numeric(dn[["latitude"]][21]), lat_crds[21])

  arr <- t2m$data()$array()
  dn <- dimnames(arr)
  expect_equal(names(dn), c("latitude", "longitude", "time"))
  expect_equal(as.numeric(dn[["latitude"]][1]), lat_crds[1])
  expect_equal(as.numeric(dn[["latitude"]][21]), lat_crds[21])

  # CORDEX NAM-44 data for a single day; auxiliary grids
  fn <- system.file("extdata", "tasmax_NAM-44_day_20410701-vncdfCF.nc", package = "ncdfCF")
  ds <- open_ncdf(fn)
  tm <- ds[["tasmax"]]
  expect_equal(sapply(tm$axes, function (x) x$length), c(x = 148, y = 140, time = 1, height = 1))
  sub <- tm$subset(X = -90:-70, Y = 40:50)
  expect_equal(sapply(sub$axes, function (x) x$length), c(lon = 31, lat = 23, time = 1, height = 1))
  expect_equal(sub$attribute("coordinates"), "height time")
})

# # Offline data, not available to CRAN or Github
# test_that("CRU data", {
#   cru <- "../testdata/cru_ts4.08.2011.2020.tmp.dat.nc"
#   if (file.exists(cru)) {
#     ds <- open_ncdf(cru)
#
#     # Test basic properties of data variables and axes
#     expect_equal(names(ds), c("tmp", "stn", "mae", "maea"))
#     tmp <- ds[["tmp"]]
#     expect_equal(class(tmp)[1L], "CFVariable")
#     expect_equal(names(tmp$axes), c("lon", "lat", "time"))
#
#     lon <- tmp$axes[["lon"]]
#     expect_equal(class(lon)[1L], "CFAxisLongitude")
#     expect_equal(lon$length, 720L)
#     expect_equal(lon$range(), c(-179.75, 179.75))
#     expect_equal(lon$coordinates, seq(from = -179.75, by = 0.5, length = 720))
#
#     lat <- tmp$axes[["lat"]]
#     expect_equal(class(lat)[1L], "CFAxisLatitude")
#     expect_equal(lat$length, 360L)
#     expect_equal(lat$range(), c(-89.75, 89.75))
#     expect_equal(lat$coordinates, seq(from = -89.75, by = 0.5, length = 360))
#     expect_null(lat$time())
#
#     time <- tmp$axes[["time"]]
#     expect_equal(class(time)[1L], "CFAxisTime")
#     expect_equal(time$length, 120L)
#     tm <- time$time()
#     expect_equal(class(tm)[1L], "CFTime")
#     expect_equal(tm$range(), c("2011-01-16", "2020-12-16"))
#
#     # Extract raw data
#     raw <- tmp[]
#     expect_equal(dim(raw), c(720, 360, 120))
#     dn <- dimnames(raw)
#     expect_equal(names(dn), c("lon", "lat", "time"))
#     ax <- attr(raw, "axis")
#     expect_equal(ax, c(lon = "X", lat = "Y", time = "T"))
#     t <- attr(raw, "time")[[1L]]
#     expect_true(inherits(t, "CFTime"))
#     expect_equal(t$calendar$origin, data.frame(year = 1900, month = 1, day = 1, hour = 0, minute = 0, second = 0, tz = "+0000", offset = 0))
#     expect_length(t, 120)
#     expect_equal(range(t$offsets), c(40557, 44179))
#     expect_equal(t$range(), range(dn$time))
#     expect_null(t$bounds)
#
#     raw <- tmp[361:380, 271:300, ]
#     expect_equal(dim(raw), c(20, 30, 120))
#     dn <- dimnames(raw)
#     expect_equal(names(dn), c("lon", "lat", "time"))
#     ax <- attr(raw, "axis")
#     expect_equal(ax, c(lon = "X", lat = "Y", time = "T"))
#     t <- attr(raw, "time")[[1L]]
#     expect_true(inherits(t, "CFTime"))
#     expect_equal(range(as.numeric(dn$lon)), c(0.25, 9.75))
#     expect_equal(range(as.numeric(dn$lat)), c(45.25, 59.75))
#
#     raw <- tmp[, 271:300, ]
#     expect_equal(dim(raw), c(720, 30, 120))
#     raw <- tmp[361:380, 271:300, 25:36]
#     expect_equal(dim(raw), c(20, 30, 12))
#     raw <- tmp[, , 4]
#     expect_equal(dim(raw), c(720, 360, 1))
#
#     # Subsetting
#     sub <- tmp$subset(lon = -65:-45, lat = -40:-30)
#     expect_true(inherits(sub, "CFArray"))
#     axes <- sub$axes
#     expect_length(axes, 3)
#     expect_equal(names(axes), c("lon", "lat", "time"))
#     expect_equal(range(axes$lon$coordinates), c(-64.75, -45.25))
#     expect_equal(range(axes$lat$coordinates), c(-39.75, -30.25))
#     expect_equal(range(axes$time$coordinates), c("2011-01-16", "2020-12-16"))
#     actual_range <- sub$attribute("actual_range")
#     expect_length(actual_range, 2)
#
#     a <- aoi(-60, -50, -35, -20)
#     sub <- tmp$subset(.aoi = a)
#     expect_true(inherits(sub, "CFArray"))
#     axes <- sub$axes
#     expect_length(axes, 3)
#     expect_equal(names(axes), c("lon", "lat", "time"))
#     expect_equal(range(axes$lon$coordinates), c(-59.75, -50.25))
#     expect_equal(range(axes$lat$coordinates), c(-34.75, -20.25))
#     expect_equal(range(axes$time$coordinates), c("2011-01-16", "2020-12-16"))
#     actual_range <- sub$attribute("actual_range")
#     expect_length(actual_range, 2)
#
#     raw <- sub$raw()
#     expect_equal(dim(raw), c(20, 30, 120))
#     dn <- dimnames(raw)
#     expect_equal(names(dn), c("lon", "lat", "time"))
#     expect_equal(range(as.numeric(dn$lon)), c(-59.75, -50.25))
#     expect_equal(range(as.numeric(dn$lat)), c(-34.75, -20.25))
#     expect_equal(range(dn$time), c("2011-01-16", "2020-12-16"))
#
#     arr <- sub$array()
#     expect_equal(dim(arr), c(30, 20, 120))
#     dn <- dimnames(arr)
#     expect_equal(names(dn), c("lat", "lon", "time"))
#     expect_equal(dn$lat[1], "-20.25")
#     expect_equal(range(as.numeric(dn$lat)), c(-34.75, -20.25))
#     expect_equal(range(as.numeric(dn$lon)), c(-59.75, -50.25))
#     expect_equal(range(dn$time), c("2011-01-16", "2020-12-16"))
#
#     # Profile data
#     prof <- tmp$profile(lon = 5, lat = 50)
#     expect_true(inherits(prof, "CFArray"))
#     expect_equal(prof$name, "location_1")
#     expect_equal(prof$axes[["lon"]]$length, 1)
#     expect_equal(prof$axes[["lat"]]$length, 1)
#     expect_equal(prof$axes[["time"]]$length, 120)
#     expect_equal(prof$axes[["lon"]]$coordinates, 5)
#     expect_equal(prof$axes[["lat"]]$coordinates, 50)
#     expect_equal(prof$attribute("coordinates"), "lon lat")
#
#     prof <- tmp$profile(lat = c(5, 10, 15), lon = c(30, 30, 30), .names = c("South Sudan", "Sudan - South Kordofan"))
#     expect_true(is.list(prof))
#     expect_true(all(sapply(prof, inherits, "CFArray")))
#     expect_equal(names(prof), c("South Sudan", "Sudan - South Kordofan", "location_3"))
#     expect_equal(sapply(prof, function(x) x$axes[["lon"]]$length), c(`South Sudan` = 1, `Sudan - South Kordofan` = 1, location_3 = 1))
#     expect_equal(sapply(prof, function(x) x$axes[["lat"]]$length), c(`South Sudan` = 1, `Sudan - South Kordofan` = 1, location_3 = 1))
#     expect_equal(sapply(prof, function(x) x$axes[["time"]]$length), c(`South Sudan` = 120, `Sudan - South Kordofan` = 120, location_3 = 120))
#
#     prof <- tmp$profile(lon = 5, .names = "Longitude 5 degrees")
#     expect_true(inherits(prof, "CFArray"))
#     expect_equal(prof$name, "Longitude 5 degrees")
#     expect_equal(prof$axes[["lon"]]$length, 1)
#     expect_equal(prof$axes[["lat"]]$length, 360)
#     expect_equal(prof$axes[["time"]]$length, 120)
#     expect_equal(prof$axes[["lon"]]$coordinates, 5)
#     expect_equal(prof$attribute("coordinates"), "lon")
#
#     # Summarise data
#     summ <- tmp$summarise("Tmean", mean, "season") # tmp = CFVariable
#     expect_equal(summ$name, "Tmean")
#     expect_equal(summ$axes[["time"]]$length, 41)
#     expect_length(summ$attribute("actual_range"), 2)
#
#     summ <- sub$summarise("Tmean", mean, "season") # sub = CFArray
#     expect_equal(summ$name, "Tmean")
#     expect_equal(summ$axes[["time"]]$length, 41)
#     expect_length(summ$attribute("actual_range"), 2)
#
#     summ <- prof$summarise("Tmean", mean, "quarter") # prof is profile in CFArray
#     expect_equal(summ$axes[["time"]]$length, 40)
#     expect_length(summ$attribute("actual_range"), 2)
#
#     summ <- sub$summarise(c("Tmin", "Tmax"), range, "year", era = list(early = 2011:2016, late = 2017:2020)) # sub = CFArray
#     expect_true(is.list(summ))
#     expect_equal(names(summ), c("early", "late"))
#     expect_equal(names(summ$early), c("Tmin", "Tmax"))
#     expect_equal(names(summ$late), c("Tmin", "Tmax"))
#     Tearly <- summ$early$Tmin$axes[["time"]]
#     expect_equal(Tearly$length, 6)
#     expect_equal(summ$late$Tmin$axes[["time"]]$length, 4)
#     expect_length(Tearly$attribute("actual_range"), 2)
#     expect_true(inherits(Tearly$time(), "CFClimatology"))
#
#     # Save to file, read back in
#     fn <- tempfile(fileext = ".nc")
#     Tmin <- summ$early$Tmin
#     Tmin$save(fn)
#     ds2 <- open_ncdf(fn)
#     expect_equal(names(ds2), "Tmin")
#     T2 <- ds2[["Tmin"]]$data()
#     expect_equal(T2$raw(), aperm(Tmin$raw(), c(2, 3, 1)))
#     expect_equal(T2$axes[["lon"]]$coordinates, Tmin$axes[["lon"]]$coordinates)
#     expect_equal(T2$axes[["lat"]]$coordinates, Tmin$axes[["lat"]]$coordinates)
#     expect_equal(T2$axes[["time"]]$coordinates, Tmin$axes[["time"]]$coordinates)
#     ds2 <- NULL
#     unlink(fn)
#   }
# })
#
# test_that("CMEMS example online", {
#   cmems <- "../testdata/cmems_mod_med_phy-tem_anfc_4.2km_P1M-m_1748594717636.nc"
#   if (file.exists(cmems)) {
#     suppressWarnings(ds <- open_ncdf(cmems))
#
#     thetao <- ds[["thetao"]]
#     tc <- thetao$axes[["time"]]$coordinates
#     dc <- thetao$axes[["depth"]]$coordinates
#     expect_length(tc, 12)
#     expect_length(dc, 141)
#
#     surfaceT <- thetao$subset(depth = c(0, 2))
#     expect_length(surfaceT$axes, 4)
#     expect_length(surfaceT$axes[["depth"]]$coordinates, 1)
#
#     # Get the longitude and latitude coordinates
#     longitude <- thetao$axes[["longitude"]]$coordinates
#     latitude <- thetao$axes[["latitude"]]$coordinates
#
#     # Make a map of the first time slice: January
#     # library(RColorBrewer)
#     # image(longitude, latitude, surfaceT$raw()[,,1], col = rev(brewer.pal(10, "RdBu")))
#
#     # library(lattice)
#     # grid <- expand.grid(lon=longitude, lat=latitude)
#     # cutpts <- c(12,13,14,15,16,17,18,19,20)
#     # levelplot(surfaceT$raw()[,,1] ~ lon * lat,
#     #           data=grid, region=TRUE,
#     #           pretty=T, at=cutpts, cuts=9,
#     #           col.regions=(rev(brewer.pal(9,"RdBu"))), contour=0,
#     #           xlab = "Longitude", ylab = "Latitude",
#     #           main = "Sea Water Potential Temperature (°C)"
#     # )
#
#     # Make a vertical profile of temperature at a single location for the month of July 2024
#     tprof <- thetao$profile(longitude = 12.71, latitude = 40.33, time = "2024-07-01", .names = "Tyrrhenian_Sea")
#     expect_equal(sapply(tprof$axes, function(ax) ax$length), c(depth = 141, longitude = 1, latitude = 1, time = 1))
#     Tjuly <- tprof$data.table()
#     expect_length(Tjuly, 5)
#     expect_equal(nrow(Tjuly), 141)
#
#     # suppressPackageStartupMessages(library(ggplot2))
#     # ggplot(Tjuly) + geom_line(aes(depth, Tyrrhenian_Sea)) + theme_bw() + xlim(0, 4000) +
#     #  labs(x = "Depth (m)", y = "Temperature (°C)", title =  "Tyrrhenian Sea", subtitle = "12.71E - 40.33N")
#   }
# })
#
# test_that("Auxiliary grids", {
#   aux <- "../testdata/pr_EUR-11_MOHC-HadGEM2-ES_rcp45_r1i1p1_CLMcom-CCLM4-8-17_v1_day_20060101-20101230.nc"
#   if (file.exists(aux)) {
#     ds <- open_ncdf(aux)
#     pr <- ds[["pr"]]
#     expect_equal(pr$attribute("coordinates"), "lat lon")
#     expect_equal(pr$attribute("grid_mapping"), "rotated_latitude_longitude")
#     sel <- pr$subset(X = -10:5, Y = 45:60, time = c("2009-01-01", "2010-01-01"))
#     expect_true(is.na(sel$attribute("coordinates")))
#     expect_true(is.na(sel$attribute("grid_mapping")))
#     lon_crds <- sel$axes[["lon"]]$coordinates
#     lon_res <- diff(lon_crds[1:2])
#     lat_crds <- sel$axes[["lat"]]$coordinates
#     lat_res <- diff(lat_crds[1:2])
#     expect_true(abs(lon_crds[1] - -10) < lon_res)
#     expect_true(abs(lat_crds[1] - 45) < lat_res)
#     expect_true(abs(lon_crds[length(lon_crds)] - 5) < lon_res)
#     expect_true(abs(lat_crds[length(lat_crds)] - 60) < lat_res)
#
#     a <- aoi(-10, 5, 45, 60, 0.2)
#     sel <- pr$subset(time = c("2009-01-01", "2010-01-01"), .aoi = a)
#     expect_true(is.na(sel$attribute("coordinates")))
#     expect_true(is.na(sel$attribute("grid_mapping")))
#     lon_crds <- sel$axes[["lon"]]$coordinates
#     lon_res <- diff(lon_crds[1:2])
#     lat_crds <- sel$axes[["lat"]]$coordinates
#     lat_res <- diff(lat_crds[1:2])
#     expect_equal(lon_res, 0.2)
#     expect_equal(lat_res, 0.2)
#     expect_equal(lon_crds[1], -9.9)
#     expect_equal(lat_crds[1], 45.1)
#     expect_equal(lon_crds[length(lon_crds)], 4.9)
#     expect_equal(lat_crds[length(lat_crds)], 59.9)
#     bnds <- sel$axes[["lon"]]$bounds
#     expect_equal(bnds$coordinates[1, 1], -10)
#     expect_equal(bnds$coordinates[2, 75], 5)
#     bnds <- sel$axes[["lat"]]$bounds
#     expect_equal(bnds$coordinates[1, 1], 45)
#     expect_equal(bnds$coordinates[2, 75], 60)
#   }
# })
#
# test_that("Multiple label sets", {
#   mme <- "../testdata/t_CMIP6_ssp245_mon_201501-210012.nc"
#   if (file.exists(mme)) {
#     ds <- open_ncdf(mme)
#     t <- ds[["t"]]
#     expect_length(t$axes, 5)
#     member <- t$axes[["member"]]
#     expect_equal(member$coordinate_names, c("member", "member_id", "gcm_institution", "gcm_model", "gcm_variant"))
#     expect_equal(member$active_coordinates, "member")
#     sel <- t$subset(time = "2020-01-01", member = 4)
#     expect_true(inherits(sel, "CFArray"))
#     expect_equal(sapply(sel$axes, function(x) x$length), c(lon = 360, lat = 180, time = 1, member = 1, height2m = 1))
#
#     member$active_coordinates <- "gcm_model"
#     expect_equal(member$coordinates[1], "ACCESS-CM2")
#     sel <- t$subset(time = "2020-01-01", member = c("CAMS-CSM1-0", "CESM2"))
#     expect_true(inherits(sel, "CFArray"))
#     expect_equal(sapply(sel$axes, function(x) x$length), c(lon = 360, lat = 180, member = 3, time = 1, height2m = 1))
#     expect_equal(sel$axes[["member"]]$coordinates, c("CAMS-CSM1-0", "CESM2-WACCM", "CESM2"))
#   }
# })
#
# test_that("Merge CFArrays", {
#   esm4 <- list.files("../testdata", "^pr_day_GFDL.*nc$", full.names = T)
#   if (length(esm4)) {
#     expect_length(esm4, 4) # Just to make sure the files are captured
#     ds <- open_ncdf(esm4[1])
#     combi <- ds[["pr"]]$subset(lon = 0:10, lat = 0:20)
#     expect_equal(names(combi$axes), c("lon", "lat", "time"))
#     crds_lon <- combi$axes[["lon"]]$coordinates
#     crds_lat <- combi$axes[["lat"]]$coordinates
#     len_time <- combi$axes[["time"]]$length
#
#     for (i in 2:4) {
#       ds <- open_ncdf(esm4[i])
#       pr <- ds[["pr"]]$subset(lon = 0:10, lat = 0:20)
#       expect_equal(pr$axes[["lon"]]$coordinates, crds_lon)
#       expect_equal(pr$axes[["lat"]]$coordinates, crds_lat)
#       len_time <- len_time + pr$axes[["time"]]$length
#       combi$append(pr, "time")
#     }
#
#     expect_equal(combi$axes[["time"]]$length, len_time)
#     expect_equal(range(diff(combi$axes[["time"]]$values)), c(1, 1))
#   } else cat("\nCheck esm4\n")
# })