Blue marble and Landsat in R 3.

botanical-collections.png


A final note on this theme. To add points to a 3d visualisation you need to know z values in order to place them on the surface. The way x,y,z are referred to by rgl is not intuitive and it took me some time to get this to work. The example uses a set of coordinates for the botanical collections and inventories of tree species that are compiled into a database. We are using these data for modelling regional tree species distributions and their potential response to climate change. A major source of this information has been the wonderful online resource provided by Missouri Botanical Gardens MOBOT, but data has also been provided from other sources including Conabio. Here I only provide the raw coordinates of collections, without species names or further details as some of the data is restricted. The interesting point to note is the spatial distribution. Notice that these points are for all trees, thus any one species will form an extremely small subset. Our knowledge of the real distribution of tree species in Mexico is thus clearly still very dependent on expert knowledge and uncompiled data from little explored areas.

The code is here. You should load rgl and sp before running it.

library(sp)
library(rgl)
load(url("https://duncanjg.files.wordpress.com/2008/03/colpoints.doc"))
load(url("https://duncanjg.files.wordpress.com/2008/03/bmmexicoapril.doc"))

coordinates(a)<-~x+y
a<-overlay(d,a)
a<-as.data.frame(a)

view.3d.rgb<-function(d,dem=4,red=1,blue=2,green=3,exag=1/30000){
 fullgrid(d)<-TRUE
  y<-d[[dem]]*exag
  y[is.na(y)]<-0
  dim(y)<-d@grid@cells.dim

x<-seq(min(d@coords[,1]),max(d@coords[,1]),length=d@grid@cells.dim[1]) 
 z<-seq(min(d@coords[,2]),max(d@coords[,2]),length=d@grid@cells.dim[2])

 colvec<-  rgb(d[[red]],d[[green]],d[[blue]],m=255)

  rgl.clear()
  rgl.surface(x, z, y, color=colvec)
}



add.points.3d<-function(d=d,x,y,elev,exag=1/30000,...){
rgl.points(x=x,z=max(d@coords[,2])-y+min(d@coords[,2]), y=elev*exag,...)}

view.3d.rgb(d)
add.points.3d(d,a$x,a$y,a$alt,col=2,size=2)


addpoints3d.doc

2 thoughts on “Blue marble and Landsat in R 3.

  1. When running this script I get an error “Error in fullgrid(d) <- TRUE : could not find function "fullgrid<-""

    Not sure why?

    • You need to load the packages sp and rgl
      library(sp)
      library(rgl)

      Then it runs. But the page needs updating to place the code on the page itself.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s