17 June 2014

[#R] How to convert coordinate system (Lat-Long to UTM)

[#R] How to convert lat-long coordinates to UTM (easting-northing)

[#R] How to convert lat-long coordinates to UTM (easting-northing)

1. Intro (rgdal installation on Mac)

This bit is part of my work in modeling the hydrology of Cikapundung Catchment. I am kind of (forced) to do the spatial analysis in R :-). But after several trial and errors, it's really not that hard.

So If you do spatial analysis using R, you would need to set your coordinate into easting and northing system. Many R packages in spatial section only read this projection system (correct me if I'm wrong). One of the is geoR package. I'll talk about this later, since I'm still figuring out how to use it to make variograms and trend maps.

Now, to do a coordinate conversion, you would need rgdal package. This is the basic R package that deals with spatial analysis. Just install it with this command:

install.packages('rgdal')

In my case, the installation went well in Mac (running Snow Leopard) and also Windows 7, but it needed a workaround to install it on my Linux machine (Ubuntu 13.10). After installation, call the package using

library(rgdal) or require(rgdal)

2. Installing rgdal on Ubuntu

As I have said, installing rgdal on Linux was not as smooth as installing the other spatial packages, say maptools, geoR, or sp. At least in my case.

It has something to do with libgdal1-dev and libproj-dev files and their version which is corresponding to the Linux kernel version.

I googled it and this terminal command did work on my machine running Ubuntu 13.10. See this post:

sudo add-apt-repository ppa:marutter/c2d4u

sudo apt-get update

sudo apt-get install r-cran-rgdal

Another way (that did not work for me, but it might work on you). According to the same link above, this command was re-ran after switching off the related update site in the Ubuntu's update manager database:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable

sudo apt-get update

sudo apt-get install libproj-dev

Another way suggested by Robin Lovelace was working on Ubuntu 13.04.

3. Converting coordinates

Then load your data, in my case was

data = read.csv('27-1997data.csv') or

data <- read.csv('27-1997data.csv')

I often use the first one with “=” sign to avoid me hitting “shift key”.

Then you must set your existing coordinate system. Because otherwise, R will treat your x and y as plain data, not coordinate. In this case we are going to convert the existing lat-long coordinates to UTM (easting-northing system).

Run this command:

cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string=CRS("+proj=longlat"))

You can change the “cord.dec”, “data$long”, and “data$lat”, based on your dataset.

Then you are going to transform the CRS or Coordinate Reference System to UTM by setting the EPSG to 32748 for WGS 84, UTM Zone 48M, southern hemisphere. This is where Bandung is located.

cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32748"))

cord.UTM

4. Coordinate conversion on GIS app

I know that this post was mainly about R. But It does not hurt to share a system using QGIS and GRASS, posted by a fellow user Dmitri, from GRASS GIS Facebook Community.

So basically we will have two options:

Option 1

Use QGIS to reproject your data from WGS-84 (EPSG:4326) to UTM 48 zone (EPSG:34728). To do that in QGIS: (a) run QGIS, (b) load the dataset, © then use menu: Layer->Save As… (select new proection in the dialog). Then you can create new GRASS mapset using your reprojected data.

Option 2

Reproject the data using GRASS: (a) Create mapset in WGS-84 (b) Import the data, © create mapset in UTM 48 zone, (d) reproj data from the first mapset to the second (r.proj and v.proj commands).

5. The complete R codes.

Do not expect fancy code with me (LOL).

# loading library

require(rgdal)
## Loading required package: rgdal
## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj
data = read.csv("27-1997data.csv")
cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string = CRS("+proj=longlat"))


# Setting existing coordinate as lat-long system
cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string = CRS("+proj=longlat"))

# Transforming coordinate to UTM using EPSG=32748 for WGS=84, UTM Zone=48M,
# Southern Hemisphere)
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32748"))
cord.UTM
## SpatialPoints:
##       coords.x1 coords.x2
##  [1,]    783993  10767111
##  [2,]    778949  10758896
##  [3,]    784001  10765562
##  [4,]    789969  10766259
##  [5,]    783890  10765672
##  [6,]    785760  10767453
##  [7,]    777100  10753022
##  [8,]    776305  10757112
##  [9,]    776193  10757222
## [10,]    785220  10765016
## [11,]    779750  10753921
## [12,]    790006  10759509
## [13,]    785175  10752843
## [14,]    780604  10759679
## [15,]    784444  10765343
## [16,]    777882  10751587
## [17,]    782471  10762234
## [18,]    782032  10761457
## [19,]    783456  10764121
## [20,]    786240  10760484
## [21,]    775304  10758102
## [22,]    785175  10752843
## [23,]    785175  10752843
## [24,]    785175  10752843
## [25,]    785175  10752843
## [26,]    782141  10761790
## [27,]    791358  10754868
## [28,]    776638  10756892
## [29,]    779824  10760892
## [30,]    783354  10762571
## [31,]    783804  10761135
## [32,]    779808  10763769
## [33,]    785175  10752843
## [34,]    785175  10752843
## [35,]    785379  10756053
## [36,]    785175  10752843
## [37,]    778412  10755906
## [38,]    785310  10748307
## [39,]    784556  10765233
## [40,]    784116  10764678
## [41,]    785175  10752843
## [42,]    768459  10755854
## [43,]    785175  10752843
## [44,]    771336  10755537
## [45,]    785257  10758266
## [46,]    785175  10752843
## [47,]    781714  10759021
## [48,]    779287  10757791
## [49,]    785247  10760147
## [50,]    790201  10764158
## [51,]    784267  10757264
## Coordinate Reference System (CRS) arguments: +init=epsg:32748
## +proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0

# Plotting points
par(mfrow = c(1, 2))
plot(cord.dec, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95)
plot(cord.UTM, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95)

plot of chunk unnamed-chunk-1

No comments: