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