## Not run:library(geomapdata)library(MBA)## for interpolation####### set up topo datadata(fujitopo)##### set up map datadata('japmap', package='geomapdata')### target regionPLOC= list(LON=c(138.3152,139.0214),LAT=c(35.09047,35.57324))PLOC$x =PLOC$LON
PLOC$y =PLOC$LAT
#### set up projectionPROJ = setPROJ(type=2, LAT0=mean(PLOC$y), LON0=mean(PLOC$x))########## select data from the topo data internal to the target topotemp = list(lon=fujitopo$lon, lat= fujitopo$lat, z=fujitopo$z)#### project target A = GLOB.XY(PLOC$LAT , PLOC$LON , PROJ)####### select toposelectionflag = topotemp$lat>+PLOC$LAT[1]& topotemp$lat<=PLOC$LAT[2]&topotemp$lon>+PLOC$LON[1]& topotemp$lon<=PLOC$LON[2]### project topo data B = GLOB.XY( topotemp$lat[selectionflag],topotemp$lon[selectionflag], PROJ)### set up out put matrix:### xo = seq(from=range(A$x)[1], to=range(A$x)[2], length=200)### yo = seq(from=range(A$y)[1], to=range(A$y)[2], length=200)####### interpolation using akima### IZ = interp(x=B$x , y=B$y, z=topotemp$z[selectionflag] , xo=xo, yo=yo)DF = cbind(x=B$x , y=B$y , z=topotemp$z[selectionflag]) IZ = mba.surf(DF,200,200, extend=TRUE)$xyz.est
xo = IZ[[1]] yo = IZ[[2]]### image(IZ)####### underwater section UZ = IZ$z
UZ[IZ$z>=0]=NA#### above sea level AZ = IZ$z
AZ[IZ$z<=-.01]=NA#### create perimeter: perim= getGEOperim(PLOC$LON, PLOC$LAT, PROJ,50)### lats for tic marks: PLAT = pretty(PLOC$LAT) PLAT = c(min(PLOC$LAT),PLAT[PLAT>min(PLOC$LAT)& PLAT<max(PLOC$LAT)],max(PLOC$LAT))PLON = pretty(PLOC$LON)### main program: DOIMG =TRUEDOCONT =TRUEPNTS =NULLBASICTOPOMAP(xo, yo , DOIMG, DOCONT, UZ, AZ, IZ, perim, PLAT, PLON,PROJ=PROJ, pnts=NULL, GRIDcol=NULL)### add in the map information plotGEOmapXY(japmap, LIM=c(PLOC$LON[1], PLOC$LAT[1],PLOC$LON[2],PLOC$LAT[2]), PROJ=PROJ, add=TRUE)## End(Not run)