How to read AMSR2 sea ice data with terra

R
Geospatial
Author

Philippe Massicotte

Published

August 11, 2022

Today I was trying to read AMSR2 sea ice data. I was surprised to discover that the files do not include coordinates or projection information. Data and coordinates are contained in different files! Maybe to save some disk space ⁉️ The sea ice data is (for example) in a file named Arc_20201010_res3.125_pyres.nc.gz whereas the coordinates are included in a file name LongitudeLatitudeGrid_3.125km_Arctic.nc. Furthermore, I found that the provided coordinates are provided in long/lat format whereas the gridded data are projected 😠

Photo by Willian Justen de Vasconcellos on Unsplash

If I read the data directly, one can see that there are no coordinates or projection information.

library(terra)
library(ggplot2)
library(tidyterra)
library(magick)

r <- rast(
  "/vsigzip//vsicurl/ftp://ftp-projects.cen.uni-hamburg.de/seaice/AMSR2/3.125km/Arc_20201010_res3.125_pyres.nc.gz",
  "sea_ice_concentration"
)

# No projection information
r
#> class       : SpatRaster 
#> dimensions  : 3584, 2432, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0.5, 2432.5, 0.5, 3584.5  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : Arc_20201010_res3.125_pyres.nc.gz:sea_ice_concentration 
#> varname     : sea_ice_concentration (daily averaged total ice concentration) 
#> name        : sea_ice_concentration 
#> unit        :                     % 
#> time        : 2020-10-10 12:00:00 UTC

Fortunately, I found this post by Michael Sumner which guided me on how to manipulate this data. What we have to do is to set the extent and the proper projection after the file is read. This can be done using the ext() and crs() function from the terra package. Based on the documentation, we can manually set the proper values.

Erratum

I was wrongly using the following extent: extent(-3837500, 3762500, -5362500, 5837500). Michael Sumner kindly informed me about the error. The proper extent to use is down below.

# Set the extent
ext(r) <- ext(-3850000, 3750000, -5350000, 5850000)

# Set the polar stereographic projection
crs(r) <- "EPSG:3413"

Now, if we take a look at the raster, we can see that it is correctly projected and has a resolution of 3.125 km as expected ✌️

r
#> class       : SpatRaster 
#> dimensions  : 3584, 2432, 1  (nrow, ncol, nlyr)
#> resolution  : 3125, 3125  (x, y)
#> extent      : -3850000, 3750000, -5350000, 5850000  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / NSIDC Sea Ice Polar Stereographic North (EPSG:3413) 
#> source      : Arc_20201010_res3.125_pyres.nc.gz:sea_ice_concentration 
#> varname     : sea_ice_concentration (daily averaged total ice concentration) 
#> name        : sea_ice_concentration 
#> unit        :                     % 
#> time        : 2020-10-10 12:00:00 UTC

Finally, visualize the raster.

# Set values of 0 to NA
NAflag(r) <- 0

ggplot() +
  geom_spatraster(data = r / 100) +
  scale_fill_viridis_c(
    na.value = "transparent",
    labels = scales::label_percent(),
    breaks = scales::breaks_pretty(n = 6),
    guide = guide_colorbar(
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(5, "cm"),
      barheight = unit(0.2, "cm")
    )
  ) +
  labs(
    fill = "Sea ice concentration"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    panel.grid = element_line(size = 0.25)
  )

