I am using the R package landscapemetrics to find forest patch centroids in a landscape. The given function to find patch centroids is get_centroids(). The centroids returned from get_centroids() are centroids of patches that are too small and/or not connected to other cells of the same class. This results in many centroid dots that may bias the results of my later analysis. See the reprex below:

Loading libraries and downloading data:

#> Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
#> Loading required package: sp

# Download and read shapefile
countyurl = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip"
  download.file(countyurl,destfile = file.path(tdir,"Counties.zip"))

Albany_county = read_sf(paste(tdir,"/cb_2018_us_county_500k.shp",sep="")) %>%
         NAME=="Albany") %>%

# Get landcover data
LC = get_nlcd(template=Albany_county,
              landmass = 'L48',
              force.redo = T,
              extraction.dir = tdir)

Recoding landcover classes to include all forest types:

values(LC)[values(LC)==42] = 41
values(LC)[values(LC)==43] = 41
LC2 = LC
values(LC2)[values(LC2)!=41] = NA

Calculating patch centroids:

centroids = get_centroids(LC2,
                          return_sp = T,
                          directions = 8,
                          cell_center = F) %>%

Map generation, focusing on examples of points in question:

         bbox = st_bbox(centroids %>%
                        filter(id %in% c(4743,
#> stars object downsampled to 841 by 1189 cells. See tm_shape manual (argument raster.downsample)

Clearly, the points in the center of the map are not located in or on large forest patches. I have seen the option of conducting a negative buffer on forest patches to remove these small patches, thus removing their centroids. The problem with this solution is that the get_centroids() function requires the landscape argument, which accepts the full raster dataset rather than the results from the separate get_patches() function. I.e., using that workflow, I cannot conduct a negative buffer on get_patches() to then supply to get_centroids().

  • What if you explicitly use the patches() function from terra to create patches before running get_centroids()? For example, p1 <- patches(LC2, directions=8, allowGaps=F); centroids = get_centroids(p1, return_sp = T, directions = 8, cell_center = F) %>% st_as_sf()
    – John Polo
    Commented Oct 19, 2022 at 1:08
  • @JohnPolo: Wouldn't you calculate centroids of patches only this way? I'd imagine this is exactly what you want to avoid.
    – dimfalk
    Commented Oct 19, 2022 at 10:46
  • I think I misunderstand @falk-env's comment. They seem to argue that my idea was the wrong approach, but then write an answer using it anyway? lol. Like I said, maybe I misunderstand. Anyway, I wonder if, with either original approach or falk-env's answer, there is an ID, or a way to make one, for each "patch" ['you keep using that word...' :D ] and for the centroids that are created. Then you can try to search and match to see if the centroids and patches really are incorrect or misaligned. May need to write the objects to file and use something like QGIS, for interactive label and search.
    – John Polo
    Commented Oct 19, 2022 at 13:27
  • ... if you can see where they are incorrectly ending up, that may give you an indication of what's going wrong.
    – John Polo
    Commented Oct 19, 2022 at 13:31
  • 1
    Sorry if I was being imprecise. My answer must look a little bit weird then, I agree :) What I meant to say - and you're right, I had a different understanding of "patches" when I wrote the upper comment and my answer - using your approach, you would also create centroids for the smallest patches (and this is what you want to avoid, at least as fast as I understand). What I actually meant in the comment above was "clutter" - undesired noise. So one could argue to calculate centroids for all patches and filter points retrospectively - or filter the data basis (as I did) and calculate centroids.
    – dimfalk
    Commented Oct 19, 2022 at 16:21

3 Answers 3


I would agree that the OP should move this analysis to terra. I tested much the same as this answer and came to exactly the same results but, gave pause in answering because, I will say that this get_centroids function seems a bit "hinky" in what it represents in regard to patch centrality. Please note that using terra you can reclassify your landcover classes by simply using: ifel(LC %in% c(41,42,43), 1, NA). I will report this bug to the landscapemetrics developers as this function should accept a SpatRaster object but, throws an error. We should be moving our analysis/functions to sf (vector) and terra (raster) because sp, rgeos, rdgal, and maptools are being retired at the beginning of the year.

  • Great, thanks a lot for reporting!
    – dimfalk
    Commented Oct 20, 2022 at 21:54
  • Or: classify(LC, cbind(c(41,42,43), 1), others=NA) Commented Oct 21, 2022 at 19:57
  • Posted this as a comment on my phone and it turned up as an answer, sorry. Thanks @RobertHijmans for the alternative using classify. Commented Oct 21, 2022 at 20:38

Not 100 % sure about this, so better check the example section at ?terra::patches yourself also, section "use zonal to remove small patches" to be precise. This seems to be the recommended sieve approach at the moment.

I have to admit, this whole conversion of spatial objects (sp/sf/SpatVector; raster/stars/SpatRaster) is quite annoying here. Your raster objects were created using {raster}, now we need a {terra} function. Should be trivial, and landscapemetrics::get_centroid() should in fact also accept SpatRaster objects as landscape argument. Perfect. But it does not and fails with some error message about the proj4string. But maybe I messed up, so better try yourself with landscape = rast(LC2).

However, the given example can be adjusted to match your code:

#> terra 1.6.17

# create a SpatRaster object to access terra functions
LC3 <- rast(LC2)

# create patches
y <- patches(LC3)

# remove patches smaller than 10 ha, save output as raster again
rz <- zonal(cellSize(y, unit = "ha"), y, sum, as.raster = TRUE)
LC4 <- ifel(rz < 10, NA, LC3) |> raster()

# calculate patch centroids after filtering by area
centroids <- get_centroids(LC4,
                           return_sp = TRUE,
                           directions = 8,
                           cell_center = TRUE) |>

# compare cells in LC2 with centroids based on LC4
# your dplyr::filter() hack does not work anymore due to changes in cell IDs
         bbox = st_bbox(centroids)) +
  tm_raster() +
  tm_shape(centroids) +

As you can see, several rather small patches in LC2 do not have a corresponding centroid anymore. Is this what you were after?

Pretty sure you should tweak the chosen threshold value in ifel() a little, since I just chose a random (but maybe plausible) value of 10 ha.

  • This works, thank you! The differences between the terra and raster package are definitely frustrating. Commented Oct 20, 2022 at 21:32

#> stars object downsampled to 841 by 1189 cells. See tm_shape manual (argument raster.downsample) -- that is the source of your confusion. The get_centroids() function works fine.

Try to use raster.downsample = FALSE and then you will see:

enter image description here

  • Wow. Great find, thank you! I think one of the goals too was to remove the single cell patches, which the answers above address. Commented Oct 20, 2022 at 22:13

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.