#dev.new(width=10,height=7) # only needed once
pdf("fs1_regions_map.pdf", width = 10, height = 7)
library(marmap)

infile <- "data_dx1000_dz5_1.1.dat"
source("read_mapped_data.R")

world <- getNOAA.bathy(lon1 = -180, lon2 = 180, lat1 = -90, lat2 = 90, resolution = 60, keep = T)
xlm <- c(-173,173)
ylm <- c(-75,85)
plot(c(0),col="white",xlim=xlm,ylim=ylm,xlab="Longitude",ylab="Latitude",yaxt="n")
axis(2,at=c(-60,-30,0,30,60))
title(main="Regions",cex.main=1.7)
plot(world, deep = 0, shallow = 0,step=0,image=T, land=T,col="grey",add = TRUE,xlim=xlm,ylim=ylm,bpal = list(c(0, max(world), "lightgrey")))
par(new=T)
plot(lonb[iEEqPacb],latb[iEEqPacb],col="red",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-120,10,"ETPac")
par(new=T)
plot(lonb[iWEqPacb],latb[iWEqPacb],col="orange",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(150,10,"WTPac")
par(new=T)
plot(lonb[iNEPacb],latb[iNEPacb],col="blue",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-150,40,"NEPac")
par(new=T)
plot(lonb[iNWPacb],latb[iNWPacb],col="brown",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(160,40,"NEPac")
par(new=T)
plot(lonb[iSWPacb],latb[iSWPacb],col="purple",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(160,-55,"SWPac")
par(new=T)
plot(lonb[iSEPacb],latb[iSEPacb],col="green",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-120,-40,"SEPac")
par(new=T)
plot(lonb[iNIndb],latb[iNIndb],col="green",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(65,5,"NInd")
par(new=T)
plot(lonb[iSIndb],latb[iSIndb],col="pink",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(65,-45,"SInd")
par(new=T)
plot(lonb[iAAb],latb[iAAb],col="blue",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(5,-75,"Antarc")
par(new=T)
plot(lonb[iSWAtlb],latb[iSWAtlb],col="pink",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-40,-40,"SWAtl")
par(new=T)
plot(lonb[iSEAtlb],latb[iSEAtlb],col="brown",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(0,-15,"SEAtl")
par(new=T)
plot(lonb[iEEqAtlb],latb[iEEqAtlb],col="red",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-35,20,"ETAtl")
par(new=T)
plot(lonb[iWEqAtlb],latb[iWEqAtlb],col="purple",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-65,15,"SWAtl")
par(new=T)
plot(lonb[iNWAtlb],latb[iNWAtlb],col="orange",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-55,30,"NWAtl")
par(new=T)
plot(lonb[iNEAtlb],latb[iNEAtlb],col="green",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(-12,46,"NEAtl")
par(new=T)
plot(lonb[iArctb],latb[iArctb],col="red",xlim=xlm,ylim=ylm,xlab="",ylab="",yaxt="n")
text(35,74,"Arct",col="black")