####################################### load packages and high resolution data #### get GEOmap from CRAN #### get WORLDMAP2 ETOPO RJML from http://www.unc.edu/~leesj/FETCH/GRAB/RPACKAGES/ library(GEOmap) library(WORLDMAP2) library(ETOPO) library(RJML) data(world2) data(ETOPO2) #################################################### #################################################### #################################################### ####################################### ####### ##### ####### ####### # # # # ##### ####### # # # # # # # ## # # # ####### # # # # # # # # # # # ####### # # ####### # ####### # # # # # ####### # # # # # # # # # # # ####### # # # # # # # # ## # # ####### ##### ####### ####### ###### # # # # ##### ############################ this is extra stuff for the map - ############################ the user can add or modify this little function iceland.extra<-function(PLOC, PROJ, PMAT) { library(WORLDMAP2) data(worldvolcs) i1 = grep("HENGIL", worldvolcs$name) i2 = grep("KRAFLA", worldvolcs$name) LLin=list(tit=worldvolcs$name[c(i1, i2)], lat=worldvolcs$lat[c(i1, i2)], lon=worldvolcs$lon[c(i1, i2)]) vees = GLOB.XY( LLin$lat, LLin$lon , PROJ) ########### JML = vees if(!is.null(PMAT)) { JML = trans3d(vees$x, vees$y, rep(0, length(vees$y)), PMAT$PMAT) } else { JML = vees } points(JML$x, JML$y, pch=2, cex=2) pos =rep(3, length=length(LLin$tit)) text(JML$x, JML$y, labels=LLin$tit, pos=pos, font=2) } LL = c(63, 67, 333.4 , 348.02 ) MORE = 'iceland.extra' iclandlocales=list() iclandlocales$'lon'=c(-16.779397,-21.22646) iclandlocales$'lat'=c(65.70905,64.06038) iclandlocales$'name'=c("Krafla", "Hengil") iclandlocales$'txt'=c("Krafla Geothermal Field, Iceland", "Hengil Geothermal Field, Iceland") iclandlocales$'img'=c("krafla_13_small.jpg", "hengill_small_1.JPG") iclandlocales$'URL'=c("http://www.unc.edu/~leesj/Iceland/SIFKRA.htm", "http://www.unc.edu/~leesj/Hengill_WEB.htm") LAT1 = as.numeric(LL[1]) LAT2 = as.numeric(LL[2]) LON1 = as.numeric(LL[3]) LON2 = as.numeric(LL[4]) gLAT = c(LAT1, LAT2 ) gLON = c(LON1, LON2 ) PROJ <- setPROJ(2, LAT0 = mean(gLAT) , LON0 = mean(gLON) ) PLOC <- list( LON=gLON, LAT=gLAT) SMAP = SELGEOmap(world2, ncut = 3, acut = c(10, Inf), proj = PROJ, LIM = c(PLOC$LON[1], PLOC$LAT[1], PLOC$LON[2], PLOC$LAT[2])) dev.new() plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=FALSE, PMAT=NULL, xlab="km", ylab="km") ############# now check the added map stuff from the function set above: iceland.extra(PLOC, PROJ, NULL) ##################################################################### ############### check this function without using the makePIXmap function ############### calcol = settopocol() PMAT = GEOTOPO(ETOPO2, PLOC , PROJ, calcol$calcol) AMAT = persp(PMAT$xo, PMAT$yo, PMAT$IZ$z, theta = 0, phi = 90, r=4000, col=PMAT$Mollist[1:(PMAT$dMOL[1]-1), 1:(PMAT$dMOL[2]-1)] , scale = FALSE, ltheta = 90, lphi=30, shade = 1, border = NA, expand=0.004, box = FALSE ) plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=TRUE, PMAT=PMAT$PMAT) ############# now check the added map stuff from the function set above iceland.extra(PLOC, PROJ, PMAT) rxy = getGEOperim(range(PLOC$LON), range(PLOC$LAT), PROJ, 10) tem = trans3d(rxy$x, rxy$y, rep(0, length(rxy$y)), PMAT$PMAT) antipolygon(tem$x, tem$y, col = "white") LONTIX = pretty(range(PLOC$LON)) LONTIX = LONTIX[LONTIX >= min(PLOC$LON) & LONTIX <= max(PLOC$LON)] LATTIX = pretty(range(PLOC$LAT)) LATTIX = LATTIX[LATTIX >= min(PLOC$LAT) & LATTIX <= max(PLOC$LAT)] TIX = list(lons = c(PLOC$LON[1], LONTIX, PLOC$LON[2]), lats = c(PLOC$LAT[1], LATTIX, PLOC$LAT[2])) ticx = diff(PLOC$LAT) * 0.01 ticy = diff(PLOC$LON) * 0.01 addLLXY(TIX$lats, TIX$lons, PROJ = PROJ, PMAT = PMAT$PMAT, col = rgb(0, 0, 0), GRIDcol = NULL, GRID = FALSE, LABS = 0, TICS = c(ticx, ticy)) #################################################### #################################################### #################################################### ######## # # ####### ####### ####### # ######## # # # # # # # # # ######## # # # # # # # # # ######## ## # # ####### ####### ####### ######## # # # # # ## # # # ######## # # # # # ## # # # ######## # ## ####### # ## ####### # # lonKR = c(123,133) latKR = c(33,44) LL = c(latKR, lonKR) LAT1 = as.numeric(LL[1]) LAT2 = as.numeric(LL[2]) LON1 = as.numeric(LL[3]) LON2 = as.numeric(LL[4]) gLAT = c(LAT1, LAT2 ) gLON = c(LON1, LON2 ) PROJ <- setPROJ(2, LAT0 = mean(gLAT) , LON0 = mean(gLON) ) PLOC <- list( LON=gLON, LAT=gLAT) SMAP = SELGEOmap(world2, ncut = 3, acut = c(10, Inf), proj = PROJ, LIM = c(PLOC$LON[1], PLOC$LAT[1], PLOC$LON[2], PLOC$LAT[2])) dev.new() plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=FALSE, PMAT=NULL, xlab="km", ylab="km") PMAT = GEOTOPO(ETOPO2, PLOC , PROJ, calcol$calcol) AMAT = persp(PMAT$xo, PMAT$yo, PMAT$IZ$z, theta = 0, phi = 90, r=4000, col=PMAT$Mollist[1:(PMAT$dMOL[1]-1), 1:(PMAT$dMOL[2]-1)] , scale = FALSE, ltheta = 90, lphi=30, shade = 1, border = NA, expand=0.004, box = FALSE ) plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=TRUE, PMAT=PMAT$PMAT) ############# now check the added map stuff from the function set above: rxy = getGEOperim(range(PLOC$LON), range(PLOC$LAT), PROJ, 10) tem = trans3d(rxy$x, rxy$y, rep(0, length(rxy$y)), PMAT$PMAT) antipolygon(tem$x, tem$y, col = "white") LONTIX = pretty(range(PLOC$LON)) LONTIX = LONTIX[LONTIX >= min(PLOC$LON) & LONTIX <= max(PLOC$LON)] LATTIX = pretty(range(PLOC$LAT)) LATTIX = LATTIX[LATTIX >= min(PLOC$LAT) & LATTIX <= max(PLOC$LAT)] TIX = list(lons = c(PLOC$LON[1], LONTIX, PLOC$LON[2]), lats = c(PLOC$LAT[1], LATTIX, PLOC$LAT[2])) ticx = diff(PLOC$LAT) * 0.01 ticy = diff(PLOC$LON) * 0.01 addLLXY(TIX$lats, TIX$lons, PROJ = PROJ, PMAT = PMAT$PMAT, col = rgb(0, 0, 0), GRIDcol = NULL, GRID = FALSE, LABS = 0, TICS = c(ticx, ticy)) #################################################### #################################################### #################################################### #### ####### ##### ####### # # # # #### ## ## # # # # # # # #### ## # # # # # # # #### ## # # ####### # # #### ## # # # # # # #### ## ## # # # # # # #### ####### ##### # # # ###### # LL = c(35.6298770230, 42.3340428297, 11.3941923118, 19.6742015957 ) LAT1 = as.numeric(LL[1]) LAT2 = as.numeric(LL[2]) LON1 = as.numeric(LL[3]) LON2 = as.numeric(LL[4]) gLAT = c(LAT1, LAT2 ) gLON = c(LON1, LON2 ) PROJ <- setPROJ(2, LAT0 = mean(gLAT) , LON0 = mean(gLON) ) PLOC <- list( LON=gLON, LAT=gLAT) SMAP = SELGEOmap(world2, ncut = 3, acut = c(10, Inf), proj = PROJ, LIM = c(PLOC$LON[1], PLOC$LAT[1], PLOC$LON[2], PLOC$LAT[2])) dev.new() plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=FALSE, PMAT=NULL, xlab="km", ylab="km") PMAT = GEOTOPO(ETOPO2, PLOC , PROJ, calcol$calcol) AMAT = persp(PMAT$xo, PMAT$yo, PMAT$IZ$z, theta = 0, phi = 90, r=4000, col=PMAT$Mollist[1:(PMAT$dMOL[1]-1), 1:(PMAT$dMOL[2]-1)] , scale = FALSE, ltheta = 90, lphi=30, shade = 1, border = NA, expand=0.004, box = FALSE ) plotGEOmapXY(SMAP, PROJ=PROJ, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2], PLOC$LAT[2] ) , add=TRUE, PMAT=PMAT$PMAT) ############# now check the added map stuff from the function set above: rxy = getGEOperim(range(PLOC$LON), range(PLOC$LAT), PROJ, 10) tem = trans3d(rxy$x, rxy$y, rep(0, length(rxy$y)), PMAT$PMAT) antipolygon(tem$x, tem$y, col = "white") LONTIX = pretty(range(PLOC$LON)) LONTIX = LONTIX[LONTIX >= min(PLOC$LON) & LONTIX <= max(PLOC$LON)] LATTIX = pretty(range(PLOC$LAT)) LATTIX = LATTIX[LATTIX >= min(PLOC$LAT) & LATTIX <= max(PLOC$LAT)] TIX = list(lons = c(PLOC$LON[1], LONTIX, PLOC$LON[2]), lats = c(PLOC$LAT[1], LATTIX, PLOC$LAT[2])) ticx = diff(PLOC$LAT) * 0.01 ticy = diff(PLOC$LON) * 0.01 addLLXY(TIX$lats, TIX$lons, PROJ = PROJ, PMAT = PMAT$PMAT, col = rgb(0, 0, 0), GRIDcol = NULL, GRID = FALSE, LABS = 0, TICS = c(ticx, ticy)) #################################################### #################################################### #################################################### ######## ##### # ####### # # # ######## # # # # # # # ## # ######## # # # # # # # # # # ######## # ####### ####### ####### # # # ######## # # # # # # # # # # ######## # # # # # # # # ## ######## #### # # # # # # # #################################################### japan.extra<-function(PLOC, PROJ, PMAT) { library(WORLDMAP2) data(worldvolcs) vlon = fmod(worldvolcs$lon, 360) vlat = worldvolcs$lat i1 = which(vlon> PLOC$LON[1] & vlon< PLOC$LON[2] & vlat> PLOC$LAT[1] & vlat< PLOC$LAT[2]) LLin=list(tit=worldvolcs$name[c(i1)], lat=worldvolcs$lat[c(i1)], lon=worldvolcs$lon[c(i1)]) vees = GLOB.XY( LLin$lat, LLin$lon , PROJ) ########### JML = vees if(!is.null(PMAT)) { JML = trans3d(vees$x, vees$y, rep(0, length(vees$y)), PMAT$PMAT) } else { JML = vees } points(JML$x, JML$y, pch=2, cex=.6, col='red') pos =rep(3, length=length(LLin$tit)) ## text(JML$x, JML$y, labels=LLin$tit, pos=pos, font=2) } #################################################### #################################################### #################################################### LL = c(30.2663646058,37.2713265518, 128.121750925,141.661881078 ) #### first time it takes a long time to extract info from the data base: jap1 = myGMT(LL) #### myGMT(LL, PROJ=jap1$PROJ, PLOC=jap1$PLOC, PMAT=jap1$PMAT, SMAP=jap1$SMAP , do.topo=TRUE, do.map=TRUE) myGMT(LL, PROJ=jap1$PROJ, PLOC=jap1$PLOC, PMAT=jap1$PMAT, SMAP=jap1$SMAP , do.topo=TRUE, do.map=TRUE, dev="") japan.extra(jap1$PLOC, jap1$PROJ, jap1$PMAT) myGMT(LL, PROJ=jap1$PROJ, PLOC=jap1$PLOC, PMAT=jap1$PMAT, SMAP=jap1$SMAP , do.topo=FALSE, do.map=TRUE, dev="") myGMT(LL, PROJ=jap1$PROJ, PLOC=jap1$PLOC, PMAT=jap1$PMAT, SMAP=jap1$SMAP , do.topo=FALSE, do.map=TRUE, dev="") japan.extra(jap1$PLOC, jap1$PROJ, NULL) #################################################### #################################################### #################################################### ### ####### # # ##### ####### ####### ####### # # ##### # ### # # # # # # # # # # # ## # # # # ### # # # # # # # # # # # # # # # # ### # ####### # # ####### # # ####### # # # # ####### ### # # # # # # # # # ## # # # # # # ### # # # # # # # # # # ## # ## # # # ### ####### # # ###### ##### # ####### # ## # # ##### # # LL = c( 31, 38, 235.5 , 248 ) cal1 = myGMT(LL) CAL.extra<-function(PLOC, PROJ, PMAT) { library(WORLDMAP2) data(eqs) vlon = fmod(eqs$lon, 360) vlat = eqs$lat i1 = which(vlon> PLOC$LON[1] & vlon< PLOC$LON[2] & vlat> PLOC$LAT[1] & vlat< PLOC$LAT[2]) LLin=list( lat=eqs$lat[c(i1)], lon=eqs$lon[c(i1)]) vees = GLOB.XY( LLin$lat, LLin$lon , PROJ) ########### JML = vees if(!is.null(PMAT)) { JML = trans3d(vees$x, vees$y, rep(0, length(vees$y)), PMAT$PMAT) } else { JML = vees } points(JML$x, JML$y, pch=".", cex=2, col='red') pos =rep(3, length=length(LLin$tit)) ## } myGMT(LL, PROJ=cal1$PROJ, PLOC=cal1$PLOC, PMAT=cal1$PMAT, SMAP=cal1$SMAP , do.topo=FALSE, do.map=TRUE, dev="") CAL.extra(cal1$PLOC, cal1$PROJ, NULL ) myGMT(LL, PROJ=cal1$PROJ, PLOC=cal1$PLOC, PMAT=cal1$PMAT, SMAP=cal1$SMAP , do.topo=TRUE, do.map=TRUE, dev="") CAL.extra(cal1$PLOC, cal1$PROJ, cal1$PMAT ) #################################################### #################################################### ####################################################