library(GEOmap)
library(ETOPO)
library(WORLDMAP2)
library(geomapdata)



data(world2)
data(worldmap)

data(ETOPO5)
data(ETOPO2)


data(USAmap)

plotGEOmap(USAmap)
L = locator(2)


L=list()
L$x=c(232.234417898,252.135508975)
L$y=c(29.6187442109,44.0562015812)

 PROJ = setPROJ(type=1, LAT0=0 , LON0=0 )

PLOC=list(lon=L$x ,lat=L$y )



nn = names(PLOC)


ilon = grep("lon", nn, ignore.case = TRUE)
ilat = grep("lat", nn, ignore.case = TRUE)
A = list(lat=PLOC[[ilat[1]]], lon=PLOC[[ilon[1]]], LAT=PLOC[[ilat[1]]], LON=fmod(PLOC[[ilon[1]]], 360) )

ALOC=A

TOPO = ETOPO2


top2 = subsetTOPO(TOPO, ALOC)
d1 = dim(top2$z)

top2$z = top2$z[  , d1[2]:1 ]

image(top2, col=terrain.colors(100) )
calcol = settopocol()

Cmat  = TOPOCOL(top2, calcol$calcol)
Dcol  = attr(Cmat, 'Dcol')

X = GLOB.XY(top2$y[1], fmod(top2$x, 360) , PROJ)

DX = X$x[2]-X$x[1]

Y = GLOB.XY(top2$y, top2$x[1], PROJ)

 any(diff(X$x)<=0)
 any(diff(Y$y)<=0)



NEWimage = list(x=X$x, y=Y$y, z=top2$z/1000)


#########    jpostscript("BIGmap")

par(mai=c(.2, .2, .2,.2))



 PMAT = persp(NEWimage$x, NEWimage$y, NEWimage$z, theta = 0, phi = 90, r=4000,
     col=Cmat[1:(Dcol[1]-1), 1:(Dcol[2]-1)] , scale = FALSE,
           ltheta = 120, lphi=75, shade = 0.5, border = NA, expand=200, box = FALSE )


sqrTICXY(NEWimage, PROJ, side=1, LLgrid=FALSE, col=grey(.7), PMAT=PMAT )
sqrTICXY(NEWimage, PROJ, side=2, LLgrid=FALSE, col=grey(.7), PMAT=PMAT )


plotGEOmapXY(world2, PROJ =PROJ, PMAT =PMAT, add = TRUE)



dev.off()