if (FALSE) {
Lat.range = c(-10, 30)
Lon.range = c(65, 117)
######
########  load up the important libraries
 library(RFOC)
 library(geomapdata)
 data(coastmap)
  ###    data(ETOPO5)
####  need to download and install ETOPO data
load(ETOPO5)
 PLOC=list(LON=Lon.range,LAT=Lat.range,lon=Lon.range,lat=Lat.range,
                 x=Lon.range, y=Lat.range )
#####   set up topography colors
COLS = settopocol()
####  set the projection ##   utm
 PROJ = setPROJ(type=2, LAT0=mean(PLOC$y) , LON0=mean(PLOC$x) )   
 NK = 300
    
   ###   extract topography from the etopo5 data base in geomapdata
     JMAT = GEOTOPO(ETOPO5, PLOC, PROJ, COLS$calcol,nx=NK, ny=NK )
#####   select relevant earthquakes
IZ = list(x=JMAT$xo, y=JMAT$yo, z=JMAT$IZ$z)
CMAT = LandSeaCol(IZ, coastmap, PROJ, calcol=NULL)
Mollist =CMAT$Cmat
dMol = attr(Mollist, "Dcol")
     ####  Under water
UZ = CMAT$UZ
    #####  above water
AZ = CMAT$AZ
    ####   blues for underwater:
   blues = shade.col(100, acol=as.vector(col2rgb("darkblue")/255),
        bcol= as.vector(col2rgb("paleturquoise")/255))
plot(x=range(IZ$x), y=range(IZ$y),
         type='n', asp=1, axes=FALSE, ann=FALSE)
image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE)
image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE)
plotGEOmapXY(coastmap,
       LIM = c(Lon.range[1],Lat.range[1] ,Lon.range[2] ,Lat.range[2]),
       PROJ =PROJ,MAPstyle =2,MAPcol ="black" ,   add=TRUE  )
}
Run the code above in your browser using DataLab