Stars: raster <-> stars migration table

[adrfantini]

Here's a first attempt at the raster <-> stars migration table we talked about in #112. Many things are still to define, and I have added my comments on what could or should be implemented (not necessarily in stars) in my humble opinion.
This migration table is currently very rough and many things are for sure imprecise or plain wrong. Any comment or correction is hugely appreciated. In a future I might add some info on xarray and GDAL too, but I'll have to do some more research since I do not use them much.

I post this here so we can get some useful discussion before it goes on the wiki.

LEGEND

? = I am not sure / I do not know how to do that in stars

* = Not present, _could_ be implemented

# = Not present, _should_ be implemented (imho)

Creating objects

raster | stars
------------ | -------------
{raster, stack, brick} (read) | read_stars or read_stars(along = ...)
{stack, brick} (concatenate layers) | c or c(along = ...)
subset | {[ ] , slice, filter}
addLayer | c() or c(along = ...)
dropLayer | {[ ] , slice, filter}
unstack | combine lapply and {[ ] , slice, filter}

Changing spatial extent and/or resolution of objects

raster | stars
------------ | -------------
merge | c (only for adjacent objects) #
mosaic | ?
crop | filter
extend | #
trim | #
aggregate | # (aggregate.stars, WIP?)
disaggregate | #
resample | {st_transform, st_warp}
projectRaster | {st_transform, st_warp}
shift | #
flip | #
rotate | *
t | ?

Cell based computation

raster | stars
------------ | -------------
calc | st_apply
overlay | c(along = , ...) %>% st_apply(...)
cover | [ ] <-
mask | [ ] <-
cut | cut
subs |
reclassify | * (forcats::fct_recode?)
init | [ ] <-
stackApply | {[ ] , slice, filter} %>% st_apply
stackSelect |

Spatial contextual computation

raster | stars
------------ | -------------
distance | #
gridDistance | *
distanceFromPoints | #
direction | *
focal | *
localFun | *
boundaries | st_as_sf(as_points=FALSE, merge=TRUE, connect8=TRUE)
clump |
adjacent | *
area | st_area
terrain | #
Moran |

Model predictions

raster | stars
------------ | -------------
Predict | ?
Interpolate | ?

Data type conversion

raster | stars
------------ | -------------
rasterize | st_as_stars
rasterToPoints | st_as_sf(as_points=TRUE)
rasterToPolygons | st_as_sf(as_points=FALSE, ...)
rasterToContour | # (st_contour)
rasterFromXYZ |
rasterFromCells |

Summarizing

raster | stars
------------ | -------------
cellStats | st_apply
summary | summary(as.vector(. %>% pull))
freq | *
crosstab |
unique | unique(as.vector(. %>% pull))
zonal | *

Accessing values of objects

raster | stars
------------ | -------------
getValues | {pull, [[ ]]}
getValuesBlock | {[ ] , slice, filter} %>% pull
getValuesFocal | {[ ] , slice, filter} %>% pull
as.matrix | is present, but does weird stuff *
as.array | is present, but does weird stuff *
extract (by cell) | {[ ] , slice, filter}
extract (by poligon) | [sf_object]
extract (by point) | #
sampleRandom | *
sampleRegular | *
minValue |
maxValue |
setMinMax |

Plotting

raster | stars
------------ | -------------
plot | plot, geom_stars
plotRGB | ?
spplot |
image |
persp |
contour |
filledContour |
text |
hist |
barplot |
density |
pairs |
boxplot |

Getting and setting dimensions

