Stars: fasterize compatibility

[noamross]

I envisioned fasterize as a link between sf and raster data types. Only one (one-way) operation is currently implemented but I had planned a fast version of raster2polygons, and potentially others. I had just been thinking about outputs besides in-memory rasters, and would be happy to try to implement stars-style input and output formats.

[edzer]

Sounds like a great idea, thanks!

[edzer]

Now that fasterize is on CRAN and stars is close, how shall we do this?

[noamross]

Yes! But FYI, I don't anticipate being able to spend significantly more time on fasterize until the summer.

I don't anticipate that it would require anything in stars - I just have to learn the data model. Fasterize basically creates an empty raster object, then the C++ code writes the data slot for the raster. So I just need a method to create an empty stars object, transfer over appropriate CRS projection ino, and then access the data, which is even easier because its an S3 object.

So, my todos:

• Add code to transfer CRS and other appropriate metadata from the sf object to an empty stars object. From a quick look this is basically st_as_star.bbox.
• Break up the rasterize() core function to separate out rasterization logic from raster/stars container-creation logic, and just pass a pointer to the data slot of the raster/stars object the rasterization function.
• Add the logic to update the stars object metadata after rasterizing (the equiv. of this. I guess this is basically adding the metadata for the third(+) dimension(s), which you don't know before rasterization because the # of unique layers is determined as you go.

So it's pretty simple, as long as we're just outputting an array. I don't anticipate fasterize() outputting a list of multiple arrays for a more complex stars object - do you think it should?

Also, would you want to import fasterize for functions like st_as_stars.sf? (and the reverse when I get polygonize() working?). I would argue against it - it makes maintenance easier to keep them separate.

As I said, I'm crunched for time at the moment. I will put this in the roadmap and issues in case someone wants to jump ahead of me, if not I _will_ get to it, just in a few months.

[edzer]

Much simpler might be to modify the R code of fasterize, to

• check if argument raster is of class stars, and if yes
• convert it to Raster
• do the usual stuff
• convert it back to stars

As long as all of this is in memory, there should be negligible overhead. Or am I overlooking something?

stars has an st_as_stars.sfc method to rasterize feature geometries, I think I'm going to use fasterize by default, for that.

Ah, I now see that you probably don't have any R code in fasterize, yes, that complicates.

[adrfantini]

Did you consider using GDAL_Polygonize and GDAL_FPolygonize, too?
Currently that is the fastest way to compute raster -> polygon transformation, followed by spex::polygonize. IIRC, @edzer, you spoke about adding them to sf::gdal_utils?

Note that going through raster might not always be the best idea, since not all stars object will play nice with raster (e.g. the curvilinear grids you just implemented).

[edzer]

... and it's sibling, GDALRasterizeGeometries. Thanks, good sugestion!

It's not one of the GDAL utils (different interface), but definitely in the scope of sf/stars. I guess it should go into the sf_as_stars.sf and st_as_sf.stars methods?

[adrfantini]

@edzer yeah, I think that would be ideal, very intuitive to put them there as optional arguments.

[mdsumner]

I got a fair way into a standard gdalwarp in Rcpp, writing to a GeoTiFF - I guess a good chunk here would be a stars/sf to MEM and MEM to stars capability and we could expose many gdal functions via that. Have you tried that @edzer ?

[edzer]

No haven't tried - good idea!

[edzer]

In the new logic, st_as_stars.sf creates a vector data cube, and doesn't rasterize; st_rasterize rasterizes, using the gdal_rasterize utility, going through temp files, and allowing all the options that has. There's a vignette explaining this.

I guess we can close here; please reopen (or comment) if one of you disagrees.