3

I am currently struggling with rasterizing and later, resampling a vector object using functions from stars and terra.

Let's start with the example dataset here: https://drive.google.com/drive/u/2/folders/1uDVQt1F64rUTfCXclwr-sM-oFjyJ2NwV

This is a dataset of waterbodies and I load them first using sf

waterBodies_nil1870 <- st_read("data\spatial\waterBodies_nil1870.shp")
waterBodies_nil1870 <- st_transform(waterBodies_nil1870, 32643)

> waterBodies_nil1870
Simple feature collection with 5 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 667536.9 ymin: 1240976 xmax: 704538.2 ymax: 1276926
Projected CRS: WGS 84 / UTM zone 43N
  id       class gridcode                       geometry
1  1 waterbodies        1 MULTIPOLYGON (((683042.8 12...
2  2 waterbodies        1 MULTIPOLYGON (((669090.5 12...
3  3 waterbodies        1 MULTIPOLYGON (((675343.5 12...
4  4 waterbodies        1 MULTIPOLYGON (((702393.8 12...
5  5 waterbodies        1 MULTIPOLYGON (((688438.4 12...

Now I would like to rasterize this simple feature and resample it so that I further compute area statistics to compare with other rasters. To do so I have been struggling with using functions from both the stars and the terra package and each time, I believe I am losing values and data associated.

Let's start with functions from the stars package:

> new <- st_rasterize(waterBodies_nil1870)
> new
stars object with 2 dimensions and 2 attributes
attribute(s):
          Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
id           1       2      2 2.547945       3    5 64790
gridcode     1       1      1 1.000000       1    1 64790
dimension(s):
  from  to  offset    delta                refsys point values x/y
x    1 259  667537  143.276 WGS 84 / UTM zone 43N FALSE   NULL [x]
y    1 251 1276926 -143.276 WGS 84 / UTM zone 43N FALSE   NULL [y]

In the above scenario, I have now lost the class attribute column and only id and gridcode are preserved for some reason. Notice that this error continues to surface even if I do: st_rasterize(waterBodies_nil1870["class"]).

Let's try functions from the terra package now:

v1 <- vect(waterBodies_nil1870)
> v1
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 5, 3  (geometries, attributes)
 extent      : 667536.9, 704538.2, 1240976, 1276926  (xmin, xmax, ymin, ymax)
 coord. ref. : WGS 84 / UTM zone 43N (EPSG:32643) 
 names       :    id       class gridcode
 type        : <num>       <chr>    <num>
 values      :     1 waterbodies        1
                   2 waterbodies        1
                   3 waterbodies        1

> r1 <- rast(v1)
> r1
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 3700.129, 3595.036  (x, y)
extent      : 667536.9, 704538.2, 1240976, 1276926  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 43N (EPSG:32643) 
> hasValues(r1)
[1] FALSE

How do I rasterize a simple features object and then resample it using R?

1 Answer 1

3

With terra you can rasterize each field to r1 and merge all the objects onto a multi-layer raster:

library(sf)


waterBodies_nil1870 <- st_read("./stackExchange/1870-nilgiris-waterBodies.shp")
waterBodies_nil1870 <- st_transform(waterBodies_nil1870, 32643)

waterBodies_nil1870

# With terra
library(terra)

v1 <- vect(waterBodies_nil1870)
r1 <- rast(v1)

# Rasterize each value of v1
nams <- names(v1)

allrast <- lapply(nams, function(x) {
  rasterize(v1, r1,
    field = x,
    touches = TRUE
  )
})

# Merge (bind) all objects
allrast <- do.call("c", allrast)

allrast

#> class       : SpatRaster 
#> dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
#> resolution  : 3700.129, 3595.036  (x, y)
#> extent      : 667536.9, 704538.2, 1240976, 1276926  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 43N (EPSG:32643) 
#> sources     : memory  
#>               memory  
#>               memory  
#> names       : id,        class, gridcode 
#> min values  :  1, water bodies,        1 
#> max values  :  5, water bodies,        1 

The object allrast has all the columns of v1.

A couple of plots just to check:

plot(allrast)

enter image description here


plot(allrast, 1)
plot(v1, add=TRUE)

enter image description here

1
  • I recreated your code with my data. However, I get very blocky dimensions. The original SpatVector has dimensions : 575, 1 (geometries, attributes), but the SpatRaster gives: dimensions : 10, 10, 1 (nrow, ncol, nlyr), resolution : 25575.62, 19973.94 (x, y). How do you obtain or a set a desired resolution through reference to another raster? Jul 18, 2022 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.