Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ESA Landcover tile mapping issue #209

Closed
karpfen opened this issue Nov 27, 2023 · 0 comments · Fixed by #211
Closed

ESA Landcover tile mapping issue #209

karpfen opened this issue Nov 27, 2023 · 0 comments · Fixed by #211
Assignees
Labels
bug Something isn't working

Comments

@karpfen
Copy link
Collaborator

karpfen commented Nov 27, 2023

I encountered a bug when calculating the ESA Landcover indicator for a portfolio that spans multiple tiles and years. The package downloads the correct tifs and builds the tileindex GPKG as expected, but upon reading the tiles within calc_indicators, it selects the wrong ones. Then it reports that the polygons and tiles do not overlap and returns NA rather than the expected landcovers.

I managed to trace the error back to this code block:

} else { # some tiles share the same bboxes, and others do not, needs proper merging
# We assume here that the tiles present in tileindex have a temporal dimension.
# Thus each timestep should end up in its own layer. Different tiles from
# the same timestep should be spatially merged. We want to avoid merging
# different tile from different timesteps. We thus assume some regularity
# in how the name of a raster file expresses its temporal dimension.
# With this assumption, we can expect the files in tindex to be ordered.
# Thus we retrive the index of all files sharing the same bbox and assume
# that they belong to different timesteps. The files in between these
# indices thus belong to the previous timestep and we can merge these
# as a vrt and later join the bands. We always assign the name of the
# first file as the layername.
unique_bboxes <- unique(unlist(all_bboxes))
layer_index <- which(all_bboxes == unique_bboxes[[1]])
temporal_gap <- layer_index[2] - layer_index[1] - 1
out <- lapply(layer_index, function(j) {
target_files <- tindex$location[j:(j + temporal_gap)]
org_filename <- basename(target_files[1])
filename <- tools::file_path_sans_ext(org_filename)
vrt_name <- tempfile(pattern = sprintf("vrt_%s.vrt", filename))
tmp <- terra::vrt(target_files, filename = vrt_name)
names(tmp) <- org_filename
tmp
})
out <- do.call(c, out)
}

However, I was not able to fix it yet. I'm also afraid to break something else in the process, so I'd prefer to look at this together once you're back in business, @goergen95.

Heres a reprex of what I mean:

library(sf)
library(mapme.biodiversity)

polys <- data.frame (
  id = 1:2,
  geom = c("Polygon ((20.62472622642718889 -17.96824385542232605, 20.62496802397372875 -17.98457342071599641, 20.63657430620789413 -17.98411345363312464, 20.6416520546853377 -17.96709382912841946, 20.62472622642718889 -17.96824385542232605))",
           "Polygon ((-80.95285139309206102 34.04564701142924577, -80.95140060781277214 34.02841503501579723, -80.93785994520627014 34.0312204790189341, -80.94172870595097891 34.04404418455595049, -80.95285139309206102 34.04564701142924577))"
           )
)

polys <- st_as_sf(polys, wkt = "geom", crs = 4326)

pf <- init_portfolio(
    polys,
    2015:2016,
    outdir = "location1",
    verbose = TRUE) %>%
  get_resources("esalandcover")

inds <- calc_indicators(pf, "landcover")
inds$landcover # landcover for one polygon is missing

pf2 <- init_portfolio(
    polys,
    2015,
    outdir = "location2",
    verbose = TRUE) %>%
  get_resources("esalandcover")

inds2 <- calc_indicators(pf2, "landcover")
inds2$landcover # correct result

sessionInfo()


R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: Europe/Berlin
tzcode source: internal

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

other attached packages:
[1] purrr_1.0.2              mapme.biodiversity_0.4.0 sf_1.0-14               

loaded via a namespace (and not attached):
 [1] s2_1.1.4           future_1.33.0      tidyr_1.3.0        utf8_1.2.4         generics_0.1.3     class_7.3-22       KernSmooth_2.23-22 stringi_1.7.12     SPEI_1.8.1        
[10] lattice_0.21-9     listenv_0.9.0      digest_0.6.33      magrittr_2.0.3     lmomco_2.4.11      grid_4.3.1         timechange_0.2.0   plyr_1.8.9         e1071_1.7-13      
[19] backports_1.4.1    reshape_0.8.9      DBI_1.1.3          httr_1.4.7         fansi_1.0.5        scales_1.2.1       codetools_0.2-19   cli_3.6.1          rlang_1.1.1       
[28] units_0.8-4        parallelly_1.36.0  munsell_0.5.0      withr_2.5.2        Lmoments_1.3-1     parallel_4.3.1     tools_4.3.1        lmom_3.0           checkmate_2.2.0   
[37] dplyr_1.1.3        colorspace_2.1-0   ggplot2_3.4.3      globals_0.16.2     curl_5.1.0         vctrs_0.6.3        R6_2.5.1           zoo_1.8-12         proxy_0.4-27      
[46] lifecycle_1.0.4    lubridate_1.9.3    classInt_0.4-10    stringr_1.5.0      MASS_7.3-60        furrr_0.3.1        pkgconfig_2.0.3    progressr_0.14.0   terra_1.7-55      
[55] pillar_1.9.0       TLMoments_0.7.5.3  gtable_0.3.4       glue_1.6.2         Rcpp_1.0.11        tibble_3.2.1       tidyselect_1.2.0   rstudioapi_0.15.0  goftest_1.2-3     
[64] wk_0.9.0           compiler_4.3.1    
@karpfen karpfen added the bug Something isn't working label Nov 27, 2023
@goergen95 goergen95 linked a pull request Nov 29, 2023 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants