LandSeaCol: Land and Sea Colors

Description

Color pixels with two palettes, one for land the other for sea.

Usage

LandSeaCol(IZ, coastmap, PROJ, calcol = NULL)

Arguments

list of x, y, z suitable for plotting with image or contour.

coastmap

coastal map from GEOmap

PROJ

projection list

calcol

color map for the land

Value

CmatMatrix of colors for each pixel
UZUnder water
AZAbove Sea Level

Details

The program uses closed polygons in the map list to separate the pixels into land versus sea. Sea is colored with a palette of blues, land is colored according to topographic color scheme extracted from palettes similar to GMT palettes.

All map and pixel coordinates are projected with the same projection parameters. calculations are done in XY coordinates.

Examples

Run this code

Lat.range = c(-10, 30)
Lon.range = c(65, 117)
######

########  load up the important libraries
 library(RFOC)
library(GEOmap)
 library(geomapdata)

 data(coastmap)
data(worldmap)
     data(ETOPO5)
 PLOC=list(LON=Lon.range,LAT=Lat.range,lon=Lon.range,lat=Lat.range,
                 x=Lon.range, y=Lat.range )


#####   set up topography colors
COLS = settopocol()

####  set the projection
 PROJ = setPROJ(type=2, LAT0=mean(PLOC$y) , LON0=mean(PLOC$x) )   ##   utm
 NK = 300
    
   ###   extract topography from the etopo5 data base in geomapdata
     JMAT = GEOTOPO(ETOPO5, PLOC, PROJ, COLS$calcol, npoints=NK)
#####   select relevant earthquakes

IZ = list(x=JMAT$xo, y=JMAT$yo, z=JMAT$IZ$z)

CMAT = LandSeaCol(IZ, coastmap, PROJ, calcol=NULL)

Mollist =CMAT$Cmat
dMol = attr(Mollist, "Dcol")

     ####  Under water
UZ = CMAT$UZ
    #####  above water
AZ = CMAT$AZ
    ####   blues for underwater:
   blues = shade.col(100, acol=as.vector(col2rgb("darkblue")/255)   , bcol= as.vector(col2rgb("paleturquoise")/255))

plot(x=range(IZ$x), y=range(IZ$y), type='n', asp=1, axes=FALSE, ann=FALSE)

image(x=IZ$x, y=IZ$y, z=(UZ), col=blues, add=TRUE)

image(x=IZ$x, y=IZ$y, z=(AZ), col=terrain.colors(100) , add=TRUE)

plotGEOmapXY(worldmap,
              LIM = c(Lon.range[1],Lat.range[1] ,Lon.range[2] ,Lat.range[2]),
              PROJ =PROJ,  add=TRUE  )