Search code examples
rspatialrgdalr-spgstat

Problems is projecting LON/LAT data in R

I have a dataframe contains PM10 concentration in air over seoul(capital city) in Korea. Please, take a look. I want to plot semivariogram from this data set. As LAT/LON data here, is in degree so I have project this data. I have projected data in this way:

library(rgdal)

seoul3112 <- read.csv("seoul3112.csv", row.name=1)
seoul3112 <- na.omit(seoul3112)

coordinates(seoul3112) <- ~LON+LAT
proj4string(seoul3112) <- "+proj=longlat +datum=WGS84"

seoul3112

after projecting I got seoul311 as like below

        coordinates PM10
1    (126.976, 37.56464)   42
2    (127.005, 37.57203)   37
3   (127.0051, 37.54031)   46
4   (127.0957, 37.54464)   47
5   (127.0411, 37.54311)   46

Q1: I found that after projecting, the value of LON/LAT show the nearly same value as previous data frame. My question is whats the actual function of this proj4string(seoul311) = "+proj=longlat +datum=WGS84" command. Here, LON/LAT(degree) is transferred to km/m or something like that?

I tried to write another code using rgdal package like below:

proj4string(seoul3112) <- "+proj=longlat +datum=WGS84" 
seoul3112 <- spTransform(seoul3112, 
                         CRS("+proj=utm +north +zone=52 +datum=WGS84"))
seoul3112

after projecting I got seoul3112 as like below

       coordinates     ID       time PM10
12      (321241, 4159438) 111121 2012030112   68
173   (323824.6, 4160203) 111123 2012030112   64
334   (323754.6, 4156684) 111131 2012030112   67
495   (331771.9, 4156998) 111141 2012030112   65
656   (326946.2, 4156927) 111142 2012030112   69

Q2. Here I can see the LON/LAT value transformed to some large value! whats the meaning of these value? m/km or something like that? in above code north means what? northern hemisphere?

Q3. As I mentioned earlier, I want to plot semivarigram over seoul in Korea (utm zone 52). So, which projection rule I should use? Should I consider utm zone? when should I consider utm zone?

I have many confusion about projecting the data. could you please answer my three question in details?

Solution

  • (FYI: It's usually bad form for 3 questions in a single SO post)

    Q1: You didn't actually "project" anything in the first block of operations. You created a "spatial" object from a plain data frame and "stated" what coordinate reference system (CRS) it was in. You also did that accurately as you just had lat/lon values. Do a str(seoul3112) to see the structure of the SpatialPointsDataFrame you ended up creating.

    Q2: You actually did "project" the coordinates to the Universal Transverse Mercator (UTM) CRS. UTM grid coordinates are expressed as a distance in meters to the east, referred to as the "easting", and a distance in meters to the north, referred to as the "northing".

    Q3: You should check with the suggested "official" government projection recommendations, but you can probably get away with something like azimuth equidistant for South Korea (and it's supported in mapproject so it's easy to work with in ggplot):

    library(ggplot2)
library(ggthemes)
library(mapdata)

seoul3112 <- read.csv("seoul3112.csv", row.name=1)
seoul3112 <- na.omit(seoul3112)

sk <- map_data("worldHires", "South Korea")

gg <- ggplot()
gg <- gg + geom_map(data=sk, map=sk,
                    aes(x=long, y=lat, map_id=region),
                    color="black", fill="white", size=0.25)
gg <- gg + geom_point(data=seoul3112, aes(x=LON, y=LAT))
gg <- gg + coord_map("azequidistant")
gg <- gg + theme_map()
gg

    enter image description here

    As you can see, depending on the projection (mapproject only supports 41) ggplot can take alleviate the need to project the points first.

    BUT, you are computing the semivariogram, so if you want to work in meters, you can do:

    coordinates(seoul3112) <- ~LON+LAT
proj4string(seoul3112) <- "+proj=longlat +datum=WGS84"
seoul3112_utm <- spTransform(seoul3112, 
                             CRS("+proj=utm +north +zone=52 +datum=WGS84"))

proj_3112 <- as.data.frame(coordinates(seoul3112_utm))
proj_3112 <- cbind.data.frame(proj_3112, seoul3112_utm@data)

    and then compute the distances (I'm assuming much here):

    dists <- dist(proj_3112[,1:2])

    and then finish your model. nlme, geoR and a few other R packages can help with the semivariogram model development and calculation (and plotting).