raster | stars
------------ | -------------
ncol | st_dimensions(.)$x %>% with(to - from + 1) (I am not sure how to generally get the x dimension name - here it is x)
nrow | st_dimensions(.)$y %>% with(to - from + 1) (I am not sure how to generally get the y dimension name - here it is y)
ncell | multiply the above
res | st_dimensions - can also not be a constant in case of irregular grids
nlayers | st_dimensions(.)$time %>% with(to - from + 1) (there is no way to generally get the z dimension name - here it is time)
names | names
xres | st_dimensions - can also not be a constant in case of irregular grids
yres | st_dimensions - can also not be a constant in case of irregular grids
xmin | st_bbox[1]
xmax | st_bbox[3]
ymin | st_bbox[2]
ymax | st_bbox[4]
extent | st_bbox
origin |
projection | st_crs
isLonLat | st_crs(.) == st_crs(4326)
filename | not available by design
bandnr |
nbands |
compareRaster | *
NAvalue |

Computing row, column, cell numbers and coordinates

raster | stars
------------ | -------------
xFromCol | st_get_dimension_values(., 'x')[col] (I am not sure how to generally get the x dimension name - here it is x)
yFromRow | st_get_dimension_values(., 'y')[row] (I am not sure how to generally get the y dimension name - here it is y)
xFromCell | not available by design
yFromCell | not available by design
xyFromCell | not available by design
colFromX | *
rowFromY | *
rowColFromCell | not available by design
cellFromXY | not available by design
cellFromRowCol | not available by design
cellsFromExtent | not available by design
coordinates | st_coordinates
validCell | not available by design
validCol | col %>% between(st_dimensions(.)$x$from, st_dimensions(.)$x$to) (I am not sure how to generally get the x dimension name - here it is x)
validRow | row %>% between(st_dimensions(.)$y$from, st_dimensions(.)$y$to) (I am not sure how to generally get the y dimension name - here it is y)
setValues | [ ] <-
writeRaster | write_stars (currently uses GDAL, somewhat limited)
KML |

The format of this table follows the raster-package entry in the raster manual.

[tim-salabim]

This looks extremely helpful!

[edzer]

How do we go about this: do you want us to edit the tables in your comment above?

[Nowosad]

Great work @adrfantini. I might suggest adding a third column with notes. This way is would be easier to use the table - "previous function", "new function", "note/comment".
Is there a plan to move it to a wiki?

[edzer]

Started one here: https://github.com/r-spatial/stars/wiki/raster---stars-comparison-table

[LDalby]

This is really helpful, thanks!
Would also be cool if one could see if a raster functionality not currently implemented in stars was planned - so a _planned_ column maybe?

[mdsumner]

Wished, maybe ;)

[LDalby]

and in that case @mdsumner, I wish there was an _extract (by points)_ method as stars[sf_object] ;-)

[adrfantini]

Thank you for the kind words and suggestions. @edzer what do you think, we keep this discussion here open and move the table fully to the wiki (where few people can edit it), or do we set up a more collaborative effort somewhere else?

and in that case @mdsumner, I wish there was an extract (by points) method as stars[sf_object] ;-)

I've been longing for that too recently tbh, it would have helped in a couple of cases.

[edzer]

In the sf realm extract is an st_intersection between an sf(c) and a stars object. It would return a vector data cube in case the stars object has more than 2 dimensions.

[LDalby]

Right, I see.
Well, I guess that would actually be helpful still. In the 2D case it would be possible to have the stars values at the point locations as a column in the sf-point object more or less directly (which at least for me is what I most often want from this kind of operation). In the case of a stars object with more than 2 dimensions, I guess you would often aim for a similar thing, that is to have the dimensions of the vector data cube assigned as columns in the sf-point object.
I guess it would not be too many extra steps for the user to go from the data cube to the desired format if the method returned two different formats depending on the number of dimensions of the stars object.

[adrfantini]

In the sf realm extract is an st_intersection between an sf(c) and a stars object. It would return a vector data cube in case the stars object has more than 2 dimensions.

You mean something like this?

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x1 = read_stars(tif)
d = data.frame(a = 1:2)
d$geom = c("POINT(295000 9115000)", "POINT(296000 9116000)")
df = st_as_sf(d, wkt = 'geom', crs = st_crs(x1))
st_intersection(df,  st_xy2sfc(x1, as_points = FALSE))

Does not seem to return data from the stars object to me.

