st_contour - mapping should be created with aes

Question

I know the "mapping should be created with aes" bit has been done to death, so I hope this is a sufficiently unusual case, since I've tried to work out why it's happening and come up short.

Per this thread, 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. The data file in question is All_Rasters_Scaled.asc.

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

Works fine. Edzer informed me that st_contour returns an sf object with the contours, and you have to pass that to geom_sf. So:

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

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. The object might have been incorrectly created somehow and have all NAs though:

st_bbox(mycontour)

xmin ymin xmax ymax

NA NA NA NA

But if I create it without any parameters specified:

mycontour <- st_contour(x = lemons.ras)

The bbox output is the same. Also the contour min/maxes in the object are oddly overlapping. Regardless of the way I create mycontour, I get the same aes error. Trying to specify aes:

geom_sf(mycontour, mapping = aes())

Coordinate system already present. Adding new coordinate system, which will replace the existing one. Error in FUN(X[[i]], ...) : object 'lon' not found

Or

geom_sf(mycontour, mapping = aes(geometry))

Warning: Ignoring unknown aesthetics: x Coordinate system already present. Adding new coordinate system, which will replace the existing one. Error in FUN(X[[i]], ...) : object 'lat' not found

Basically it looks like the ggmap / ggplot framework doesn't understand that mycontour, created by geom_contour, is an sf object, and that its aes mapping is the geometry column.

Any suggestions very much appreciated. I'm not sure if there's an underlying issue whereby geom_contour is failing to create the object properly and therefore all values are NA, and I don't know how to extract the geometry data from an sf object.

Code for the mymap object for ggmap is:

myLocation <- st_bbox(lemons.ras) %>% as.vector()
myMap <- get_map(location = myLocation)

But sf objects typically work fine with a blank ggplot()+

Thanks in advance!

Edit: adding more details:

class(mycontour)

"sf" "data.frame"

mycontour$geometry[1]

Geometry set for 1 feature

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 ...

More edits:

ggplot() +
geom_stars(lemons.ras)

Works fine, creates:

st_bbox(lemons.ras)

xmin ymin xmax ymax

-79.28970 25.68331 -79.21143 25.78939

i.e. looks fine also.

summary(lemons.ras)

                   Length Class  Mode   

All_Rasters_Scaled.asc 36270 -none- numeric

Including or not including x= in st_contour() shouldn't (and, having tested, doesn't) make a difference. I know it's unnecessary but I like doing it for two reasons:

1. It rarely hurts to be explicit, especially when sharing code & teaching others.

2. With RStudio's code completion macros, you can hit tab within a function and it'll give you the list of parameters to choose from. I try to use those in order. Also if you explicitly use x (or whatever the default first parameter is) then the macro will stop giving you that parameter option. All minor semantics though :)

lemons.ras$All_Rasters_Scaled.asc is a matrix arrray, 155x234 in size, with values ranging from 0 to 1 (it's been scaled, hence the name). So conceptually none of the values in mycontour from st_contour should be NA.

st_contour's param na.rm is default TRUE thus should be removing these NAs in any case, I feel.

I'm increasingly thinking that there's an issue with st_contour...

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

@lovelery 's code to test:

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)[, 1:50, 1:50, 1:2]
x[[1]] = round(x[[1]]/5)
l = st_contour(x, contour_lines = TRUE, breaks = 11:15)
plot(l[1], key.pos = 1, pal = sf.colors, lwd = 2, key.length = 0.8)

Works fine:

l is created with st_contour(x, contour_lines = TRUE, ...) which makes a multi linestring feature instead of a multipolygon feature. If I do the same to lemons.ras, the resulting mycontour has all NA values, which is probably related to the problem. l:

mycontour:

The files that created them with st_contour: x:

lemons.ras (note there are values here, from 0 to 1):

Differences:

• x$L7_ETMs.tif values are integers. But having upscaled lemons.ras values to integers, mycontour still fails.

• x$L7_ETMs.tif is 3 dimensions [50 x 50 x 2], lemons.ras$All_Rasters_Scales.asc is 2 dimensions [155 x 234]. But the code works with plot(l[1]) and plot(l) so the extra band/dimension shouldn't matter, right?!

Edit: if I remove the 1:2 from the x = read_stars, plot(l) works the same. BUT: x is now Type: list [50 x 50 x 6]. Before was [50 x 50 x 2]. Ok this is because the input file, L7_ETMs.tif, is [349 x 352 x 6] and the indexing read_stars(tif)[, 1:50, 1:50, 1:2] constrains this to 2 bands.

st_dimensions(x) gives:

st_dimensions(lemons.ras) gives:

So this might be the root of the issue. Even though lemons.ras plots perfectly fine having been imported as a raster and read into stars, it seems to be missing band attributes, which honestly don't look to be adding much of anything, but maybe the st_contour code expects them to be present and fails without them.

Actually looking at the direct output of st_dimensions is maybe even more illuminating: For x:

For lemons.ras:

Compared to x, lemons.ras has NA offset, delta and point. x has NULL values, so I guess the x/y geometry data for x is within the band, whereas for lemons.ras it's in the values. LAWD raster data is confusing.

Reading the stars intro again, it looks like splitting & merging attributes requires the band to be present in the first place. So the problem looks increasingly to be due to the lemons.ras input raster not having a band dimension.

The file was created with:

writeRaster(x, name, format = "ascii", datatype = "FLT4S", bylayer = T)

I don't know if the choice of ascii means no band info? I doubt it? After creation of the original raster earlier in the code, it has properties:

class : RasterLayer

dimensions : 3270, 3245, 10611150 (nrow, ncol, ncell)

resolution : 50, 50 (x, y)

extent : 594830.7, 757080.7, 2763971, 2927471 (xmin, xmax, ymin, ymax)

crs : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs

source : memory

names : layer

values : 1, 1 (min, max)

DURING creation of this file, it looks like band or layer parameters were NOT used. This is necessarily going to be a single layer/band raster, and evidently the output is successful enough that it survives multiple rounds of processing, editing, joining, applying, etc etc, before hitting this issue. Values were added after raster creation with setValues. And there is the "names : layer" entry above.

intermediate Conclusion this morning:

my input raster doesn't have a band attribute and is just a plain/flat raster, which is fine for simple plotting but evidently isn't the rich format that stars expects. However stars will work with this file and convert it to a (weak) stars object, but then it will fail in st_contour, likely due to the lack of the band.