This vignette shows how stars object can be moved from vector and raster representations.

Rasterizing an sf vector object

The algorithm used is the GDAL rasterize utility, all options of this utility can be passed to st_rasterize. The geometry of the final raster can be controlled by passing a target bounding box and either the raster dimensions nx and ny, or pixel size by the deltax and deltay parameters.

Vectorizing a raster object to an sf object

stars objects can be converted into an sf object using st_as_sf. It has a number of options, depending on whether pixels represent the point value at the pixel center, or small square polygons with a single value.

We will work again with the landsat-7 6-band image, but will select the first band and round the values:

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)

Polygonizing

In case raster cells reflect point values and we want to get a vector representatoin of the whole field, we can draw contour lines and export the contour sets (only available when the GDAL version is at least 2.4.0):

l =  st_contour(x,contour_lines=TRUE, breaks = 8:14)
plot(l, key.pos = 1, pal = sf.colors(6), lwd = 2, key.length = 0.8)

Exporting to points

Alternatively, we can simply export all the pixels as points, and get them either as a wide table with all bands per point, and no replicated POINT geometries:

or as a long table with a single attribute and all points replicated:

as we can see, an additional attribute band now indicates which band is concerned.

Exporting to polygons

Alternatively, we can export to polygons and either get a single polygon per pixel, as in

or merge polygons that have identical pixel values;

p = st_as_sf(x, as_points = FALSE, merge = TRUE)

When plotted with boundaries, we see the resolved boundaries of areas with the same pixel value:

plot(p)

A further option connect8 can be set to TRUE to use 8 connectedness, rather than the default 4 connectedness algorithm. In both cases, the polygons returned will often be invalid according to the simple feature standard, but can be made valid using lwgeom::st_make_valid.

Reprojecting a raster

If we accept that curvilinear rasters are rasters too, and that regular and rectilinear grids are special cases of curvilinear grids, reprojecting a raster is no longer a “problem”, it just recomputes new coordinates for every raster cell, and generally results in a curvilinear grid (that sometimes can be brought back to a regular or rectilinear grid). If curvilinear grid cells are represented by coordinates of the cell center, the actual shape of a grid cell gets lost, and this may be a larger effect if grid cells are large or if the transformation is stronger non-linear.

An example of the reprojection of the grid created above is

where it should be noted that the dimensionality of the grid didn’t change: the same set of raster cells has been replotted in the new CRS, but now in a curvilinear grid.

Warping a raster

Warping a raster means creating a new regular grid in a new CRS, based on a (usually regular) grid in another CRS. We can do the transformation of the previous section by first creating a target grid:

nc %>% st_transform("+proj=laea +lat_0=34 +lon_0=-60") %>% st_bbox() %>%
    st_as_stars() -> newgrid

and then warping the old raster to the new

This new object has a regular grid in the new CRS, aligned with the new x- and y-axes.