[edzer]

@adrfantini obviously not, because st_xy2sfc gives an sfc: geometries without attributes. But do try

st_intersection(df,  st_as_sf(x1, as_points = FALSE))

as your last command.

[edzer]

(A problem is that this is very expensive solution because it ignores the grid topology; for regular/rectilinear/rotate/shear grids a cheap approach is not hard to implement but remains TBD)

[adrfantini]

@adrfantini obviously not, because st_xy2sfc gives an sfc: geometries without attributes. But do try

st_intersection(df,  st_as_sf(x1, as_points = FALSE))

as your last command.

Oh I see, but then that's a sf-only solution, I was thinking of something directly with stars.

Anyway, what do you suggest we do with the above table? Obviously a place where there is some kind of collaborative editing would be handy.

[edzer]

Yes, we need a stars solution for this.

Re: table: move it to the wiki and we'll give edit rights to anyone who wants it.

[edzer]

Interestingly,

st_join(df,  st_as_sf(x1, as_points = FALSE))

gives you the same table as st_intersection().

[edzer]

The stars version:

> aggregate(x1, df, function(x) x[1], as_points = FALSE)
stars object with 2 dimensions and 1 attribute
attribute(s):
  L7_ETMs.tif   
 Min.   :46.00  
 1st Qu.:49.00  
 Median :63.00  
 Mean   :64.08  
 3rd Qu.:76.00  
 Max.   :88.00  
dimension(s):
         from to offset delta                       refsys point
geometry    1  2     NA    NA +proj=utm +zone=25 +south...  TRUE
band        1  6     NA    NA                           NA    NA
                                                 values
geometry POINT (295000 9115000), POINT (296000 9116000)
band                                               NULL

[adrfantini]

A mockup is up now:
https://github.com/r-spatial/stars/wiki/raster-stars-comparison-table

Please give it a read and edit (once you have editing rights) what you feel is incorrect or misleading. Thanks.

[edzer]

Great! Some edits here.

[adrfantini]

Great! Some edits here.

Thank you, excellent changes. Some things I blatantly missed (e.g. forgetting st_crop... which was added after an issue I created myself! Getting old...) are now present.

Looks great.

[edzer]

As I tweeted, the table is a raster index with suggestions how to do things with stars, not a comparison - things you can do with stars but not with raster are missing. I refrained from calling it a raster to stars migration table, as I don't want to suggest people should make that migration. But as it is now is also partially misleading. Any ideas? Add a table with features in stars, not in raster?

[adrfantini]

I find it hard to convey why one should use stars over raster. It's not about individual functions (raster is much more complete), but more about design differences. stars works well with sf, deals well with projections, handles multiple variables at a time, supports more than 3 dimensions, is more modern in design (dplyr support etc.) and more promising in the long run.

On the other hand, raster too has some advantages. For example, I personally find stars's proxy objects more cumbersome and less mature than the (completely transparent) raster counterpart.

These differences cannot be conveyed well in a table, I believe.

[Nowosad]

@edzer I think this table is good as it is, because it is useful for two groups:

1. Users - people who use raster and maybe want to try stars, or somebody who uses stars but used raster in the past.
2. Developers - to see what more can be implemented, or where the implementation could be improved (e.g. extract).

This table shows stars in the context of raster data.
Probably, it would be nice to have other tables, for example with stars functions in a context of spatiotemporal data, etc.

[edzer]

I also think the table is good, but the title isn't. I'll change it into "How raster functions map to stars functions".

[statnmap]

Are there any plan to update the Wiki page ? https://github.com/r-spatial/stars/wiki/How-%60raster%60-functions-map-to-%60stars%60-functions
At least with some links to issues to refer to like https://github.com/r-spatial/stars/issues/176
Also, wouldn't this table be better placed in a vignette in the package? At least, this would allow other dev to help update the content.

[edzer]

Agreed - I'd be happy to give you write access to the wiki, let me know; a vignette is also a great idea!