r-spatial
Support for `stars` with `units`
At the moment, the mapview function doesn’t seems to work with stars objects when raster values are of class units. When converting to Raster*, mapview does work, since the units information is dropped during the conversion. Perhaps stars+units objects may be supported, dropping the units information if any. Thanks!
library(mapview)
library(stars)
## Loading required package: abind
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(units)
## udunits system database from /usr/share/xml/udunits
# Data
x = read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
x
## stars object with 3 dimensions and 1 attribute
## attribute(s):
## L7_ETMs.tif
## Min. : 1.00
## 1st Qu.: 54.00
## Median : 69.00
## Mean : 68.91
## 3rd Qu.: 86.00
## Max. :255.00
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 349 288776 28.5 UTM Zone 25, Southern Hem... FALSE NULL [x]
## y 1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE NULL [y]
## band 1 6 NA NA NA NA NULL
# Set units
x[[1]] = set_units(x[[1]], "m")
x
## stars object with 3 dimensions and 1 attribute
## attribute(s):
## L7_ETMs.tif [m]
## Min. : 1.00
## 1st Qu.: 54.00
## Median : 69.00
## Mean : 68.91
## 3rd Qu.: 86.00
## Max. :255.00
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 349 288776 28.5 UTM Zone 25, Southern Hem... FALSE NULL [x]
## y 1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE NULL [y]
## band 1 6 NA NA NA NA NULL
# Map
mapview(x)
## Error in Ops.units(lim, prop * d * c(-1, 1)): both operands of the expression should be "units" objects
mapview(as(x, "Raster"))
Session info:
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] units_0.6-7 stars_0.4-4 sf_0.9-6 abind_1.4-5 mapview_2.9.4
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.17 purrr_0.3.4
## [4] lattice_0.20-41 colorspace_1.4-1 vctrs_0.3.4
## [7] generics_0.0.2 leaflet.providers_1.9.0 htmltools_0.5.0
## [10] stats4_4.0.2 yaml_2.2.1 base64enc_0.1-3
## [13] rlang_0.4.7 e1071_1.7-3 pillar_1.4.6
## [16] glue_1.4.2 DBI_1.1.0 sp_1.4-2
## [19] lifecycle_0.2.0 stringr_1.4.0 munsell_0.5.0
## [22] raster_3.3-13 htmlwidgets_1.5.1 codetools_0.2-16
## [25] evaluate_0.14 knitr_1.29 callr_3.4.4
## [28] ps_1.3.4 crosstalk_1.1.0.1 parallel_4.0.2
## [31] class_7.3-17 leafem_0.1.3 Rcpp_1.0.5
## [34] KernSmooth_2.23-17 scales_1.1.1 classInt_0.4-3
## [37] satellite_1.0.2 lwgeom_0.2-5 jsonlite_1.7.1
## [40] leaflet_2.0.3 webshot_0.5.2 farver_2.0.3
## [43] png_0.1-7 digest_0.6.25 stringi_1.5.3
## [46] processx_3.4.4 dplyr_1.0.2 grid_4.0.2
## [49] rgdal_1.5-16 tools_4.0.2 magrittr_1.5
## [52] tibble_3.0.3 crayon_1.3.4 pkgconfig_2.0.3
## [55] ellipsis_0.3.1 rmarkdown_2.3 R6_2.4.1
## [58] compiler_4.0.2
Thanks! Dropping the units should be fine as we don't really care for those when plotting?
Guess that they can be dropped. If there should be any effect, at all, then perhaps just adding a label (such as "[m]") in the legend title.
That is of course the easiest. The whole point of passing through measurement units however is that the automatically show up in plots; users could use the drop_units methods for stars objects to suppress this.
Right now the title of the legend is the name of the object, not the name of the first attribute shown. That, along with the units, would make more sense IMO.
When doing
x = c(
"avhrr-only-v2.19810901.nc",
"avhrr-only-v2.19810902.nc",
"avhrr-only-v2.19810903.nc",
"avhrr-only-v2.19810904.nc",
"avhrr-only-v2.19810905.nc",
"avhrr-only-v2.19810906.nc",
"avhrr-only-v2.19810907.nc",
"avhrr-only-v2.19810908.nc",
"avhrr-only-v2.19810909.nc"
)
# see the second vignette:
# install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", type = "source")
file_list = system.file(paste0("netcdf/", x), package = "starsdata")
(y = read_stars(file_list, quiet = TRUE))
then
st_crs(y)=4326
mapview(drop_units(adrop(y)))
gives drop_units(adrop(y)) - 1 as legend title. When I do
plot(adrop(y)[,,,1])
the title is sst [° C]. Makes more sense.
@edzer thanks for the example! I totally agree. Installing starsdata right now.................... (takes a while) :-)
This is a bit more involved, as it means we need to re-think how we use layer names, groups, layerId, etc. I am currently using the layer name as the layerId (if not provided) which is also used as the file name when writing to disk and attaching to the html. Including the units currently means that writing to disk may error as we may have some special character such as '°' in the file name... Thus, we need to separate the names used for disk I/O and for display, which means adjustments in many places I fear... But I agree that this is what should happen.
@edzer I am looking for a way to somehow consistently extract a valid 2d (x/y) "slice" from a stars object. I think your flatten function is what I am looking for. Though, with your example above (4 attributes along 9 time slices) this still leaves the question of how to consistently handle these data or how to expose control over the additional dimensions via the API.
Should we expect (or at least allow) the user to provide the object in the form you do above (adrop(y)[, , , 1]) or have additional argument(s) (somehow) to allow for this specification? Should we e.g. provide an argument along similar to read_stars?
I think I need to re-implement the stars method from scratch as things are sufficiently different from the *Raster approach. My vision is also to use stars internally for all raster based methods and simply coerce *Raster objects that come in.