You really are forced to define a spatial extent to bound the data. You can draw a bounding box on a plot screen using
raster::drawExtent(). One issue with this data is that it comes as MULTIPART geometry, which is one attribute row representing multiple spatial features. This is why a simple indexing approach or a geometry intersection operator will not work using an sp class object.
I would recommend exploding the polygon geometry into SINGLEPART.
Add libraries, data and subset
incheon <- rgeos::readOGR(getwd(), "gadm36_KOR_1")
incheon <-incheon[incheon$NAME_1 == "Incheon",]
Here you can explode the MULTIPART geometry. If you look at
dim(incheon) you will see 1 feature with 10 attributes and 111 features with 10 attributes after applying
incheon <- spatialEco::explode(incheon, sp = TRUE)
Here is where we create an extent
e <- as(raster::extent(125.6698, 126.8242, 36.76623, 37.90437),
And, finally intersect the data to create the spatial subset.
incheon <- rgeos::gIntersection(incheon, e, byid=TRUE)