3

I wish to transform a set of South Korean province shapes from CRS (Coordinate Reference System) GRS80, but st_transform returns empty geometries.

I have downloaded the dataset from this page (in Korean) which specifies the CRS as GRS80. I attach the file, as it is not accessible outside Korea.

Next is the R code resulting in empty geometries:

library(sf)

#Read shapes into R
sca_nsdi <- st_read("Z_NGII_N3A_G0010000.shp", stringsAsFactors=F)

#CRS is not specified
is.na(st_crs(sca_nsdi))

#Assign CRS as EPSG code of GRS80 per this reference: epsg.io/42310
st_crs(sca_nsdi) <- 42310 

#All geometries full
sum(!st_is_empty(sca_nsdi))==nrow(sca_nsdi)

#Transform CRS
sca_nsdi <- st_transform(sca_nsdi, 5179)

#Resulting geometries all empty
sum(st_is_empty(sca_nsdi))==nrow(sca_nsdi)
Improve this question

1 Answer 1

Reset to default
4

That data is not in that coordinate system.

If I read the data in, I see:

> sca_nsdi <- st_read("Z_NGII_N3A_G0010000.shp", stringsAsFactors=F)
Reading layer `Z_NGII_N3A_G0010000' from data source `/home/rowlings/Downloads/SO/korea/Z_NGII_N3A_G0010000.shp' using driver `ESRI Shapefile'
Simple feature collection with 38 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 22534.67 ymin: -42034.99 xmax: 632509 ymax: 545440.1
epsg (SRID):    NA
proj4string:    NA

and when I plot it I see a map of South Korea. But look at this:

bbox:           xmin: 22534.67 ymin: -42034.99 xmax: 632509 ymax: 545440.1

The ymin is showing a negative number, but for that coordinate system, according to epsg.io, that should be south of the equator. The equator does not run through South Korea. Hence your proposition that the data is in EPSG 42310 is wrong.

A little detective work with http://projfinder.com/ provides evidence that it is actually ESPG code 2097. Let's try that. Starting with the original shapefile data:

Edit and Note:

a little more detective work means its more likely epsg:5181, which applies to the whole of Korea (2097 says "central belt" and is mostly the west of the country) and also is based on the GRS 1980 ellipsoid which might be where you picked up the idea that any "GRS 80" CRS would do. That only refers to the modelled shape of the earth and not to the projection type, units, offset and scaling etc.

Doing the same as I've done below with 5181 doesn't make a visible difference to the maps and is probably only a few metres different to 2097. But I'd bet 5181 is more likely the projection used to create the data. Or maybe 5186, the year-2010 revision of 5181.

-- end edit note.

> st_crs(sca_nsdi)=2097   # probably should be 5181

Let's see where this is in lat-long, since we can verify that easily:

> st_transform(sca_nsdi,4326)
Simple feature collection with 38 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 125.0763 ymin: 33.11198 xmax: 131.8734 ymax: 38.40004
epsg (SRID):    4326

which looks like the correct location for South Korea. If we transform to your desired system, 5179, we see no empty polygons:

> sca_nsdi_transform = st_transform(sca_nsdi, 5179)
> any(st_is_empty(sca_nsdi_transform))
[1] FALSE

So your problems were two-fold. First, the Shapefile should have had a .prj file with the projection information. Ideally you should contact the source and find out why. Secondly your assumption (based on what?) that it was in 42310. A little inspection showed me this could not possibly be the case, and I used projfinder to match the location of Seoul to show it was possibly EPSG code 2097.

Overlaying your data (in 2097) with OpenStreetMap backgrounds shows very good alignment, so I'm confident this is good enough, even if it is not perfectly correct (it might be a metre or ten out). But if you can find out from authority what the CRS really is, use that.

enter image description here

Improve this answer
4
  • Fantastic work, thank you! Indeed, this post came after much unsuccessful wrangling to find out the CRS number, even contacting the publisher (but not the authoring organization) of that data. My assumption of GRS80 was because the download page of the data specifies this as the reference system, as indicated in the OP.
    – syre
    May 9, 2020 at 5:36
  • What was more specifically your "little detective work with projfinder.com"? Was it entering the coordinates of the centroid of that bbox into that page?
    – syre
    May 9, 2020 at 5:38
  • 1
    I plotted your map in R with its original coordinates and clicked on the location of Seoul to get Seoul's X-Y in the original coordinates. Then I zoomed projfinder to Seoul and put those coordinates in projfinder's X-Y box. Projfinder then goes over all its thousands of coordinate systems to see if any have lat-long on that map. Not sure why it didn't find 5181 or 5186.
    – Spacedman
    May 9, 2020 at 6:21
  • 1
    Google-translating the web page I see the coordinate system is "New GRS 1980 Central". EPSG 5186 is the newest that has "Central" and "Korea" in the description, and uses the GRS1980 ellipsoid. 99.999% sure that's the one!
    – Spacedman
    May 9, 2020 at 6:27

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.