st_contour causes "non-numeric argument" error. 'See Also' suggestion not supported.

[SimonDedman]

Hi folks,

I can plot my stars object in ggplot & ggmap fine as a surface using geom_stars.
I'm trying to plot that same object as contours, and it's causing issues:
All_Rasters_Scales.asc

library(stars)
lemons.ras <- read_stars("All_Rasters_Scales.asc")) %>%
  st_set_crs(2958) %>%
  st_transform(4326)
ggplot() +
geom_stars(data = lemons.ras)

Is fine.

ggplot() +
st_contour(x = lemons.ras)

Warning: Incompatible methods ("+.gg", "Ops.data.frame") for "+" Error in ggmap(myMap) + st_contour(lemons.ras) : non-numeric argument to binary operator

(same error for ggmap & ggplot opening calls). I've seen various help online relating to R profiles' default themes, but I've not edited mine ever. Also this fails on those 2 lines, i.e. before theme is called.

The advice in the 'See Also' section of st_contour is:

For polygonizing rasters following grid boundaries, see st_as_sf with arguments as_points=FALSE and merge=TRUE; contour plots contour lines using R's native algorithm (which also plots contour levels)

But:

lemons.sf <- st_as_sf(lemons.ras, as_points = FALSE, merge = TRUE)

Error in st_as_sf.stars(lemons.ras, as_points = FALSE, merge = TRUE) : merge=TRUE and as_points=TRUE unsupported; consider using st_contour for generating contour lines

If true, this suggests that in Adriano's table, only the bottom option is available ("POLYGONS, not merged").

According to getGDALVersionInfo() is have "GDAL 3.2.2, released 2021/03/05"

Thanks in advance.

[edzer]

st_contour() returns an sf object with the contours; you'd have to pass that to geom_sf().

[SimonDedman]

Aaah, ok, thanks Edzer. I don't suppose there's any chance you could add your answer to the st_contour man page? Might just be me, but it seemed intuitive that the returned object would be directly compatible into a ggplot chain, but I can see how that's not the case.

Thanks again. Finally getting into rasters and stars is a real gift to the community.

Edit: sorry to bother you again. Perhaps I've misunderstood something. If I pass to geom_sf directly:

ggmap(myMap) +
  geom_sf(st_contour(x = lemons.ras,
                     contour_lines = FALSE,
                     breaks = c(0.5, 0.95)))

Or via an intermediary

mycontour <- st_contour(as above)
ggmap(myMap) + geom_sf(mycontour)

I get:

Error: mapping must be created by aes()

I was under the impression that this error shouldn't be seen when using sf objects, indeed that this is somewhat the point of the sf ecosystem, since the geometry column is present. Querying the object:

st_geometry(mycontour)

Geometry set for 3 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
MULTIPOLYGON (((NA NA, NA NA, NA NA, NA NA, NA ...
MULTIPOLYGON (((NA NA, NA NA, NA NA, NA NA, NA ...
MULTIPOLYGON (((NA NA, NA NA, NA NA, NA NA, NA ...

Looks 'normal' i.e. the same structure as a coastline sf object I created which plots fine in geom_sf; of possible concern is the NA bounding box. The original stars raster had one:

st_bbox(lemons.ras)

xmin ymin xmax ymax
-79.28970 25.68331 -79.21143 25.78939

But it doesn't look like I can set one:

st_set_bbox(mycontour, st_bbox(lemons.ras))

Error in UseMethod("st_set_bbox") : no applicable method for 'st_set_bbox' applied to an object of class "c('sf', 'data.frame')"

I would assume that NA values in the multipolygon geometries shouldn't cause the mapping error, they'd just not plot?

Thanks again for your thoughts.

[gavg712]

The outputs of st_transform() do not resample the data to new locations but a curvilinear grid is returned, and st_contour() is not set up to work with curvilinear grids. As I have seen, there are two alternatives in this case

1. Use st_warp() to resample the original data into a new regular grid.
2. Transform the contours after st_contour() with the original data as input.

The second option is less loss of information in my point of view. Here an example with your data:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
#> Registered S3 methods overwritten by 'stars':
#>   method             from
#>   st_bbox.SpatRaster sf  
#>   st_crs.SpatRaster  sf
library(ggplot2)
f <- tempfile(fileext = ".asc")
curl::curl_download("https://raw.githubusercontent.com/SimonDedman/dBBMM_HomeRange/main/data/All_Rasters_Scaled.asc",
                    f)

lemons.ras <- read_stars(f) %>%
  st_set_crs(2958)

stars:::is_curvilinear(lemons.ras)
#> [1] FALSE

# without transforming
ggplot() +
  geom_stars(data = lemons.ras, 
             alpha = 0.5) +
  geom_sf(data = st_contour(x = lemons.ras,
                            contour_lines = FALSE), 
          fill = NA, color = "orange")

# Transforming
stars:::is_curvilinear(lemons.ras %>% 
                         st_transform(4326))
#> [1] TRUE

# I will transform on the fly to manage the plotting with ggplot2. You can make permanent transforms.

ggplot() +
  geom_stars(data = lemons.ras %>% st_transform(4326), # Raster layer transform 
             alpha = 0.5) +
  geom_sf(data = st_contour(x = lemons.ras,
                            contour_lines = FALSE) %>%
            st_transform(4326), # Vector transform after st_contour() 
          fill = NA, color = "orange")

Cheers.

^{Created on 2021-11-08 by the reprex package (v2.0.1)}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.2 (2021-11-01)
#>  os       Ubuntu 20.04.3 LTS          
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/Guayaquil           
#>  date     2021-11-08                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date       lib source        
#>  abind       * 1.4-5   2016-07-21 [2] CRAN (R 4.0.2)
#>  assertthat    0.2.1   2019-03-21 [2] CRAN (R 4.0.2)
#>  backports     1.2.1   2020-12-09 [2] CRAN (R 4.0.3)
#>  class         7.3-19  2021-05-03 [4] CRAN (R 4.0.5)
#>  classInt      0.4-3   2020-04-07 [2] CRAN (R 4.0.2)
#>  cli           3.0.1   2021-07-17 [2] CRAN (R 4.1.0)
#>  colorspace    2.0-2   2021-06-24 [2] CRAN (R 4.1.0)
#>  crayon        1.4.1   2021-02-08 [2] CRAN (R 4.0.4)
#>  curl          4.3.2   2021-06-23 [2] CRAN (R 4.1.0)
#>  DBI           1.1.1   2021-01-15 [2] CRAN (R 4.0.3)
#>  digest        0.6.28  2021-09-23 [2] CRAN (R 4.1.1)
#>  dplyr         1.0.7   2021-06-18 [2] CRAN (R 4.1.0)
#>  e1071         1.7-9   2021-09-16 [2] CRAN (R 4.1.1)
#>  ellipsis      0.3.2   2021-04-29 [2] CRAN (R 4.0.5)
#>  evaluate      0.14    2019-05-28 [2] CRAN (R 4.0.2)
#>  fansi         0.5.0   2021-05-25 [2] CRAN (R 4.1.0)
#>  farver        2.1.0   2021-02-28 [2] CRAN (R 4.0.5)
#>  fastmap       1.1.0   2021-01-25 [2] CRAN (R 4.0.3)
#>  fs            1.5.0   2020-07-31 [2] CRAN (R 4.0.2)
#>  generics      0.1.0   2020-10-31 [2] CRAN (R 4.0.3)
#>  ggplot2     * 3.3.5   2021-06-25 [2] CRAN (R 4.1.0)
#>  glue          1.4.2   2020-08-27 [2] CRAN (R 4.0.2)
#>  gtable        0.3.0   2019-03-25 [2] CRAN (R 4.0.2)
#>  highr         0.9     2021-04-16 [2] CRAN (R 4.0.5)
#>  htmltools     0.5.2   2021-08-25 [2] CRAN (R 4.1.1)
#>  httr          1.4.2   2020-07-20 [2] CRAN (R 4.0.2)
#>  KernSmooth    2.23-20 2021-05-03 [4] CRAN (R 4.0.5)
#>  knitr         1.36    2021-09-29 [2] CRAN (R 4.1.1)
#>  labeling      0.4.2   2020-10-20 [2] CRAN (R 4.0.3)
#>  lifecycle     1.0.1   2021-09-24 [2] CRAN (R 4.1.1)
#>  lwgeom        0.2-8   2021-10-06 [2] CRAN (R 4.1.1)
#>  magrittr      2.0.1   2020-11-17 [2] CRAN (R 4.0.3)
#>  mime          0.12    2021-09-28 [2] CRAN (R 4.1.1)
#>  munsell       0.5.0   2018-06-12 [2] CRAN (R 4.0.2)
#>  pillar        1.6.4   2021-10-18 [2] CRAN (R 4.1.1)
#>  pkgconfig     2.0.3   2019-09-22 [2] CRAN (R 4.0.2)
#>  proxy         0.4-26  2021-06-07 [2] CRAN (R 4.1.0)
#>  purrr         0.3.4   2020-04-17 [2] CRAN (R 4.0.2)
#>  R.cache       0.15.0  2021-04-30 [2] CRAN (R 4.0.5)
#>  R.methodsS3   1.8.1   2020-08-26 [2] CRAN (R 4.0.3)
#>  R.oo          1.24.0  2020-08-26 [2] CRAN (R 4.0.3)
#>  R.utils       2.11.0  2021-09-26 [2] CRAN (R 4.1.1)
#>  R6            2.5.1   2021-08-19 [2] CRAN (R 4.1.1)
#>  Rcpp          1.0.7   2021-07-07 [2] CRAN (R 4.1.0)
#>  reprex        2.0.1   2021-08-05 [2] CRAN (R 4.1.1)
#>  rlang         0.4.12  2021-10-18 [2] CRAN (R 4.1.1)
#>  rmarkdown     2.11    2021-09-14 [2] CRAN (R 4.1.1)
#>  rstudioapi    0.13    2020-11-12 [2] CRAN (R 4.0.3)
#>  s2            1.0.7   2021-09-28 [2] CRAN (R 4.1.1)
#>  scales        1.1.1   2020-05-11 [2] CRAN (R 4.0.2)
#>  sessioninfo   1.1.1   2018-11-05 [2] CRAN (R 4.0.5)
#>  sf          * 1.0-3   2021-10-07 [2] CRAN (R 4.1.1)
#>  stars       * 0.5-3   2021-06-08 [2] CRAN (R 4.1.1)
#>  stringi       1.7.5   2021-10-04 [2] CRAN (R 4.1.1)
#>  stringr       1.4.0   2019-02-10 [2] CRAN (R 4.0.2)
#>  styler        1.6.2   2021-09-23 [2] CRAN (R 4.1.1)
#>  tibble        3.1.5   2021-09-30 [2] CRAN (R 4.1.1)
#>  tidyselect    1.1.1   2021-04-30 [2] CRAN (R 4.0.5)
#>  units         0.7-2   2021-06-08 [2] CRAN (R 4.1.0)
#>  utf8          1.2.2   2021-07-24 [2] CRAN (R 4.1.0)
#>  vctrs         0.3.8   2021-04-29 [2] CRAN (R 4.0.5)
#>  withr         2.4.2   2021-04-18 [2] CRAN (R 4.0.5)
#>  wk            0.5.0   2021-07-13 [2] CRAN (R 4.1.0)
#>  xfun          0.27    2021-10-18 [2] CRAN (R 4.1.1)
#>  xml2          1.3.2   2020-04-23 [2] CRAN (R 4.0.2)
#>  yaml          2.2.1   2020-02-01 [2] CRAN (R 4.0.2)
#> 
#> [1] /home/gabo/R/x86_64-pc-linux-gnu-library/4.1
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library

[SimonDedman]

Thanks for this Gabo. This does indeed work, and I'm closing it as such. However, trying to run it with ggmap instead of ggplot results in:

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Error in FUN(X[[i]], ...) : object 'lon' not found

I presume the error is related to:

dkahle/ggmap#160 and dkahle/ggmap#162, however I've tried the solutions to these, with no success. I can't work out what the CRS of the bbox of myMap is. Or what the crs of myMap is in general. But I get the "Coordinate system already present" warning with the geom_stars line alone:

geom_stars(data = lemons.ras %>% st_transform(4326), inherit.aes = FALSE) +

But not when not transforming:

geom_stars(data = lemons.ras, inherit.aes = FALSE) +

suggesting the ggmap object myMap crs is the same as lemons.ras i.e. 2958. But

geom_stars(data = lemons.ras %>% st_transform(2958), inherit.aes = FALSE) +

Gives the "coordinate system..." warning, and plots, but plots nothing. This seems super weird given st_crs(lemons.ras) returns EPSG:2958, i.e. it's already in that CRS... Somewhat regardless, st_crs(mycontour) gives "EPSG:2958" as well, but it still gives the "coordinate system" warning and also the "object 'lon' not found" error.

I suspect I'll just have to add this onto the existing ggmap issue above, but I figured I'd share it with you in case you had any experience/ideas.

Cheers!