Dividing Individual Spatial Polygons Equally in R

Question

I have a shapefile of polygons that are the townships in the state of Iowa.I'd like to divide each element (ie each township) into 9 equal parts (i.e. a 3 x 3 grid for each township). I've figured out how to do this, but am having trouble forming a new dataframe out of the new polygons. My code is below. The data can be downloaded here: https://ufile.io/wi6tt

library(sf)
library(tidyverse)
setwd("~/Desktop")

iowa<-st_read( dsn="Townships/iowa", layer="PLSS_Township_Boundaries", stringsAsFactors = F) # import data

## Make division 
r<-NULL
for (row in 1:nrow(iowa)) {
  r[[row]]<-st_make_grid(iowa[row,],n=c(3,3))
}

# Combine together
region<-NULL 
for (row in 1:nrow(iowa)) {
  region<-rbind(region,r[[row]])
}

region<-st_sfc(region,crs=4326) #convert to sfc
reg_id<-data.frame(reg_id=1:length(region)) #make ID for dataframe

# Make SF 
region_df<-st_sf(reg_id,region)

The last line gives the following error:

Error in `[[<-.data.frame`(`*tmp*`, all_sfc_names[i], value = list(list( : replacement has 1644 rows, data has 14796

1664 is the number of rows in the initial Iowa dataframe. Clearly the number of rows does not match the number of elements.

This might be a general r thing, rather than a spatial one, but I figured I'd post the whole thing in case someone had an idea on how to do the entirety of this a little cleaner