2

A source shapefile contains clear outlying points/area. How can they be isolated and removed from the spatial polygon in R?

library(rgeos) #required by maptools
library(maptools)

incheon <- readShapeSpatial("KOR_adm1.shp") #upload Korean provinces
incheon <- incheon[incheon@data$NAME_1=="Incheon",]
plot(incheon) #plot reveals clear outlier in the southeast corner of the map

enter image description here

5
  • 1
    It's usually better to investigate why you have an outlier and fix it in your dataset rather than just delete it from the analysis
    – Ian Turton
    Apr 24, 2020 at 9:38
  • @IanTurton I downloaded the shapefiles for South Korea from GADM, which looks like a reliable source (UC Davis). I've never ever heard that the city of Incheon has outlying possessions (apart from islands to its west). It seems highly unlikely. My idea is to fix the dataset in R. Are you saying it should be fixed elsewhere? Is it easier or more reliable to do in QGIS for example?
    – syre
    Apr 24, 2020 at 9:51
  • 1
    I would check with a domain expert (i.e. a South Korean) - also fixing it in QGIS is much easier (for me) than in R. Likely explanations are miss entered coordinates moving an area.
    – Ian Turton
    Apr 24, 2020 at 10:01
  • 1
    @IanTurton Misentered coordinates: you mean that this could be an actual piece of Incheon which got misplaced, in which case the solution is not deletion but correction of the coordinates. Correct?
    – syre
    Apr 24, 2020 at 10:04
  • 3
    It might be an error in the NAME_1 column. All the polygons that surrounds the "outliers" have the NAME_1="Busan". But as @IanTurton mentioned, it is better to check with a domain expert before you assume something.
    – Dataform
    Apr 24, 2020 at 10:52

2 Answers 2

1

First to the issue of your data. As hinted at by @Dataform in the comments, it appears that the Jung district of Busan has accidentally been coupled with the Jung district of Incheon in the GADM data. This is illustrated in the code below using the level 2 GADM data.

library(sf) # I'll use sf for this

sk <- st_read("gadm36_KOR_2.shp") # Read in South Korea
incheon <- sk[sk$NAME_1 == 'Incheon',] # Retrieve Incheon
busan <- sk[sk$NAME_1 == 'Busan',] # Retrieve Busan

incheon_jung <- incheon[incheon$NAME_2 == 'Jung',] # Jung district in Incheon
busan[busan$NAME_2 == 'Jung',] # Subsetting for Jung in Busan returns 0 rows

# Plotting only Busan
plot(st_geometry(busan))

enter image description here

Adding the Jung district from Incheon shows that the outlying area is supposed to be the Jung district of Busan

plot(st_geometry(incheon_jung), col = 'red', add = TRUE)

enter image description here

As per the comments from @Ian Turton I'd still be careful, but as the issue appears to be that the Jung district in Busan has accidentally been added to Incheon, and not erroneous coordinates, you should be able to just remove the outlying area and be left with Incheon. This can certainly also be done in QGIS, but I'll add a possible approach using R.

library(sf)

# Read in the level 1 South Korea GADM data that you are using
sk <- st_read("gadm36_KOR_1.shp")

# Use a projected coordinate system for clipping the data. 
# (A quick Google search suggested Korea 2000 / Unified CS (EPSG:5179)) 
sk <- st_transform(sk, 5179)

# Retrieve Incheon, including the outlying area
incheon <- sk[sk$NAME_1 == 'Incheon',]

# I manually inspected coordinates making up the bounding box for Incheon
(bb <- st_bbox(incheon))
xmin      ymin      xmax      ymax 
844186.2 1677674.0 1141187.7 1981910.8  

# And added 5000 to ymin to make sure the outlying area is outside the 
# extent of the coordinates
bb[2] <- bb[2] + 5000

# Make this a polygon
bpoly <- st_as_sfc(bb)

# Crop the Incheon data with this polygon
inchcor <- st_intersection(incheon, bpoly)

# And plot
plot(st_geometry(inchcor))

enter image description here

0

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

library(raster)
library(sp)
library(rgdal)
library(spatialEco)
library(rgeos)

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 explode.

incheon <- spatialEco::explode(incheon, sp = TRUE)

Here is where we create an extent SpatialPolgons object.

e <- as(raster::extent(125.6698, 126.8242, 36.76623, 37.90437), 
        "SpatialPolygons")

And, finally intersect the data to create the spatial subset.

incheon <- rgeos::gIntersection(incheon, e, byid=TRUE)
  plot(incheon)     

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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