Session info 
#> ═ Session info ═══════════════════════════════════════════════════════════════════════════════════════════════════════
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  ! package     * version date (UTC) lib source
#>  P cachem        1.0.8   2023-05-01 [?] RSPM
#>  P class         7.3-22  2023-05-03 [?] CRAN (R 4.3.1)
#>  P classInt      0.4-10  2023-09-05 [?] RSPM
#>  P cli           3.6.2   2023-12-11 [?] RSPM
#>  P codetools     0.2-19  2023-02-01 [?] CRAN (R 4.2.2)
#>  P colorspace    2.1-0   2023-01-23 [?] RSPM
#>  P data.table    1.15.4  2024-03-30 [?] RSPM
#>  P DBI           1.2.2   2024-02-16 [?] RSPM
#>  P devtools      2.4.5   2022-10-11 [?] RSPM (R 4.4.0)
#>  P digest        0.6.35  2024-03-11 [?] RSPM
#>  P dplyr         1.1.4   2023-11-17 [?] RSPM
#>  P e1071         1.7-14  2023-12-06 [?] RSPM
#>  P ellipsis      0.3.2   2021-04-29 [?] RSPM
#>  P evaluate      0.23    2023-11-01 [?] RSPM
#>  P fansi         1.0.6   2023-12-08 [?] RSPM
#>  P farver        2.1.1   2022-07-06 [?] RSPM
#>  P fastmap       1.1.1   2023-02-24 [?] RSPM
#>  P fs            1.6.4   2024-04-25 [?] CRAN (R 4.4.0)
#>  P generics      0.1.3   2022-07-05 [?] RSPM
#>  P ggplot2     * 3.5.1   2024-04-23 [?] RSPM
#>  P glue          1.7.0   2024-01-09 [?] RSPM
#>  P gtable        0.3.5   2024-04-22 [?] RSPM
#>  P htmltools     0.5.8.1 2024-04-04 [?] RSPM
#>  P htmlwidgets   1.6.4   2023-12-06 [?] RSPM
#>  P httpuv        1.6.15  2024-03-26 [?] RSPM
#>  P jsonlite      1.8.8   2023-12-04 [?] RSPM
#>  P KernSmooth    2.23-22 2023-07-10 [?] CRAN (R 4.3.1)
#>  P knitr         1.46    2024-04-06 [?] RSPM
#>  P later         1.3.2   2023-12-06 [?] RSPM
#>  P lifecycle     1.0.4   2023-11-07 [?] RSPM
#>  P magick      * 2.8.3   2024-02-18 [?] RSPM
#>  P magrittr      2.0.3   2022-03-30 [?] RSPM
#>  P memoise       2.0.1   2021-11-26 [?] RSPM
#>  P mime          0.12    2021-09-28 [?] RSPM
#>  P miniUI        0.1.1.1 2018-05-18 [?] RSPM (R 4.4.0)
#>  P munsell       0.5.1   2024-04-01 [?] RSPM
#>  P pillar        1.9.0   2023-03-22 [?] RSPM
#>  P pkgbuild      1.4.4   2024-03-17 [?] RSPM (R 4.4.0)
#>  P pkgconfig     2.0.3   2019-09-22 [?] RSPM
#>  P pkgload       1.3.4   2024-01-16 [?] RSPM (R 4.4.0)
#>  P processx      3.8.4   2024-03-16 [?] RSPM
#>  P profvis       0.3.8   2023-05-02 [?] RSPM (R 4.4.0)
#>  P promises      1.3.0   2024-04-05 [?] RSPM
#>  P proxy         0.4-27  2022-06-09 [?] RSPM
#>  P ps            1.7.6   2024-01-18 [?] RSPM
#>  P purrr         1.0.2   2023-08-10 [?] RSPM
#>  P quarto      * 1.4     2024-03-06 [?] RSPM
#>  P R.cache       0.16.0  2022-07-21 [?] RSPM
#>  P R.methodsS3   1.8.2   2022-06-13 [?] RSPM
#>  P R.oo          1.26.0  2024-01-24 [?] RSPM
#>  P R.utils       2.12.3  2023-11-18 [?] RSPM
#>  P R6            2.5.1   2021-08-19 [?] RSPM
#>  P Rcpp          1.0.12  2024-01-09 [?] RSPM
#>  P remotes       2.5.0   2024-03-17 [?] RSPM (R 4.4.0)
#>  P renv          1.0.7   2024-04-11 [?] RSPM (R 4.4.0)
#>  P rlang         1.1.3   2024-01-10 [?] RSPM
#>  P rmarkdown     2.26    2024-03-05 [?] RSPM
#>  P rstudioapi    0.16.0  2024-03-24 [?] RSPM
#>  P scales        1.3.0   2023-11-28 [?] RSPM
#>  P sessioninfo   1.2.2   2021-12-06 [?] RSPM (R 4.4.0)
#>  P sf            1.0-16  2024-03-24 [?] RSPM
#>  P shiny         1.8.1.1 2024-04-02 [?] RSPM (R 4.4.0)
#>  P stringi       1.8.3   2023-12-11 [?] RSPM
#>  P stringr       1.5.1   2023-11-14 [?] RSPM
#>  P styler      * 1.10.3  2024-04-07 [?] RSPM
#>  P terra       * 1.7-71  2024-01-31 [?] RSPM
#>  P tibble        3.2.1   2023-03-20 [?] RSPM
#>  P tidyr         1.3.1   2024-01-24 [?] RSPM
#>  P tidyselect    1.2.1   2024-03-11 [?] RSPM
#>  P tidyterra   * 0.6.0   2024-04-22 [?] RSPM
#>  P tinytex       0.50    2024-03-16 [?] RSPM
#>  P units         0.8-5   2023-11-28 [?] RSPM
#>  P urlchecker    1.0.1   2021-11-30 [?] RSPM (R 4.4.0)
#>  P usethis       2.2.3   2024-02-19 [?] RSPM (R 4.4.0)
#>  P utf8          1.2.4   2023-10-22 [?] RSPM
#>  P vctrs         0.6.5   2023-12-01 [?] RSPM
#>  P viridisLite   0.4.2   2023-05-02 [?] RSPM
#>  P withr         3.0.0   2024-01-16 [?] RSPM
#>  P xfun          0.43    2024-03-25 [?] RSPM
#>  P xtable        1.8-4   2019-04-21 [?] RSPM (R 4.4.0)
#>  P yaml          2.3.8   2023-12-11 [?] RSPM
#> 
#>  [1] /tmp/RtmpbmRvXN/renv-use-libpath-25a668312ef1d2
#>  [2] /tmp/RtmpbmRvXN/renv-sandbox
#> 
#>  P ── Loaded and on-disk path mismatch.
#> 
#> ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────