# ETOPO2_NEpacific.dat was downloaded from http://www.ngdc.noaa.gov/mgg/gdas/gd_designagrid.html
# using the ETOPO2 2-minute Global Relief and saved in XYZ format with comma separation, and renamed.
# created by Ian Taylor June 21, 2006
# share and improve at will

bath = read.csv('C:/dogfish/maps/ETOPO2_NEpacific.dat', header=F) #obviously you'll need to change the path
names(bath) = c('lat','lon','dep')

bath = round(bath,4)                  #rounding off to the same digit so that they will match exactly
latvec = round(seq(-135,-120,1/30),4) #these are the values that I used to download the data
lonvec = round(seq(30,55,1/30),4)

# putting the data into matrices
mat1 = matrix(bath$lat,length(lonvec),length(latvec),byrow=T)
mat2 = matrix(bath$lon,length(lonvec),length(latvec),byrow=T)
mat3 = matrix(bath$dep,length(lonvec),length(latvec),byrow=T)

# using mat1 and mat2 to check that the matrices came out right
summary(mat1[1,]==latvec)
#   Mode    TRUE
#logical     451
summary(mat2[,1]==lonvec)
#   Mode   FALSE    TRUE
#logical     750       1
summary(mat2[nrow(mat2):1,1]==lonvec) #shows that lonvec is reversed campared to matrix values
#   Mode    TRUE
#logical     751

depmat = t(mat3[nrow(mat3):1,]) #flipping upside down and transposing to match lonvec in the contour command

#
dep100 = contourLines(latvec,lonvec,depmat, levels=-100)
dep200 = contourLines(latvec,lonvec,depmat, levels=-200)
dep500 = contourLines(latvec,lonvec,depmat, levels=-500)
dep1000 = contourLines(latvec,lonvec,depmat, levels=-1000)

listLen100=NULL
for(i in 1:length(dep100)) listLen100 = c(listLen100,length(dep100[[i]]$x))
listLen100
#  [1]  138 1988    5    7   17    5    9    5    5    5   43    5   13    7    5    7    5    7   19    5    7
# [22]   13    5    5    5    5    9    5    5    5    5    5   11    5    9    5   11    4    5   11    5   99
# [43]    5    5    7    5   11   13   11    9    5   13    5    5    5   37    7    7   29    5    5   11   13
# [64]    9   23   11    5    5    5    5    5   13    5    5    5    7    5    5  225    5    5    5    5   11
# [85]    7    5    5    5    5    5    5   11   29    7    5    9    5   13   11    5    7   27    5    5   13
#[106]    5   13    5    7    5    5    5   11    5    9    7    7    7    9    7    7    7    5    5    5    5
#[127]   19    5    5    5   21    5    5   49    3
contour100 = cbind(dep100[[2]]$x, dep100[[2]]$y) # 2nd element in the this is the biggest one

listLen200=NULL
for(i in 1:length(dep200)) listLen200 = c(listLen200,length(dep200[[i]]$x))
listLen200
# [1] 1858    5   11    5    9   23    7    5    5    7    5    5    9   27   10   11    5   13    5    5   83    5
#[23]    5    4   25   11    7    5    5    5    9    5    5    5    9   25    7   11    5   13    5    7    5   11
#[45]    5   59    9   13   11    5   13   19    5    5    5    5    5  115    5    7    5    5   21   15    5    5
#[67]   21    7    5   11    5   64    5
contour200 = cbind(dep200[[1]]$x, dep200[[1]]$y) # 1st element in the this is the biggest one

listLen500=NULL
for(i in 1:length(dep500)) listLen500 = c(listLen500,length(dep500[[i]]$x))
listLen500
# [1] 1425    5    9    7   15    5    5    9    7   33   23   11    9    7    7    9    5   11    9    9    7    9
#[23]    5   43    7    7    7   13    9    9   17    9
contour500 = cbind(dep500[[1]]$x, dep500[[1]]$y) # 1st element in the this is the biggest one

listLen1000=NULL
for(i in 1:length(dep1000)) listLen1000 = c(listLen1000,length(dep1000[[i]]$x))
listLen1000
# [1]    5 1563    9    5    9   17    7   13   21    5   15    5   17    5    7    5   19    7   15    5    5    7
#[23]    5   13    7    7   11    5    7    7    9   43    7
contour1000 = cbind(dep1000[[2]]$x, dep1000[[2]]$y) # 2nd element in the this is the biggest one

library(maps)
library(mapdata)

windows(5,5)
map('worldHires', xlim=c(-128,-121),ylim=c(46,50), border=0, fill=TRUE,col='navajowhite3')
box();axis(1, seq(-128,-120,2,),paste(seq(128,120,-2),'°W', sep=''));axis(2, seq(45,50,1),paste(seq(45,50,1),'°N', sep=''))
# add thicker line for just the main coastal contour segment
lines(contour100,col=4, lwd=2)
lines(contour200,col=2, lwd=2)
lines(contour500,col=3, lwd=2)
lines(contour1000,col=5, lwd=2)
box()
legend('bottomleft',c('100 m','200 m','500 m', '1000 m'), lwd=2, col=c(4,2,3,5))
# add a thin line for all contour segments (just to show what could be thrown away)
contour(latvec,lonvec,depmat, levels=-100, add=T, col=4, drawlabels=F)
contour(latvec,lonvec,depmat, levels=-200, add=T, col=2, drawlabels=F)
contour(latvec,lonvec,depmat, levels=-500, add=T, col=3, drawlabels=F)
contour(latvec,lonvec,depmat, levels=-1000, add=T, col=5, drawlabels=F)

# whole US coast
windows(5,7)
map('worldHires', xlim=c(-128,-122),ylim=c(40,50), border=0, fill=TRUE,col='navajowhite3')
box();axis(1, seq(-128,-120,2,),paste(seq(128,120,-2),'°W', sep=''));axis(2, seq(40,50,2),paste(seq(40,50,2),'°N', sep=''))
lines(contour100,col=4)
lines(contour200,col=2)
lines(contour500,col=3)
lines(contour1000,col=5)
box()