What is GDAL and RGDAL?

GDAL stands for Geospatial Data Abstraction Library and is a popular open source GIS library originally developed and maintained by Frank Warmerdam with contributions from many. It is used in GIS for processing various types of GIS datasources - both Vector and Raster data. The GDAL part generally refers to RASTER functionality while the complementary OGR refers to Vector data. The RASTER functionality is very useful in analysis such as analysing Weather satellite raster data and Digital Elevation Models (DEMS) or converting between RASTER formats.

RGDAL is an R-statistical package that makes these library functions available from the R statistical environment (both RASTER and Vector functions) - thus making them available within PostgreSQL PL-R environment. In this exercise we shall install and play with this library first from R and then from within PostgreSQL.

Installing RGDAL

For this exercise, we assume you have already installed PostgreSQL and PL/R. If you haven’t already, read our PLR Part 1: Up and Running with PL/R (PLR) in PostgreSQL: An almost Idiot's Guide - http://www.bostongis.com/PrinterFriendly.aspx?content_name=postgresql_plr_tut01

RGDAL relies on the package sp which contains the SpatialGridDataFrame class among other things, so you'll need to install that first. To install sp do the following:

1. Launch your R GUI environment
2. chooseCRANmirror() - not all CRAN Mirrors have the sp package so try the one closest to you and if that one doesn't have it, then try another.
3. utils:::menuInstallPkgs()
4. SELECT sp from list of options

1. Download the RGDAL library and manual from here - http://cran.r-project.org/web/packages/rgdal/index.html.
   
   NOTE: For Windows users, you can use the binary provided if you do not want to compile yourself.
2. Once compiled or binary unzipped - copy the resutling rgdal folder into library folder of your R install.

Testing out RGDAL in R Environment

Launch RGUI. A lot of these excercises are similar to ones documented in the RGDAL help file. Now lets test drive the library. Locate a small image on one of your drives. E.g. for windows users, you can copy an image to root of C of server or something of that sort for linux usual /usr/local/whatever will work fine. Then do something like.

library(rgdal)
TestLogo <- readGDAL("C:/test.jpg")
TestLogo

The above should read the jpeg into TestLogo variable and typing TestLogo will spit out all sorts of details about the Jpeg file. That will look something like below:

test.jpg has GDAL driver JPEG 
and has 9 rows and 9 columns
Closing GDAL dataset handle 0x048f1230...  destroyed ... done.
> TestLogo
Object of class SpatialGridDataFrame
Object of class SpatialGrid
Grid topology:
  cellcentre.offset cellsize cells.dim
x               0.5        1         9
y               0.5        1         9
SpatialPoints:
        x   y
 [1,] 0.5 8.5
 [2,] 1.5 8.5
 [3,] 2.5 8.5
 [4,] 3.5 8.5
 [5,] 4.5 8.5
 [6,] 5.5 8.5
 [7,] 6.5 8.5
 [8,] 7.5 8.5
:
:
Coordinate Reference System (CRS) arguments: NA 

Data summary:
     band1            band2           band3      
 Min.   : 40.00   Min.   : 59.0   Min.   :117.0  
 1st Qu.: 51.00   1st Qu.: 71.0   1st Qu.:124.0  
 Median : 52.00   Median : 72.0   Median :125.0  
 Mean   : 86.99   Mean   :103.7   Mean   :147.6  
 3rd Qu.: 61.00   3rd Qu.: 79.0   3rd Qu.:132.0  
 Max.   :255.00   Max.   :255.0   Max.   :255.0 

For cleanup and destroy TestLogo variable

Now lets try with some real Geospatial data.

Download Digital Elevation Model (DEM) International sample from http://www.mapmart.com/DEM/DEM.htm and extract files into a folder.

Displaying data from RASTER file

The below R code will open up the CustomArea file and start at position Y=8, x=6 of the file and show 5 rows and 3 columns on the R graphics device.

custarea <- GDAL.open("/path/to/CustomAreaSRTM90.dem")
displayDataset(custarea,offset=c(6,8), region.dim=c(5,3), reduction=1, band=1)

To show the whole file you would do

If you had multiple color bands and you wanted to display all bands, then you would do

When we are done with the file, we should close it to clean up loose-ends

Reading Data From Raster file

Here are a couple of useful exercises that you would commonly do on a DEM or GeoTiff

CustomAreaSRTM90.dem has GDAL driver USGSDEM 
and has 3525 rows and 6275 columns

This is slower as it reads the whole file into memory
custarea <- readGDAL("/path/to/CustomAreaSRTM90.dem")
proj4string(custarea)

This is faster as it just opens the file and , but file needs to be closed

custarea <- open.SpatialGDAL("/path/to/CustomAreaSRTM90.dem", read.only=TRUE)
proj4string(custarea)
close(custarea)

custarea <- open.SpatialGDAL("/path/to/CustomAreaSRTM90.dem", read.only=TRUE)
bbox(custarea)
close(custarea)

GDALinfo("/path/to/CustomAreaSRTM90.dem")
cust_summary <- GDALinfo("/path/to/CustomAreaSRTM90.dem")
cust_summary["rows"]
cust_summary["ll.x"]
cust_summary["ll.y"]

custarea <- GDAL.open("/path/to/CustomAreaSRTM90.dem")
custarea_gt <- getRasterTable(custarea, band = 1, offset = c(0,0), region.dim = c(10,10))
GDAL.close(custarea)
custarea_gt

The result of the above looks something like this

            x        y band1
1   -127.4554 50.23167   438
2   -127.4546 50.23167   410
3   -127.4537 50.23167   382
4   -127.4529 50.23167   356
5   -127.4521 50.23167   329
6   -127.4512 50.23167   306
7   -127.4504 50.23167   288
8   -127.4496 50.23167   272
9   -127.4487 50.23167   262
10  -127.4479 50.23167   257
11  -127.4554 50.23084   399
12  -127.4546 50.23084   378
13  -127.4537 50.23084   355
14  -127.4529 50.23084   333
15  -127.4521 50.23084   310
16  -127.4512 50.23084   290
:
:

Using RGDAL from PostgreSQL

Now switch over to PostgreSQL and use PgAdmin or psql to write the function. Below is a simple PL/R function that loads the R library and reads the meta data from a file passed in.

CREATE or REPLACE FUNCTION getTileInfo(atilefile varchar, element varchar) returns text AS
$$
    library(rgdal)
    tile_summary <-GDALinfo(atilefile)
    return(tile_summary[element])

$$
language 'plr';

Here is how we can use the above function from within PostgreSQL that returns a bit of info from the file we pass in as argument and the element piece we want. Keep in mind the path is relative to the PostgreSQL server since R is installed on the server and that this function will work for any GDAL supported file format such as DEM or GeoTiff or JPEG etc..