I have an sf-object (shp downloadable here: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT) with multiple glaciers.

I want to delete glaciers smaller than 0.05km² by using a for-loop. I think, I just got something wrong in the last loop, but to be reproducable: my general code so far is

glaciers <- readOGR("C:/***", "Glacier_Clip_Polygon")
glaciers_sf <- st_as_sf(glaciers)
glaciers_comb <- st_union(glaciers_sf) # some glaciers are separated into different polygons so I merged them
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")

Then I access the area of the different polygons and store them in a vector to be used in a loop.

area <- st_area(glaciers_comb)
area_v <- as.vector(area)

area_cutoff <- vector("numeric")

for(i in 1:length(area_v)) {
    area_cutoff[i] <- area_v[i]

area_cutoff[is.na(area_cutoff)] <- 0

Now I want R to loop through this vector and only save the big polygons to the original object.

for(i in 1:length(area_cutoff)){
  glaciers_comb[i] <- glaciers_comb[i]

However, R glaciers_combdoes not change when I run this loop. Where did I go wrong?


Because in this:

for(i in 1:length(area_cutoff)){
   glaciers_comb[i] <- glaciers_comb[i]

If the if condition is true you are replacing element i with element i - which is no replacement. And when the if condition isn't true then element i is unchanged. So nothing changes.

You can do what you want with a one-liner:

glaciers_large <- glaciers_comb[as.numeric(st_area(glaciers_comb)) > 50000]

leaving 27 features:

> glaciers_large
Geometry set for 27 features 
geometry type:  POLYGON

This is basic R use of square brackets to subset elements of a vector by a true or false condition. The as.numeric is only needed because the area has units attached and otherwise you get:

> glaciers_large <- glaciers_comb[st_area(glaciers_comb) > 50000]
Error in Ops.units(st_area(glaciers_comb), 50000) : 
  both operands of the expression should be "units" objects

So really you should compare against a size in metres^2 using the units package:

> glaciers_large <- glaciers_comb[st_area(glaciers_comb) > units::set_units(50000,"m^2")]
  • Thanks a lot for the detailed answer! – Florian Mlehliv Jul 3 '19 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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