get.GTOPO30<-function(filename= "W100N40",shrink=FALSE){

#
# this function will rad in a tile of the 
# global elevation data set GTOPO30
#

fileDEM<- paste(filename,".DEM", sep="")
fileHDR<- paste(filename,".HDR", sep="")

temp<- read.table(file=fileHDR, header=FALSE,stringsAsFactors=FALSE)
info<- as.double( temp[,2])
names( info)<- temp[,1]

nr<- info["NROWS"]
nc<- info["NCOLS"]

a<- info["ULXMAP"] + info["XDIM"]/2
b<- a+ (nc-1)* info["XDIM"]
xtemp<- seq(a, b,,nc)

b<- info["ULYMAP"] - info["YDIM"]/2
a<- b - (nr-1)* info["YDIM"]
ytemp<- seq(a, b,,nr)



handle<-file(fileDEM, "rb")
open( handle)

ztemp<-  readBin(handle, what=as.integer(1) ,nc*nr, size=2, endian="big")
close( handle)

#permute into R image format  (1,1) element is bottom left 
# and set to R missing value code
ztemp<- matrix( ztemp, nrow=nc)[,nr:1]
ztemp[ ztemp== info["NODATA"]] <- NA

if( !shrink ){
# usually you want this one -- full resolution.
  list(x= xtemp, y=ytemp, z= ztemp)}
else{
   half.image( half.image( list(x= xtemp, y=ytemp, z= ztemp)))
  cat("Note: Raw image averaged by factor of 4  (4X4)-> 1 \n
           using half.image function aplied twice ", fill=TRUE)
}


##ind<- 3000 + (1::4000
#image( out$x[ind], out$y[ind], out$z[ind,ind])
#world( add=TRUE, col="grey", lwd=2)
}