GEOmap (version 2.1)

BASICTOPOMAP: Basic Topogrpahy Map

Description

Basic Topogrpahy Map

Usage

BASICTOPOMAP(xo, yo, DOIMG, DOCONT, UZ, AZ, IZ, perim, PLAT, PLON, PROJ = PROJ, pnts = NULL, GRIDcol = NULL)

Arguments

xo

vector of x-coordinates

yo

vector of y-coordinates

DOIMG

logical, add image

DOCONT

logical, add contours

UZ

matrix of image values under sea level

AZ

matrix of image values above sea level

IZ

matrix of image values

perim

perimeter vectors

PLAT

latitudes for tic-marks

PLON

longitude for tic-marks

PROJ

projection list

pnts

points to add to plot

GRIDcol

color for grid

Value

• Graphical Side effects

Details

Image is processed prior to calling

See Also

DOTOPOMAPI, GEOTOPO

Examples

library(geomapdata)
library(MBA) ##  for interpolation
#######  set up topo data
data(fujitopo)
#####  set up map data
data(japmap)


###  target region
PLOC= list(LON=c(138.3152, 139.0214), 
LAT=c(35.09047, 35.57324))

PLOC$x =PLOC$LON
PLOC$y =PLOC$LAT



####  set up projection
PROJ = 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 topo
selectionflag = 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 = TRUE
DOCONT = TRUE
PNTS  = NULL

BASICTOPOMAP(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 )