Paragon Corpoation PostGIS Spatial Database Engine OSGeo.org The Open Source Geospatial Foundation UMN Mapserver Boston Geographic Information Systems       PostGreSQL Object Relational Database Management System
Home   About Boston GIS   Consulting Services  Boston GIS Blog  Postgres OnLine Journal  Planet PostGIS  PostGIS Funding

Purpose of BostonGIS

BostonGIS is a testbed for GIS and Web Mapping solutions utilizing open source, freely available and/or open gis technologies. We will be using mostly Boston, Massachusetts data to provide mapping and spatial database examples.

If you have some thoughts or comments on what you would like to see covered on this site, drop us a line on our Feed Back page.


GIS Tutorials on Opensource and OpenGIS technologies Tutorials
GIS Article comments Article and Tutorial Comments
Boston GIS BLog Rss FeedBoston GIS blog

PDF HTML All BostonGIS tutorials packaged together in an E-Book.


Tutorial and Tip Sites
Desktop GIS
External Data
GIS Events and Groups
GIS SDKs and Frameworks
External Resources
Glossary
GIS Blogs Around Boston
External GIS Blogs
External Papers Articles
GIS Quick Guides and References
OpenStreetMap and OpenLayers Tutorials
PostGIS, pgRouting, and PostgreSQL Tutorials
Part 1: Getting Started With PostGIS: An almost Idiot's Guide more ...
pgRouting: Loading OpenStreetMap with Osm2Po and route querying more ...
Part 1: Getting Started With PostGIS: An almost Idiot's Guide (PostGIS 2.0) more ...
OSCON 2009: Tips and Tricks for Writing PostGIS Spatial Queries more ...
PGCon2009: PostGIS 1.4, PostgreSQL 8.4 Spatial Analysis Queries, Building Geometries, Open Jump more ...
PLR Part 3: PL/R and Geospatial Data Abstraction Library (GDAL) RGDAL

Printer Friendly

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

GDAL.close(TestLogo)

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

displayDataset(custarea,offset=c(0,0), region.dim=dim(custarea), reduction=1, band=1)

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

displayDataset(custarea,offset=c(0,0), region.dim=dim(custarea), reduction=1, band=1:3)

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

GDAL.close(custarea)

Reading Data From Raster file

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

Load into a spatial dataframe table
custdf <-open.SpatialGDAL("/path/to/CustomAreaSRTM90.dem")
You should see output such as this:

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



How to read the projection information

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)



How to read bounding box

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

How to read the projection information and other meta data
GDALinfo("/path/to/CustomAreaSRTM90.dem")
cust_summary <- GDALinfo("/path/to/CustomAreaSRTM90.dem")
cust_summary["rows"]
cust_summary["ll.x"]
cust_summary["ll.y"]


Return a RasterTable containing lat long and band data

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
:
:
custarea_gt["band1"]

Will give you just the band1 column which gives the elevation levels.

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..

SELECT getTileInfo('/path/to/CustomAreaSRTM90.dem', 'll.x') As x_pos, getTileInfo('/path/to/CustomAreaSRTM90.dem', 'll.y') As y_pos, getTileInfo('/path/to/CustomAreaSRTM90.dem', 'rows') As num_rows, getTileInfo('/path/to/CustomAreaSRTM90.dem', 'columns') As num_cols



Post Comments About PLR Part 3: PL/R and Geospatial Data Abstraction Library (GDAL) RGDAL
PostGIS Nearest Neighbor: A Generic Solution - Much Faster than Previous Solution more ...
Solving the Nearest Neighbor Problem in PostGIS more ...
PLR Part 2: PL/R and PostGIS more ...
PLR Part 1: Up and Running with PL/R (PLR) in PostgreSQL: An almost Idiot's Guide more ...
Part 3: PostGIS Loading Data from Non-Spatial Sources more ...
Part 2: Introduction to Spatial Queries and SFSQL with PostGIS more ...
Miscellaneous Tutorials/Cheatsheets/Examples
SpatiaLite Tutorials
Boston External Map Examples
SQL Server 2008 Tutorials
UMN Mapserver Tutorials
General Commentary
Boston GIS      Copyright 2024      Paragon Corporation