I am looking for a reliable shapefile of the Seoul Capital Area (SCA), comprising Seoul, Incheon and Gyeonggi-do provinces. One common source is GADM. However, the coastline is not up-to-date so I have found another source in Korea's NSDI (in Korean) which includes a more up-to-date coastline. (I attach the files, as they are not accessible outside Korea.)

However, I find that the datasets are substantially non-aligned, so I wonder if I have made a mistake or whether the problem lies in one of the sources. As can be seen on the map below, the GADM regions (black outline) display a significant shift to the south of the NSDI regions (red outline). The shift appears to increase in the eastward direction, westernmost Incheon being well aligned, which makes me suspect that this is a projection issue. Perhaps, the wrong CRS is being applied somewhere. enter image description here

Here's my code:


#***1 GADM data processing***

#Load GADM shapefile
sca_gadm <- readRDS("gadm36_KOR_1_sf.rds") #class "sf" "data.frame"

#Subset SCA
sca_gadm <- sca_gadm[sca_gadm$NAME_1 %in% c("Seoul", "Incheon", "Gyeonggi-do"),]

#Flatten coordinate system
sca_gadm <- st_transform(sca_gadm, 5179)

#Remove Incheon outlying polygon error (courtesy of https://gis.stackexchange.com/a/359722/88704)

#Inspect coordinates of bounding box
(bb <- st_bbox(sca_gadm)) #class "bbox"
#xmin      ymin      xmax      ymax 
#844186.2 1677674.0 1141187.7 2024763.7

#Add 50000 to ymin to make sure outlying area is outside coordinates
bb[2] <- bb[2] + 50000

#Crop SCA with bb polygon
sca_gadm <- st_intersection(sca_gadm, st_as_sfc(bb))

#***2 NSDI data processing***

#Dissolve NSDI sub-regions for comparison with GADM (courtesy of https://gis.stackexchange.com/a/361510/88704)

#Load, buffer, dissolve, reverse buffer: Seoul
sel_nsdi <- cbind(st_read("LARD_ADM_SECT_SGG_11.shp"), Name="Seoul")
sel_nsdi <- st_buffer(sel_nsdi, 35)
sel_nsdi$geometry <- st_union(sel_nsdi$geometry)
sel_nsdi <- st_buffer(sel_nsdi, -35)

#Load, buffer, dissolve, reverse buffer: Incheon
inc_nsdi <- cbind(st_read("LARD_ADM_SECT_SGG_28.shp"), Name="Incheon")
inc_nsdi <- st_buffer(inc_nsdi, 25)
inc_nsdi$geometry <- st_union(inc_nsdi$geometry)
inc_nsdi <- st_buffer(inc_nsdi, -25)

#Load, buffer, dissolve, reverse buffer: Gyeonggi-do
ggd_nsdi <- cbind(st_read("LARD_ADM_SECT_SGG_41.shp"), Name="Gyeonggi-do")
ggd_nsdi <- st_buffer(ggd_nsdi, 50)
ggd_nsdi$geometry <- st_union(ggd_nsdi$geometry)
ggd_nsdi <- st_buffer(ggd_nsdi, -50)

#Combine sub-regions into SCA
sca_nsdi <- rbind(sel_nsdi, inc_nsdi, ggd_nsdi)

#***3 Overlay GADM and NSDI regions***

#Both datasets share same EPSG:5179
st_crs(sca_nsdi, parameters=T)$epsg==st_crs(sca_gadm, parameters=T)$epsg

#Plot Overlaying regions
plot(sca_nsdi$geometry, border="red")
plot(sca_gadm$geometry, add=T)


Your Answer

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

Browse other questions tagged or ask your own question.