sf::st_read() and terra::vect() reading in different extents from the same shapefile

Question

I am currently working with multiple shapefiles in R and working towards rasterizing them using a combination of functions from sf,terra and stars. However, I noticed while reading in a shapefile that terra::vect() and sf::st_read() are providing me with different extents for the same shapefile. Please see below.

 terra::vect(shapefile1)
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 538, 3  (geometries, attributes)
 extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :    id       class gridcode
 type        : <int>       <chr>    <int>
 values      :     1 settlements        2
                   2 settlements        2
                   3 settlements        2

sf::st_read(shapefile1)
Reading layer `1870-settlements' from data source 
  `data\spatial\1870-settlements.shp' 
  using driver `ESRI Shapefile'
replacing null geometries with empty geometries
Simple feature collection with 538 features and 3 fields (with 1 geometry empty)
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 76.492 ymin: 11.22461 xmax: 77.00837 ymax: 11.53395
Geodetic CRS:  WGS 84

Why am I getting this discrepancy? Such an error surely impacts future analyses.

Edit: Adding a link to the file - https://drive.google.com/drive/u/2/folders/1uDVQt1F64rUTfCXclwr-sM-oFjyJ2NwV

Answer 1

The incorrect extent shown by terra does not appear to affect anything of importance, but this bug has now been fixed in terra version 1.5-28 (currently the development version, which you may be able to install with install.packages('terra', repos='https://rspatial.r-universe.dev'))

Records that do not have an associated geometry now get an extent that is NA (as in sf).

library(terra)
# terra 1.5.28
v <- vect("1870-nilgiris-settlements.shp")
v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 538, 3  (geometries, attributes)
# extent      : 76.492, 77.00837, 11.22461, 11.53395  (xmin, xmax, ymin, ymax)
# source      : 1870-nilgiris-settlements.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :    id       class gridcode
# type        : <int>       <chr>    <int>
# values      :     1 settlements        2
#                   2 settlements        2
#                   3 settlements        2

i <- emptyGeoms(v)
i
#[1] 371
v[i]
# class       : SpatVector 
# geometry    : null 
# dimensions  : 1, 3  (geometries, attributes)
# extent      : NaN, NaN, NaN, NaN  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :    id       class gridcode
# type        : <int>       <chr>    <int>
# values      :   371 settlements        2

Answer 2

The 371st feature has a null geometry, which terra assigns the extent of "anywhere on the globe":

> v[371]
 class       : SpatVector 
 geometry    : unknown 
 dimensions  : 1, 3  (geometries, attributes)
 extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :    id       class gridcode
 type        : <int>       <chr>    <int>
 values      :   371 settlements        2

Whereas sf says "Its... no idea":

> s[371,]
Simple feature collection with 1 feature and 3 fields (with 1 geometry empty)
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
     id       class gridcode      geometry
371 371 settlements        2 POLYGON EMPTY

There's an argument for either approach, I guess. You could also argue that sf should return NA for the bounds of anything with at least one feature with NA bounds, but it doesn't:

> s[371:372,]
Simple feature collection with 2 features and 3 fields (with 1 geometry empty)
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 76.87192 ymin: 11.49424 xmax: 76.87287 ymax: 11.49516
Geodetic CRS:  WGS 84
     id       class gridcode                       geometry
371 371 settlements        2                  POLYGON EMPTY
372 372 settlements        2 POLYGON ((76.87235 11.49516...

Interesting different behaviour between packages though, might be worth reporting to the maintainers.

I don't know what the Simple Features standard says about bounds for missing geometries. In the end I guess you'll want to fix this point anyway, either dropping it from analyses or giving it a location.