I come to these situations where I have to work with spatial datasets frequently. Sometimes the datasets is in lat-long format and sometimes on UTM (Universal transverse Mercator) coordinate or state plane coordinate system.

I extensively use EFDC (Environmental Fluid Dynamic Code) model to do environmental modeling and have to convert the data into UTM coordinate system. I had used CORPSCON in the past but as an R enthusiast I wanted to dig deeper in R and finad a way to do the transformation in R environment. Thankfully, I figured out that if we use the packages sp, rgdal, and raster, we can be able to convert the dataset (tiff or any x, y, z file) into lat-long / eastings-northings format.

Data download

To demonstrate the spatial transformation feature, we will be using the Lake Victoria dataset. Let’s go ahead and download lake victoria bathymetry data from Harvard using this link

I have saved tehe data as LV_Bathy_V6_1.tif and the size of the file is 48.9 MB. If you want to download the data with R uncomment the first couple lines in the code below.

Please install the packages if you don’t have these packages in your machine. You can use install.packages(c("raster", "sp", "rgdal")) to bulk download the packages.

library(raster)
## Loading required package: sp
library(sp)
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/DSI/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/DSI/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
#download.file(url = "https://dataverse.harvard.edu/api/access/datafile/3074951", 
#              destfile = "../../data/spatial_analysis/LV_Bathy_V6_1.tif")

filepath <- "../../data/spatial_analysis/LV_Bathy_V6_1.tif"
victoria <- raster(filepath)

Now that we have read the file in a raster victoria, let’s plot and see how the bathymetry looks like.

plot(victoria)

In the above figure, you can see that the deeper portion of the lake has dark green shade and the shallow region has grey color contours.

Convert raster to points

The raster dataset can be converted to SpatialPointsDataFrame using the rasterToPoints function in raster package.

victoria.pts <- raster::rasterToPoints(x = victoria, spatial = TRUE)
victoria.pts
## class       : SpatialPointsDataFrame 
## features    : 5788449 
## extent      : 683898.6, 1020899, -310550.6, 48449.37  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       :     LV_Bathy_V6_1 
## min values  : 0.839264631271362 
## max values  :  77.8775787353516

Change projection to geographic coordinate system

The projection of the new dataset is lcc i.e, Projected Coordinate System: Africa Lambert Conformal Conic. Now, let’s go ahead and convert it to geographic coordinate system.

geo.prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
victoria.pts <- spTransform(victoria.pts, CRS(geo.prj))
victoria.pts
## class       : SpatialPointsDataFrame 
## features    : 5788449 
## extent      : 31.60411, 34.85023, -3.002646, 0.4870742  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       :     LV_Bathy_V6_1 
## min values  : 0.839264631271362 
## max values  :  77.8775787353516
head(coordinates(victoria.pts))
##             x         y
## [1,] 33.36222 0.4869828
## [2,] 33.36319 0.4869866
## [3,] 33.36415 0.4869904
## [4,] 33.36512 0.4869942
## [5,] 33.36608 0.4869980
## [6,] 33.36704 0.4870018

Create x, y, z dataframe

Now let’s create x, y and z dataframe from the transformed data victoria.pts. The first column of victoria.pts is longitude and second column in latitude.

victoria.pts@data <- data.frame(victoria.pts@data, long = coordinates(victoria.pts)[,1], 
                                lat = coordinates(victoria.pts)[,2])
# View head of the dataset
head(victoria.pts@data)
##   LV_Bathy_V6_1     long       lat
## 1      5.511920 33.36222 0.4869828
## 2      5.517967 33.36319 0.4869866
## 3      5.515999 33.36415 0.4869904
## 4      5.514107 33.36512 0.4869942
## 5      5.512265 33.36608 0.4869980
## 6      5.510451 33.36704 0.4870018
dim(victoria.pts@data)
## [1] 5788449       3

Let’s look at the structure of the dataframe.

str(victoria.pts@data)
## 'data.frame':    5788449 obs. of  3 variables:
##  $ LV_Bathy_V6_1: num  5.51 5.52 5.52 5.51 5.51 ...
##  $ long         : num  33.4 33.4 33.4 33.4 33.4 ...
##  $ lat          : num  0.487 0.487 0.487 0.487 0.487 ...

If we want to convert lat long to UTM for lake victoria, we can use the following projection

utm.prj <- "+proj=utm +zone=16 ellps=WGS84"

One important thing to remember is that when you are specifying the options for the projection, don’t use any spacing. For example if you use the following the program / transformation would throw as error.

# Bad example
utm.prj <- "+proj = utm +zone = 16 ellps = WGS84"