r-spatial
Different outputs from stars::st_contour() and graphics::contour in R
I have a raster object (link to file here) that plots contours 'wrong' in ggplot2 compared to base graphics.
Using ggplot2, the bottom area is missing the two 50% contour circles that are present in the base graphics plot. The total areas produced in the base graphics plot align to the calculated volume area sizes produced by the move::getVolumeUD(), which I'm inclined to trust/believe, and also align with the output generated using the fishtrack3d package: in terms of both calculated volume areas and UD contours.
EDIT: 2 packages were not loaded; hence the could would not run. Fixed now.
library(ggplot2)
library(stars)
library(magrittr)
library(move)
test <- stars::read_stars(file.path("All_Rasters_Scaled_Weighted.asc")) %>% sf::st_set_crs(4326)
test %<>% starsExtra::trim2()
test <- as(test, "Raster")
test <- new(".UD", test)
is.na(test[[1]]) <- test[[1]] == 0
plot(test, xlab="location_long", ylab="location_lat", asp = 1)
contour(test, levels=c(.5, .95), col=c(6,2), add=TRUE, lwd=2)
That produces the following plot:
Obviously the asp ratio is different b/w contours and plot and looks ugly, but the key point is that it produces the 2 bottom 50% contour areas. This is like the plot produced by the fishtrack3d R package shown below:
Below is the ggplot code:
test <- stars::read_stars(file.path("All_Rasters_Scaled_Weighted.asc")) %>% sf::st_set_crs(4326)
contour1colour = "red" # colour for contour 1, typically 95%.
contour2colour = "orange"
legendtitle = "Percent UD Contours"
ggplot() +
geom_stars(data = test) +
ggplot2::geom_sf(data = stars::st_contour(x = test, contour_lines = TRUE, breaks = max(test[[1]], na.rm = TRUE) * 0.05), fill = NA, inherit.aes = FALSE,
ggplot2::aes(colour = "95% UD")) +
ggplot2::geom_sf(data = stars::st_contour(x = test, contour_lines = TRUE, breaks = max(test[[1]], na.rm = TRUE) * 0.5), fill = NA, inherit.aes = FALSE, ggplot2::aes(colour = "50% UD")) +
ggplot2::scale_colour_manual(name = legendtitle, values = c("50% UD" = contour2colour, "95% UD" = contour1colour))
which produces the following 'ugly' plot:
In this plot you can see that the bottom contours are absent. Given that the object is the same in both approaches, can anyone propose why the contour outputs are different?
Possibly relevant: in the ggplot2 approach we used 0.5 and 0.05 proportion of the object's max value as breaks values in st_contour(). In the base R approach, we set levels at 0.5 and 0.95. HOWEVER, the max value of this object is 0.008. According to the manual page for graphics::contour(), levels values seem like they should be absolute, rather than relative. Yet they behave as relative... or at least plot near identically to what we'd expect.
A SO thread reporting this issue exists here.
Thanks in advance for any help.
Please provide a working reprex:
> library(ggplot2)
> library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
> test <- stars::read_stars(file.path("All_Rasters_Scaled_Weighted.asc")) %>% sf::st_set_crs(4326)
> test %<>% starsExtra::trim2()
Error in test %<>% starsExtra::trim2() : could not find function "%<>%"
> test <- as(test, "Raster")
Warning message:
multiple methods tables found for ‘area’
> test <- new(".UD", test)
Error in getClass(Class, where = topenv(parent.frame())) :
“.UD” is not a defined class
> is.na(test[[1]]) <- test[[1]] == 0
> plot(test, xlab="location_long", ylab="location_lat", asp = 1)
Error in as.double(y) :
cannot coerce type 'S4' to vector of type 'double'
> contour(test, levels=c(.5, .95), col=c(6,2), add=TRUE, lwd=2)
Error in contour.default(test, levels = c(0.5, 0.95), col = c(6, 2), add = TRUE, :
no proper 'z' matrix specified
It's not a surprise that contour() doesn't work after the plot before it, as it doesn't know the plot settings. It should (only) work after a call to graphics::image().
My apologies, the provided code was missing the loading of two packages, which have now been included in the OP. The example should now be reproducible. Do you have any thoughts regarding the discrepancy? Thank you.
As I mentioned, the first plot is a matter of two base plot methods (plot and contour) that make different assumptions, so don't match, but also have nothing to do with stars. The fishtrack3d vs the ggplot2 plot seems a matter of aspect ratio; for longlata not on the equator both approaches seem to choose a different value for this. Which one do you find appropriate?
Hi Edzer,
The base plot methods were shown simply to demonstrate that these approaches take the same input object and both (plot/contour & fishtrack3d) produce an output that shows two 50%UD (inner) contours inside the bottom 'blob' of 95%UD contour. For all three plots, the aspect ratios are a mess, the colours are janky or absent, but this is not our concern - these are just quick & dirty outputs to demonstrate the titular issue.
Given stars::st_contour is using the same input file as graphics::contour & fishtrack3d (read in by stars::read_stars for all), one would expect the same contours for each of them, no? I.e.: there should be two yellow circles inside the bottom red blob in the third plot. We were hoping you'd be able to reproduce this discrepancy, and be similarly surprised that the ggplot2+stars::st_contour approach doesn't match the graphics::contour and fishtrack3d outputs.
Hi @edzer , hope you're well. I don't suppose you've had a few moments to look into this? With this issue, and an unrelated problem whereby converting a raster from raster::raster using stars::st_as_stars possibly loses key information, we are unable to advance our project any further.
Many thanks in advance.
Given stars::st_contour is using the same input file as graphics::contour & fishtrack3d (read in by stars::read_stars for all), one would expect the same contours for each of them, no?
I don't - both use different algorithms and come from different code bases with different history. st_contour uses a GDAL algorithm, base::contour uses an R internal algorithm. I'm not a specialist on contour algorithms.
With this issue, and an unrelated problem whereby converting a raster from raster::raster using stars::st_as_stars possibly loses key information
I put a lot of effort in stars <--> terra compatibility, I think raster is legacy and you should consider converting raster to terra (or migrate) and then look at conversion from that to stars. Happy to look into issues with terra.
Thanks Edzer. The intel on the underlying differences in algorithms is very interesting. In my tests I found that the st_contour / GDAL implementation works correctly, while the graphics::contour implementation doesn't.
Re: reading in a stars object and assigning a CRS, we solved this, though the successful approach was unintuitive and surprising.