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

Rasterizing an sf vector object

library(stars)
## Loading required package: abind
system.file("gpkg/nc.gpkg", package = "sf") %>%
  read_sf() %>%
  st_transform(32119) -> nc
nc$dens = nc$BIR79 / units::set_units(st_area(nc), km^2)
(nc.st = st_rasterize(nc["dens"], dx = 5000, dy = 5000))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  dens [1/km^2]   
##  Min.   : 0.255  
##  1st Qu.: 1.226  
##  Median : 1.932  
##  Mean   : 3.346  
##  3rd Qu.: 3.826  
##  Max.   :21.248  
##  NA's   :4808    
## dimension(s):
##   from  to offset delta                       refsys point values    
## x    1 162 123829  5000 PROJCS["NAD83 / North Car... FALSE   NULL [x]
## y    1  61 318260 -5000 PROJCS["NAD83 / North Car... FALSE   NULL [y]
plot(nc.st)

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 dx and dy 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 representation 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:

st_as_sf(x, as_points = TRUE, merge = FALSE)
## Simple feature collection with 2500 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 288790.5 ymin: 9119350 xmax: 290187 ymax: 9120747
## CRS:            PROJCS["UTM Zone 25, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
## First 10 features:
##    L7_ETMs.tif.V1 L7_ETMs.tif.V2                 geometry
## 1              14             11 POINT (288790.5 9120747)
## 2              14             11   POINT (288819 9120747)
## 3              13             10 POINT (288847.5 9120747)
## 4              12              9   POINT (288876 9120747)
## 5              12             10 POINT (288904.5 9120747)
## 6              12             10   POINT (288933 9120747)
## 7              12             10 POINT (288961.5 9120747)
## 8              12             10   POINT (288990 9120747)
## 9              13             10 POINT (289018.5 9120747)
## 10             13             10   POINT (289047 9120747)

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

st_as_sf(x, as_points = TRUE, merge = FALSE, long = TRUE)
## Simple feature collection with 5000 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 288790.5 ymin: 9119350 xmax: 290187 ymax: 9120747
## CRS:            PROJCS["UTM Zone 25, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
## First 10 features:
##    band L7_ETMs.tif                 geometry
## 1     1          14 POINT (288790.5 9120747)
## 2     1          14   POINT (288819 9120747)
## 3     1          13 POINT (288847.5 9120747)
## 4     1          12   POINT (288876 9120747)
## 5     1          12 POINT (288904.5 9120747)
## 6     1          12   POINT (288933 9120747)
## 7     1          12 POINT (288961.5 9120747)
## 8     1          12   POINT (288990 9120747)
## 9     1          13 POINT (289018.5 9120747)
## 10    1          13   POINT (289047 9120747)

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

st_as_sf(x[1], as_points = FALSE, merge = FALSE)
## Simple feature collection with 2500 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 288776.3 ymin: 9119336 xmax: 290201.3 ymax: 9120761
## CRS:            PROJCS["UTM Zone 25, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
## First 10 features:
##    L7_ETMs.tif.V1 L7_ETMs.tif.V2                       geometry
## 1              14             11 POLYGON ((288776.3 9120761,...
## 2              14             11 POLYGON ((288804.8 9120761,...
## 3              13             10 POLYGON ((288833.3 9120761,...
## 4              12              9 POLYGON ((288861.8 9120761,...
## 5              12             10 POLYGON ((288890.3 9120761,...
## 6              12             10 POLYGON ((288918.8 9120761,...
## 7              12             10 POLYGON ((288947.3 9120761,...
## 8              12             10 POLYGON ((288975.8 9120761,...
## 9              13             10 POLYGON ((289004.3 9120761,...
## 10             13             10 POLYGON ((289032.8 9120761,...

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.

Switching between vector and raster in stars objects

We can convert a raster dimension to a vector dimension while keeping other dimensions as they are in a stars object by

x.sf = st_xy2sfc(x, as_points = TRUE)
x.sf
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##   L7_ETMs.tif   
##  Min.   : 7.00  
##  1st Qu.: 9.00  
##  Median :11.00  
##  Mean   :11.25  
##  3rd Qu.:12.00  
##  Max.   :28.00  
## dimension(s):
##          from   to offset delta                       refsys point
## geometry    1 2500     NA    NA PROJCS["UTM Zone 25, Sout...  TRUE
## band        1    2     NA    NA                           NA    NA
##                                                       values
## geometry POINT (288790.5 9120747),...,POINT (290187 9119350)
## band                                                    NULL

which also requires setting the as_points arguments as in st_as_sf.

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

nc.st %>% st_transform("+proj=laea +lat_0=34 +lon_0=-60") -> nc.curv
nc.curv
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  dens [1/km^2]   
##  Min.   : 0.255  
##  1st Qu.: 1.226  
##  Median : 1.932  
##  Mean   : 3.346  
##  3rd Qu.: 3.826  
##  Max.   :21.248  
##  NA's   :4808    
## dimension(s):
##   from  to offset delta                       refsys point
## x    1 162     NA    NA +proj=laea +lat_0=34 +lon... FALSE
## y    1  61     NA    NA +proj=laea +lat_0=34 +lon... FALSE
##                           values    
## x [162x61] -2210936,...,-1371611 [x]
## y    [162x61] 90650.2,...,538204 [y]
## curvilinear grid
plot(nc.curv, border = NA, graticule = TRUE)

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

nc.st %>% st_warp(newgrid) -> nc.new
nc.new
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  dens [1/km^2]  
##  Min.   : 0.25  
##  1st Qu.: 1.23  
##  Median : 1.93  
##  Mean   : 3.34  
##  3rd Qu.: 3.83  
##  Max.   :21.25  
##  NA's   :36155  
## dimension(s):
##   from  to   offset    delta                       refsys point values    
## x    1 380 -2188110  2098.09 +proj=laea +lat_0=34 +lon...    NA   NULL [x]
## y    1 171   494924 -2098.09 +proj=laea +lat_0=34 +lon...    NA   NULL [y]
plot(nc.new)

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