Duncan Golicher’s weblog

R is of course a statistical environment, so it really shouldn’t be expected to be doing the job of GRASS (Nviz), Google Earth, WorldWind, Qgis or Udig as well. However some nice visualizations can be done with rgl, including rgb composites from satellite images. They can even be more useful for quick regional overviews than the alternatives mentioned above.

The code to produce these images be found here bluemarble3d.doc. The first line will download the imported spatial object from this site. R does need a lot of memory for images of this size, so I guess you need a minimum of 1GB RAM (I am using a Toshiba Tecra laptop with Nvidia running Ubuntu Feisty 7.04 with 2GB). Given that, you should get very own zoomable 3d image of Mexico and Central America to play with in the same time it takes the screenshot below to run. It shows the steps in real time. (For R beginners, just open the file, paste all the code into the R console and wait a little for the download). The resolution is quite coarse (2 minute) so this is only suitable for a regional overview. You will need the rgl and sp packages installed first of course.

bluemarble3d.doc

Exactly he same can be done with Landsat imagery. landsat.doc

If you want to get your own coverages for this for your own region the starting point is to run R from GRASS in a lat-long location. For example

g.region -p

projection: 3 (Latitude-Longitude)
datum: wgs84
ellipsoid: a=6378137 es=0.006694379990141317
north: 32:51:47.86586N
south: 5:28:57.975592N
west: 118:30:20.609692W
east: 76:26:43.294846W
nsres: 0:02:00.060768
ewres: 0:01:59.982024
rows: 821
cols: 1262
cells: 1036102

Then get some of the Blue Marble (Modis) imagery that can be downloaded from the NASA wms strait into GRASS using r.in.onearth. From within R you can use system to run GRASS. This does work with GRASS under Cygwin in Windows with the shell command.

system(”r.in.onearth -b file=/tmp/test month=Apr time=2005-3-24 ‘wgetopt=-c -t 5 –user-agent=MSIE5.5′ “)

The layers can then be imported into R with spgrass6.

library(spgrass6)

d<-readRAST6(c(”BMNG_Apr.red”,”BMNG_Apr.blue”,”BMNG_Apr.green”,”alt”))

Notice that this line imported a digital elevation model of my own as well, as one is of course needed for the 3d terrain effect.

The animation above first shows an accessibility index calculated as a cost surface with the entrance to the watershed as a starting point. There is an excellent “How To” on calculating cost surfaces in GRASS available here. The cost surface (warm colours accessible, cold colours inaccessible) is followed by an animation of deforestation estimated by supervised classification of a series of landsat images taken in1975, 1987, 1999 and 2003. The area is a focus case study, the Tablon watershed in the Sepaltura biosphere reserve. The set of images in Google Earth format are available here.

Landsat imagery is the most frequently used for tracing patterns of deforestation over time. Landsat has been acquiring coverage of the Earth’s surface since 1972 when Landsat 1 was launched. Since then, four other satellites have been in operation. Landsat 1, 2, and 3 flew in a circular orbit 913 kilometers (570 miles) above the Earth’s surface and circled the Earth every 103 minutes (14 times a day). Landsat 4 and 5 fly about 705 kilometers (440 miles) above the Earth and circle every 98 minutes.

Landsat 4 and 5 are still operating. Landsat 6 was launched in 1993. Thus the sensors used to provide images have changed over time and are not easily comparable. The current Thematic Mapper (TM) class of Landsat sensors began with Landsat-4 (1982). This series continued with the nearly identical sensor on Landsat-5 (1984). Landsat-7 Enhanced Thematic Mapper Plus (ETM+) was carried into orbit in 1999. Currently, both the Landsat-5 TM and the Landsat-7 ETM+ are operational and providing data. Landsat-7 ETM+ experienced a failure of its Scan Line Corrector mechanism in May 2003. Data products have been developed to fill these gaps using other ETM+ scenes.

The images used for this analysis were downloaded from the global land cover facility and processed in GRASS. We are currently engaged in a classifying deforestation over a much larger area using similar imagery. Quantifying patterns of deforestation from these historical satellite images is never straightforward, for many reasons. Differences in the conditions when the images where taken can lead to misclassification of pixels leading to errors when documenting the overall spatial pattern. However even though it can be difficult to tell with certainty if a specific small area has changed, the general regional trends detected are quite reliable.

Recently high resolution images of this particular study region have become available and we are acquiring more. It is encouraging to find that at least the classification into “forest-non-forest” that has been achieved using Landsat imagery provides a very acceptable match with the visual impression provided by high resolution images. This can be seen below. You can check your own visual impression of the match by downloading the KML files and changing the transparency in Google Earth. The following lines load the layers into R.

urlstring<-”ftp://anonymous@200.23.34.16/simulacion/Datos”
a<-paste(urlstring,”tablondefor.rob”,sep=”/”)
load(url(a))