Package 'GEOmap'

Description

Topographic and Geologic Mapping

Details

Set of routines for making Map Projections (forward and inverse), Topographic Maps, Perspective plots, geological databases, interactive plotting and selection of focus regions.

Note

High level plotting:

BASICTOPOMAP DOTOPOMAPI geoLEGEND GEOsymbols locworld plotGEOmap plotGEOmapXY linesGEOmapXY rectGEOmapXY textGEOmapXY pointsGEOmapXY insideGEOmapXY plotUTM plotworldmap XSECDEM

PLOTTING:

circle addLLXY addTIX antipolygon zebra demcmap setXMCOL shade.col

Geological Map Symbols:

bcars faultdip faultperp horseshoe normalfault OverTurned perpen teeth thrust SynAnticline SSfault

Data manipulation:

getGEOmap boundGEOmap SELGEOmap geoarea GEOTOPO getGEOperim GETXprofile Lintersect LOCPOLIMAP pline selectPOLImap setplotmat SETPOLIMAP settopocol subsetTOPO

Misc:

getgreatarc ccw difflon DUMPLOC getsplineG inpoly inside PointsAlong polyintern

Projections:

setPROJ projtype GLOB.XY XY.GLOB MAPconstants GCLCFR lambert.cc.ll lambert.cc.xy lambert.ea.ll lambert.ea.xy lcgc merc.sphr.ll merc.sphr.xy utmbox utm.elps.ll utm.elps.xy utm.sphr.ll utm.sphr.xy stereo.sphr.ll stereo.sphr.xy equid.cyl.ll equid.cyl.xy

Author(s)

Jonathan M. Lees<jonathan.lees.edu> Maintainer:Jonathan M. Lees<[email protected]>

References

Snyder, John P., Map Projections- a working manual, USGS, Professional Paper, 1987.

Lees, J. M., Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences, Computers & Geosciences, 26, 7, 751-761, 2000.

See Also

RSEIS

Examples

################  projections
proj = setPROJ(type = 2, LAT0 =23, LON0 = 35)

### get lat-lon
LL = XY.GLOB(200, 300, proj)


##  find x-y again, should be the same
XY = GLOB.XY(LL$lat, LL$lon, proj)
XY
################
library(geomapdata)
data(worldmap)
  KAMlat = c(48.5,  65)
    KAMlon = c(150, 171)
    
    PLOC=list(LON=KAMlon,LAT=KAMlat)


    PLON = seq(from=KAMlon[1], to=KAMlon[2], by=2)
    PLAT = seq(from=KAMlat[1], to=KAMlat[2], by=2)
    
    proj = setPROJ(2, LON0=mean(KAMlon), LAT0=mean(KAMlat))

 xy = GLOB.XY(KAMlat,  KAMlon , proj)
 kbox=list(x=range(xy$x, na.rm=TRUE), y=range(xy$y, na.rm=TRUE))

 plot(kbox$x,kbox$y, type='n', axes=FALSE, xlab="", ylab="", asp=1)
   plotGEOmapXY(worldmap, LIM=c(KAMlon[1], KAMlat[1], KAMlon[2],
KAMlat[2]),  add=TRUE, PROJ=proj, axes=FALSE, xlab="", ylab="" )

sqrTICXY(kbox , proj, side=c(1,2,3,4), LLgrid=TRUE, col=grey(.7) )
title("Crude Map of Kamchatka")

Description

Add Lat-Lon points using projection

Usage

addLLXY(lats, lons, PROJ = PROJ, PMAT = NULL,
col = gray(0.7), GRID = TRUE, GRIDcol = 1, LABS = NULL,
LABcol = 1, BORDER = NULL, TICS = c(1, 1), xpd=TRUE)

Arguments

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

plotGEOmapXY, sqrTICXY

Examples

library(geomapdata)


data('fujitopo', package='geomapdata')
data('japmap', package='geomapdata')

PLOC=list(LON=range(c( japmap$STROKES$LON1,japmap$STROKES$LON2) ),
LAT=range(c( japmap$STROKES$LAT1,japmap$STROKES$LAT2) ))
PLOC$x = PLOC$LON
PLOC$y = PLOC$LAT

PROJ = setPROJ(type=2, LAT0=mean(PLOC$y) , LON0=mean(PLOC$x) )
isel1 = which( japmap$STROKES$code != "i" & japmap$STROKES$num>120 )


plotGEOmapXY(japmap, PROJ=PROJ,SEL=isel1,  add=FALSE, axes=FALSE, xlab="", ylab="")
A = PLOC

   PLAT =  pretty(A$LAT)
    PLAT = c(min(A$LAT),  PLAT[PLAT>min(A$LAT) & PLAT<max(A$LAT)],max(A$LAT)) 
  PLON  = pretty(A$LON)
        PLON = c(min(A$LON), PLON[PLON>min(A$LON) & PLON<max(A$LON)],
 max(A$LON))


addLLXY(PLAT,  PLON, PROJ=PROJ, LABS=TRUE, PMAT=NULL, TICS=c(.1,.1) )

###############

Description

Add Tic marks to map

Usage

addTIX(lats, lons, PROJ = list(), PMAT = NULL,
col = gray(0.7), TICS = c(1, 1), OUTER = TRUE,
sides = c(1, 2, 3, 4))

Arguments

Details

attempts to make correct default values

Value

Graphical Side Effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

addLLXY

Examples

##########3  this program is run internally


PLOC=list(LON=c(137.008, 141.000),
LAT=c(34.000, 36.992),
x=c(137.008, 141.000),
y=c(34.000, 36.992))

PROJ = setPROJ(type=2, LAT0=mean(PLOC$y) , LON0=mean(PLOC$x) )

gxy = GLOB.XY(PLOC$LAT, PLOC$LON, PROJ)

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)

PLON = c(min(PLOC$LON), PLON[PLON>min(PLOC$LON)&PLON<max(PLOC$LON)], max(PLOC$LON))
     

plot(gxy$x, gxy$y,  asp=TRUE)

 addTIX(PLAT, PLON, PMAT=NULL, col='red', TICS=c(.1,.1), PROJ=PROJ)

Description

Calculate points along a great arc

Usage

along.great(phi1, lam0, c, Az)

Arguments

Details

All input and output is radians

Value

List:

Author(s)

Jonathan M. Lees<[email protected]>

Examples

lat1 <- 48.856578
lon1 <- 2.351828	

A = along.great(lat1*pi/180, lon1*pi/180, 50*pi/180, -63*pi/180)

lat=A$phi*180/pi
lon = A$lam*180/pi

Description

Fill a plot with a color outside the confines of a polygon.

Usage

antipolygon(x, y, col = 0, corner=1, pct=.4)

Arguments

Details

antipolygon uses par("usr") to determine the external bounds of plotting region. Corners are labels from bottom left counter-clockwise, 1-4.

Value

List:

Used for graphical side effect

Note

If the figure is resized after plotting, filling may not appear correct.

Author(s)

Jonathan M. Lees <[email protected]>

See Also

polygon, par

Examples

set.seed(2018)
x = runif(100)
y = runif(100)

#########  some data points to plot:

plot(x,y)
###########   create polygon:
pp =list(x=c(0.231,0.316,0.169,0.343,0.311,0.484,0.757,
         0.555,0.800,0.563,0.427,0.412,0.203),
     y=c(0.774,0.622,0.401,0.386,0.138,0.312,0.200,0.459,
        0.658,0.624,0.954,0.686,0.813))

polygon(pp)

antipolygon(x=pp$x, y=pp$y,col='blue')
####  where as this does not look so good
plot(x,y)
antipolygon(x=pp$x, y=pp$y,col='blue', corner=2)

Description

Basic Topogrpahy Map

Usage

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

Arguments

Details

Image is processed prior to calling

Value

Graphical Side effects

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

DOTOPOMAPI, GEOTOPO

Examples

## Not run: 


library(geomapdata)
library(MBA) ##  for interpolation
#######  set up topo data
data(fujitopo)
#####  set up map data
data('japmap', package='geomapdata' )


###  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 )


## End(Not run)

Description

Add Box Cars to a line.

Usage

bcars(x, y, h1 = 1, h2 = 0.3, rot, col = "black", border = "black")

Arguments

Details

Used for plotting detachment faults in USGS format.

Value

Graphical Side effects

Author(s)

Jonathan M. Lees<[email protected]>

Examples

G=list()
G$x=c(-1.0960,-0.9942,-0.8909,-0.7846,-0.6738,-0.5570,-0.4657,-0.3709,
-0.2734,-0.1740,-0.0734, 0.0246, 0.1218, 0.2169, 0.3086, 0.3956, 0.4641, 
0.5293, 0.5919, 0.6530, 0.7131)
G$y=c(-0.72392,-0.62145,-0.52135,-0.42599,-0.33774,-0.25896,-0.20759,
-0.16160,-0.11981,-0.08105,-0.04414,-0.00885, 0.02774, 0.06759, 0.11262, 
0.16480, 0.21487, 0.27001, 0.32895, 0.39044, 0.45319)



 g = PointsAlong(G$x, G$y, N=6)

 sk = 3

###############
plot(G$x, G$y, type='n',asp=1, axes=FALSE, xlab='', ylab='')

 lines(G$x,G$y,col='blue')
 bcars(g$x,g$y,h1=sk,h2=sk*.5, rot=g$rot , col='blue')


###############
plot(G$x, G$y, type='n',asp=1, axes=FALSE, xlab='', ylab='')
 lines(G$x,G$y,col='blue')
 bcars(g$x,g$y,h1=sk,h2=sk*.5, rot=g$rot , col=NA, border='blue')

Description

Given a GEOmap strucutre, set the bounds for the strokes.

Usage

boundGEOmap(MAP, NEGLON = FALSE, projtype = 2)

Arguments

Details

Used to rectify a new map after reading in from ascii file. Can take GMT map ascii map files and convert to GEOmap.

Value

List structure:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

worldmap

Examples

library(geomapdata)
data(worldmap)
worldmap = boundGEOmap(worldmap)

Description

Check for counter-clockwise orientation for polygons. Positive is counterclockwise.

Usage

CCcheck(Z)

Arguments

Details

Uses sign of the area of the polygon to determine polarity.

Value

Note

Based on the idea calculated area of a polygon.

Author(s)

Jonathan M. Lees<[email protected]>

Examples

Y=list()
Y$x=c(170,175,184,191,194,190,177,166,162,164)
Y$y=c(-54,-60,-60,-50,-26,8,34,37,10,-15)

plot(c(160, 200),c(-85, 85), type='n')
points(Y)
lines(Y)

CCcheck(Y)


Z = list(x=rev(Y$x), y=rev(Y$y))

CCcheck(Z)

Description

Used for determining if points are in polygons.

Usage

ccw(p0, p1, p2)

Arguments

Value

returns 1 or 0 depending on position of points

Author(s)

Jonathan M. Lees <[email protected]>

See Also

Lintersect

Examples

l1 = list(p1=list(x=0, y=0), p2=list(x=1,y=1))
l2 = list(p1=list(x=6, y=4), p2=list(x=-1,y=-12))

ccw(l1$p1, l1$p2, l2$p1)

Description

Global Maps of Coast

Usage

data(coastmap)

Format

List structure:

STROKES

list(nam, num, index, col, style, code, LAT1, LAT2, LON1, LON2)

POINTS

list(lat, lon)

PROJ

list(type, LAT0, LON0, LAT1, LAT2, LATS, LONS, DLAT, DLON, FE, FN, name)

Details

This map list is used for filling in coastal lines for global maps. The style=3 is for filling in polygons. The strokes are named for easier access to particular parts ofthe globe. Asia and Africa are one stroke, as are North and South America. there are currently three codes: C=major coast, c=smaller coasts, L=interior lakes.

Examples

data(coastmap)
#######  see the codes:
unique(coastmap$STROKES$code)
#########  see the different names:
unique(coastmap$STROKES$nam)

#########  change the colors based on code
coastmap$STROKES$col[coastmap$STROKES$code=="C" ] = rgb(1, .6, .6)
coastmap$STROKES$col[coastmap$STROKES$code=="c" ] = rgb(1, .9, .9)
coastmap$STROKES$col[coastmap$STROKES$code=="L" ] = rgb(.6, .6, 1)

plotGEOmap(coastmap , border='black' , add=FALSE, xaxs='i')





##

Description

Draw acircular arc from angle 1 to angle 2 at a given location.

Usage

darc(rad = 1, ang1 = 0, ang2 = 360, x1 = 0, y1 = 0, n = 1)

Arguments

Details

If angle1 > angle2 arc is drawn in opposite direction

Value

list(x,y)

Author(s)

Jonathan M. Lees<[email protected]>

Examples

plot(c(0,1), c(0,1), type='n', ann=FALSE, asp=1)

A = darc(.3, 23, 47, .5, .5, n=1)

lines(A$x, A$y)

Description

Return a small data base of Datum values for use in UTM projections.

Usage

DATUMinfo()

Details

The function just return a list with the relavent information.

Value

List:

Author(s)

Jonathan M. Lees<[email protected]>

References

stevedutch.net

See Also

UTM.xy, UTM.ll, setPROJ

Examples

h = DATUMinfo()
data.frame(h)

Description

create a color map from a DEM (Digital Elevation Map)

Usage

demcmap(ZTOPO, n = 100, ccol = NULL)

Arguments

Value

vector of rgb colors

Author(s)

Jonathan M. Lees<[email protected]>

See Also

rgb, settopocol

Description

Difference between Longitudes

Usage

difflon(LON1, LON2)

Arguments

Details

takes into account crossing the zero longitude

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

difflon( 34 , 67)

###  here we cross the zero line
difflon( 344 , 67)

Description

Calculate distance, Azimuth and Back-Azimuth from two points on Globe.

Usage

distaz(olat, olon, tlat, tlon)

Arguments

Details

Program is set up for one origin (olat, olon) pair and many target (tlat, tlon) pairs given as vectors.

If multiple olat and olon are given, the program returns a list of outputs for each.

If olat or any tlat is greater than 90 or less than -90 NA is returned and error flag is 0.

If any tlat and tlon is equal to olat and olon, the points are coincident. In that case the distances are set to zero, but the az and baz are NA, and the error flag is set to 0.

Value

List:

Author(s)

Jonathan M. Lees<[email protected]>

See Also

along.great, getgreatarc

Examples

####  one point
d = distaz(12, 23, -32,    -65)
d

####  many random target points
org = c(80.222, -100.940)
targ = cbind(runif(10, 10, 50), runif(10, 20, 100))


distaz(org[1], org[2], targ[,1], targ[,2])

############  if origin and target are identical
#####  the distance is zero, but the az and baz are not defined
distaz(80.222, -100.940, 80.222, -100.940)


########################   set one of the targets equal to the origin
targ[7,1] = org[1]
targ[7,2] = org[2]

distaz(org[1], org[2], targ[,1], targ[,2])

####  put in erroneous latitude data

targ[3,1] = -91.3


distaz(org[1], org[2], targ[,1], targ[,2])

Description

Convert decimal degrees to degree, minutes, seconds

Usage

dms(d1)

Arguments

Value

list

Author(s)

Jonathan M. Lees<[email protected]>

Examples

dms(33.12345)

H = dms(-91.8765)

print(H)

newH = H$d+H$m/60+H$s/3600
print(newH)

Description

For saving vectors to a file after the locator function has been executed.

Usage

DUMPLOC(zloc, dig = 12)

Arguments

Value

Side effects: print to screen

Author(s)

Jonathan M. Lees<[email protected]>

Examples

G=list()
G$x=c(-1.0960,-0.9942,-0.8909,-0.7846,-0.6738,-0.5570,-0.4657,-0.3709,
-0.2734,-0.1740,-0.0734, 0.0246, 0.1218, 0.2169, 0.3086, 0.3956, 0.4641, 
0.5293, 0.5919, 0.6530, 0.7131)
G$y=c(-0.72392,-0.62145,-0.52135,-0.42599,-0.33774,-0.25896,-0.20759,
-0.16160,-0.11981,-0.08105,-0.04414,-0.00885, 0.02774, 0.06759, 0.11262, 
0.16480, 0.21487, 0.27001, 0.32895, 0.39044, 0.45319)

g = PointsAlong(G$x, G$y, N=3)
DUMPLOC(g, dig = 5)

Description

Global Earthquake catalog locations from Engdahl, et al.

Usage

data(EHB.LLZ)

Format

lat

Latitude

lon

Longitude

z

depth in km

Source

Data is extrcted from an earthquake data base of relocated events provided by Robert Engdahl.

References

Engdahl, E. R., R. D. van der Hilst, S. H. Kirby, G. Ekstrom, K. M. Shedlock, and A. F. Sheehan (1998), A global survey of slab structures and internal processes using a combined data base of high-resolution earthquake hypocenters, tomographic images and focal mechanism data, Seismol. Res. Lett., 69, 153-154.

Examples

data(EHB.LLZ)
## maybe str(EHB.LLZ) ; plot(EHB.LLZ) ...

Description

Ellipsoidal Distance given Latitude and Longitude

Usage

Ellipsoidal.Distance(olat, olon, tlat, tlon, a = 6378137, b = 6356752.314, tol=10^(-12))

Arguments

Details

Uses Vincenty's formulation to calculate the distance along a great circle on an ellipsoidal body.

If a and be are not provided, they are set by default to a=6378137.0 , b=6356752.314, the WGS-84 standard.

Only one pair of (olat, olon) and (tlat, tlon) can be given at a time. The program is not vectorized.

Quoting from the wiki page this algorithm was extracted from:

"Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of an spheroid, developed by Thaddeus Vincenty in 1975. They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth.

The first (direct) method computes the location of a point which is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (.020 sec) on the Earth ellipsoid"

Value

list

Note

Latitudes >90 and < -90 are not allowed. NA's are returned.

If points are identical, a distance of zero is returned and NA for the azimuths. If there is some problems with convergence or division by zero, NA's are returned and error message is printed.

A couple of known cases that do not work are, e.g.: (olat=0; olon=0; tlat=0; tlon=-180) and (olat=0; olon=0; tlat=0; tlon=180). They will return NA's to avoid division by zero.

I am not sure how to deal with these cases yet.

The reverse azimuth is the angle from the meridian on the target point to the great circle from the origin to the target (as far as I can tell). If distaz and Ellipsoidal.Distance are compared, they give the same azimuth, and the absolute angles of baz (from distaz) and revaz (from Ellipsoidal.Distance) will add to 180 degrees.

Author(s)

Jonathan M. Lees<[email protected]>

References

http://en.wikipedia.org/wiki/Vincenty%27s_formulae

Vincenty, T. (April 1975). Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Survey Review XXIII (misprinted as XXII) (176): 88.201393. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf. Retrieved 2009-07-11.

See Also

distaz

Examples

####    compare to spheroidal calculation distaz
####


R.MAPK = 6378.2064
N =20

OUT = list(dadist=0, ed2dist=0, ed1dist=0, dif2=0, dif1=0, pct1=0)
for( i in 1:N)
  {

    olat = runif(1, -90, 90)
    olon = runif(1, 0, 180)

     tlat = runif(1, -90, 90)
    tlon = runif(1, 0, 180)

##########  older spherical calculation
    da = distaz(olat, olon, tlat, tlon)
#####  ed1 = elliposidal earth
    ed1 = Ellipsoidal.Distance(olat, olon, tlat, tlon)
#####   ed2 spherical earth using
############      ellipsoidal calculations, compare with
distaz
   ed2 =  Ellipsoidal.Distance(olat, olon, tlat, tlon, a=R.MAPK*1000, b=R.MAPK*1000)

    dif1 =  da$dist-ed1$dis
    dif2 =  da$dist-ed2$dis

    pct1 = 100*dif1/ed1$dist

##############   OUT = format( c(da$dist, ed2$dist, ed1$dist, dif2, dif1, pct1) , digits=10)

    OUT$dadist[i] =da$dist
    OUT$ed2dist[i] =ed2$dist
OUT$ed1dist[i]=ed1$dist
OUT$dif2[i]= dif2
OUT$dif1[i]=dif1
OUT$pct1[i]=pct1

###cat(paste(collapse=" ", OUT), sep="\n")

  }


print( data.frame(OUT) )



###############   some extreme cases can cause problems
#######  here compare  Ellipsoidal.Distance with spherical program distaz

Alat = c(90,   90,  90,   90,  45,  45,  45,  45,   0,    0,    0, 0)
Alon = c(180, 180,-180, -180,  45,  45,  45,  45,   0,    0,    0, 0)
Blat = c(-90, -45,   0,   45, -45,  0,    0,  -80,  45,   0,    0, 0)
Blon = c(180,-180, 180,    0, -45,  0, -180,  100, -60, -180, 180, 0)


BOUT = list(olat=0, olon=0, tlat=0, tlon=0, dadist=0, ed2dist=0, daaz=0, ed2az=0, dabaz=0, ed2baz=0)

R.MAPK = 6378.2064
for(i in 1:length(Alat))
{

  olat = Alat[i]
  olon  = Alon[i]
  tlat  = Blat[i]
  tlon  = Blon[i]

 da = distaz(olat, olon, tlat, tlon)
  ed2 = Ellipsoidal.Distance(olat, olon, tlat, tlon, a=R.MAPK*1000, b=R.MAPK*1000)
 cat(paste("i=", i), sep="\n")
 

 BOUT$olon[i] =olon
 BOUT$olat[i] =olat
 BOUT$tlat[i] =tlat
 BOUT$tlon[i] =tlon


      BOUT$dadist[i] =da$dist
   BOUT$ed2dist[i] =ed2$dist

BOUT$daaz[i]= da$az
BOUT$dabaz[i]= da$baz

BOUT$ed2az[i]= ed2$az
BOUT$ed2baz[i]=  ed2$revaz


}

print(data.frame(BOUT))

Description

Extract a set of eathquakes in swath along a cross sectional line

Usage

eqswath(x, y, z, L, width = 1, PROJ = NULL)

Arguments

Details

All units should be the same.

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

XSECwin, XSECEQ

Examples

#############  create data
x = runif(100, 1, 100)  
y = runif(100, 1, 100) 
z = runif(100, 1, 10) 
plot(x,y, asp=1)
## L = locator()

L=list()
L$x=c( 5.42328560757,64.62879777806)
L$y=c(89.843266449785,-0.174423911329)

J = eqswath(x, y, z, L, width = 10, PROJ = NULL)

##########   show box:
plot(x,y, asp=1)
lines(J$InvBox$x, J$InvBox$y)


############  show cross section with events plotted
plot(J$r, -J$depth)

Description

Select sections of a MAP-list structure based on stroke index

Usage

ExcludeGEOmap(MAP, SEL, INOUT = "out")

Arguments

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

getGEOmap, plotGEOmap, SELGEOmap, boundGEOmap

Examples

data(coastmap)

###  extract (include)  the first 6 strokes from world map


A1 = ExcludeGEOmap(coastmap, 1:6, INOUT="in")
print(A1$STROKES$nam)

Description

Calculate an expanded bounding region based on a percent of the existing boundaries

Usage

expandbound(g, pct = 0.1)

Arguments

Details

uses the range of the exising vector to estimate the expanded bound

Value

vector, new range

Author(s)

Jonathan M. Lees<[email protected]>

Examples

i = 5:10
exi = expandbound(i, pct = 0.1)
range(i)
range(exi)

Description

Explode a set of points away from a center point

Usage

explode(fxy, dixplo=1, mult=1, cenx=0, ceny=0, PLOT=FALSE)

Arguments

Details

If cenx and ceny is missing it is assumed to be the mean of the coordinates. Program calculates the new locations radiating away from the central point. No protection against overlapping symbols is included.

Value

list of new x,y values

Author(s)

Jonathan M. Lees<[email protected]>

See Also

ExplodeSymbols, NoOverlap

Examples

############  random data
x = rnorm(20)
y = rnorm(20)

NEW = explode(list(x=x,y=y), dixplo =1)

plot(range(c(x,NEW$x)), range(c(y,NEW$y)), asp=1, type='n')
segments(x,y,NEW$x, NEW$y)
points(x,y, pch=3, col='red')
points(NEW$x, NEW$y, pch=6, col='blue', cex=2)


###  try a larger radius:
NEW2 = explode(list(x=x,y=y), dixplo =1.3)
points(NEW2$x, NEW2$y, pch=7, col='brown', cex=2, xpd=TRUE)
arrows(NEW$x, NEW$y,NEW2$x, NEW2$y, col='green' )



#####   try with a different center
cenx=-1; ceny=-1
NEW = explode(list(x=x,y=y), dixplo =1, cenx=cenx, ceny=ceny)
plot(range(c(x,NEW$x)), range(c(y,NEW$y)), asp=1, type='n')
points(x,y, pch=3, col='red')
segments(x,y,NEW$x, NEW$y)
points(NEW$x, NEW$y, pch=6, col='blue', cex=2)
points(cenx, ceny, pch=8, col='purple')
text(cenx, ceny, labels="Center Point", pos=1)

Description

Interactive program for redistributing symbols for later plotting. Used for Focal Mechanisms.

Usage

ExplodeSymbols(XY, fsiz = 1, STARTXY = NULL, MAP = NULL)

Arguments

Details

The program is interactive. It starts by plotting the points as symbols. A number of buttons are provided for exploding the points semi automatically. To move each point click near its current point, then click at the destination followed by a click on the HAND button. several symbols can be moved at the same time.

You must click on the screen and on the buttons to get this code working - the program will not work in batch mode or run as a script You click in the active screen area and then press a button on top (or bottom) - the button takes your clicks and does something Here are some hints:

Buttons:Buttons appear on top and bottom of the plotting region.

HAND: If you want to move only one symbol (focal mech) click near it and then click where you want it to go. Then click the HAND button You may click several at once, but for each click oin a symbol there has to be a click somewhere to relocate it. (i.e. there must be an even number of clicks on the screen before hitting the HAND button)

SEL: If you want to explode several symbols at once, first select them: click lower left, then upper right of rectangle enclosing the selection. Once a selection is made it remains active until another selection is made so you can keep changing the radius and center for different explosions Then click CIRC.

RECT Choose a rectangle (lower left and upper right), then click RECT for an explosion

RECT2 After selecting, choose a center and a distance. symbols will be moved to a rectangular perimeter defined by the two points

CIRC After selection, click once for the circle center, and a second time for the radius, then click CIRC

LINE After selection,will explode the events away from a line, a given distance away. The line is given by 2 points and the distance by a third perpendicular distance.

Value

list of new x,y values

Note

For now the map is given in lat-lon coordinates- the same as the points being moved. There is no map projection used.

Author(s)

Jonathan M. Lees<[email protected]>

See Also

rekt2line

Examples

## Not run: 
F1 = list(x=rnorm(43), y=rnorm(43))
SMXY = ExplodeSymbols(F1, 0.03)




## End(Not run)

Description

Show Fault dip

Usage

faultdip(x, y, rot = 0, h = 1, lab = "")

Arguments

Value

Graphical Side effect

Author(s)

Jonathan M. Lees<[email protected]>

See Also

perpen, PointsAlong, getsplineG

Examples

ff=list()
ff$x=c(-4.850,-4.700,-3.934,-2.528, 0.603, 2.647, 3.861, 2.626)
ff$y=c(-4.045,-2.087,-0.710, 0.172, 1.291, 2.087,-0.753,-4.131)

G =getsplineG(ff$x, ff$y, kdiv=20)
g = PointsAlong(G$x, G$y, N=5)


plot(c(-5,5), c(-5,5), asp=1, type='n' )
lines(G)

angs = 180*atan(g$rot$sn/g$rot$cs)/pi
faultdip(g$x , g$y , rot=angs, h=.5, lab='')

Description

Draw perpendicular marks on fault trace

Usage

faultperp(x, y, N = 20, endtol = 0.1, h = 1, col = "black")

Arguments

Value

Graphical Side effect

Author(s)

Jonathan M. Lees<[email protected]>

See Also

OverTurned

Examples

ff=list()
ff$x=c(-4.850,-4.700,-3.934,-2.528, 0.603, 2.647, 3.861, 2.626)
ff$y=c(-4.045,-2.087,-0.710, 0.172, 1.291, 2.087,-0.753,-4.131)

G =getsplineG(ff$x, ff$y, kdiv=20)
g = PointsAlong(G$x, G$y, N=5)


plot(c(-5,5), c(-5,5), asp=1, type='n' )
lines(G)

faultperp(G$x, G$y, N = 10, endtol = 0.1, h = .3, col = "black")

Description

Correct wrapping for GEOmaps

Usage

fixCoastwrap(Z, maxdis = 100)

Arguments

Details

Based on mapswrap program

Value

List:

Author(s)

Jonathan M. Lees<[email protected]>

Examples

data(coastmap)
SEL = which(coastmap$STROKES$nam=="AFROASIA")

A = ExcludeGEOmap(coastmap, SEL, INOUT="in")

plot(A$POINTS$lon, A$POINTS$lat, type='n')

points(A$POINTS$lon, A$POINTS$lat, pch='.')

######  note that the map wraps around.

B = fixCoastwrap(list(x=A$POINTS$lon, y=A$POINTS$lat), 100)
 which(is.na(B$x))



lines(B)


polygon(B, col=rgb(.8,1, .8))

Description

OLD projection sometimes used in Lees' tomography. No need for projection data, it is included in the code.

Usage

gclc(phiorg, lamorg, phi, lam)

Arguments

Details

This may be defunct now.

Value

Note

Orignally from R. S. Crosson

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

lcgc

Examples

gclc(23, 35, 23.5, 35.6)

Description

vector of areas of polygons in map

Usage

geoarea(MAP, proj=NULL, ncut=10)

Arguments

Details

Uses sf function. If proj is NULL then the project is reset to UTM spherical for each element seperately to calculate the area in km. ncut is used to eliminate area calculations with strokes less than the specified number.

Value

vector of areas

Note

areas smaller than a certain tolerance are NA

Author(s)

Jonathan M. Lees<[email protected]>

See Also

sf::st_area

Description

Create and add Geological legend from GEOmap Structure

Usage

geoLEGEND(names, shades, zx, zy, nx, ny, side=1, cex=0.5)

Arguments

Details

Adds geological legend based on information provided. Legend is placed in margin.

Value

Graphical Side Effects

Note

If plot is resized, should re-run this as the units depend on the screen size information and the transformation of user coordinates.

Author(s)

Jonathan M. Lees<[email protected]>

Examples

## Not run: 

library(RPMG)
library(RSEIS)
library(GEOmap)
library(geomapdata)

data(cosogeol)
data(cosomap)
     data(faults)
     data(hiways)
     data(owens)

proj = cosomap$PROJ
   
XMCOL = setXMCOL()

newcol = XMCOL[cosogeol$STROKES$col+1]
cosocolnums = cosogeol$STROKES$col
cosogeol$STROKES$col = newcol
ss = strsplit(cosogeol$STROKES$nam, split="_")     

geo = unlist(sapply(ss  , "[[", 1))


UGEO = unique(geo)


mgeo = match( geo, UGEO )

gcol = paste(sep=".", geo, cosogeol$STROKES$col)


ucol = unique(gcol)

N = length(ucol)


spucol = strsplit(ucol,split="\.")     

       
names = unlist(sapply(spucol  , "[[", 1))

shades = unlist(sapply(spucol  , "[[", 2))

ORDN = order(names)
### example:


par(mai=c(0.5, 1.5, 0.5, 0.5) )

 plotGEOmapXY(cosomap, PROJ=proj,  add=FALSE, ann=FALSE, axes=FALSE)
  

     plotGEOmapXY(cosogeol, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  
geoLEGEND(names[ORDN], shades[ORDN], .28, .14, 4, 16, side=2)

####
par(mai=c(0.5, 0.5, 1.0, 0.5) )

 plotGEOmapXY(cosomap, PROJ=proj,  add=FALSE, ann=FALSE, axes=FALSE)
  

     plotGEOmapXY(cosogeol, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  
geoLEGEND(names[ORDN], shades[ORDN], .28, .14, 16, 6, side=3)


####
par(mai=c(0.5, 0.5, 0.5, 1) )

 plotGEOmapXY(cosomap, PROJ=proj,  add=FALSE, ann=FALSE, axes=FALSE)
  

     plotGEOmapXY(cosogeol, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  
geoLEGEND(names[ORDN], shades[ORDN], .28, .14, 3, 16, side=4)


####
par(mai=c(1.5, 0.5, 0.5, 0.5) )

 plotGEOmapXY(cosomap, PROJ=proj,  add=FALSE, ann=FALSE, axes=FALSE)
  

     plotGEOmapXY(cosogeol, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  
geoLEGEND(names[ORDN], shades[ORDN], .28, .14, 16, 3, side=1)



## End(Not run)

Description

Break a line at specified indices into a list

Usage

GEOmap.breakline(Z, ww)

Arguments

Value

Author(s)

Jonathan M. Lees<[email protected]>

Examples

Y=list()
Y$x=c(170,175,184,191,194,190,177,166,162,164)
Y$y=c(-54,-60,-60,-50,-26,8,34,37,10,-15)

GEOmap.breakline(Y, 5)

Description

Break up a polygon

Usage

GEOmap.breakpoly(Z, ww)

Arguments

Details

The NA values in Z represent breaks. GEOmap.breakpoly breaks the polygon up into individual strokes. The beginning and the ending of the stroke are combined.

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

fixCoastwrap, GEOmap.breakline

Examples

x=1:100
y = 1:100

ww = c(25, 53, 75)


A = list(x=x, y=y)

W = GEOmap.breakpoly(A , ww)

Description

Combine Two GEOmaps into one

Usage

GEOmap.cat(MAP1, MAP2)

Arguments

Details

Maps are combine consecutively.

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

GEOmap.Extract, GEOmap.CombineStrokes, list.GEOmap

Examples

data(coastmap)
CUBA = GEOmap.Extract(coastmap,90, INOUT="in" )

NSAMER =  GEOmap.Extract(coastmap,2, INOUT="in" )
AMAP = GEOmap.cat(CUBA, NSAMER)
plotGEOmap(AMAP )

Description

Combine strokes in a GEOmap list

Usage

GEOmap.CombineStrokes(MAP, SEL)

Arguments

Details

Stokes are combined in the order designated by the SEL index vector. The direction of the strokes is not modified - this may have to be fixed so that strokes align properly.

Value

GEOmap list

Author(s)

Jonathan M. Lees<[email protected]>

See Also

GEOmap.cat, GEOmap.Extract, GEOmap.CombineStrokes, list.GEOmap

Examples

data(coastmap)
SEL = which(coastmap$STROKES$nam=="Caribbean")


CAR =  GEOmap.Extract(coastmap, SEL, INOUT="in" )

plotGEOmap(CAR, MAPstyle=3, NUMB=TRUE)

CAR2 = GEOmap.CombineStrokes(CAR, SEL =c(6:15)  )

plotGEOmap(CAR2, MAPstyle=3, MAPcol='red' , add=TRUE)

Description

Extract or Exclude parts of a GEOmap list.

Usage

GEOmap.Extract(MAP, SEL, INOUT = "out")
fastExtract(MAP, SEL, INOUT = "out")
GEOmap.limit(MAP, LLlim )

Arguments

Details

Can either extract from the GEOmap data list with in, or exclude with out. fastExtract is the same but may be faster since it does not process all the strokes in the base GEOmap.

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

GEOmap.cat, GEOmap.Extract, GEOmap.CombineStrokes, list.GEOmap, getGEOmap, plotGEOmap, SELGEOmap, boundGEOmap,

Examples

data(coastmap)
SEL=which(coastmap$STROKES$nam=="AMERICAS")
NSAMER =  GEOmap.Extract(coastmap,SEL, INOUT="in" )
plotGEOmap(NSAMER)

Description

Inverse of list.GEOmap.

Usage

GEOmap.list(MAP, SEL = 1)

Arguments

Details

Returns the GEOmap strokes and instead of a long vector for the points they are broken down into a list of strokes.

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

GEOmap.cat, GEOmap.Extract, GEOmap.CombineStrokes, list.GEOmap

Examples

data(coastmap)
SEL=which(coastmap$STROKES$nam=='CUBA')
G = GEOmap.list(coastmap, SEL=SEL )

###  Lat-Lon of Cuba
G$LL

Description

Plot a set of Geological Symbols

Usage

GEOsymbols()

Details

Currently the choices in symbols are:

contact anticline syncline OverTurned-ant OverTurned-syn perp thrust normal dextral sinestral detachment bcars

Value

Graphical Side effect

Author(s)

Jonathan M. Lees<[email protected]>

See Also

bcars, thrust, teeth, SynAnticline, SSfault, horseshoe, strikeslip, OverTurned, normalfault, PointsAlong

Examples

GEOsymbols()

Description

Extract subset of a topographic database, interpolate and plot using the persp program.

Usage

GEOTOPO(TOPO, PLOC, PROJ, calcol=NULL, nx=500, ny=500, nb = 4, mb = 4, hb = 8, PLOT=TRUE)

Arguments

Details

The return matrix PMAT is a rotation matrix used for adding geographic (projected) data onto the perspective plot.

ETOPO5 or ETOPO2 can be downloaded from and installed using these links: http://leesj.sites.oasis.unc.edu/FETCH/GRAB/RPACKAGES/ETOPO2.RData and http://leesj.sites.oasis.unc.edu/FETCH/GRAB/RPACKAGES/ETOPO5.RData

Value

Note

If PLOT is false the transform matrix PMAT and the color mapping matrix Cmat will be returned as NA. To create these for future plotting, use TOPOCOL or LandSeaCol functions. TOPOCOL simply assigns values above sea level with one color scale and those below with under water colors. LandSeaCol requires a coastal map and fills in land areas with terrain colors and sea areas with blue palette colors.

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

See Also

subsetTOPO, TOPOCOL, LandSeaCol, settopocol, subsetTOPO, persp, DOTOPOMAPI

Examples

## Not run: 

library(geomapdata)

####  need to download and install ETOPO data
###   data(ETOPO5)
load(ETOPO5)
PLOC=list(LON=c(137.008, 141.000),LAT=c(34.000, 36.992),
            x=c(137.008, 141.000), y=c(34.000, 36.992) )

 PROJ = setPROJ(type=2, LAT0=mean(PLOC$y) , LON0=mean(PLOC$x) )
COLS = settopocol()
JMAT = GEOTOPO(ETOPO5, PLOC, PROJ, COLS$calcol, nx=1000, ny=1000, nb=8, mb=8, hb=12, PLOT=TRUE)

############  this plot can be duplicated by using the output or GEOTOPO


 PMAT = persp(JMAT$xo, JMAT$yo, JMAT$IZ$z, theta = 0, phi = 90, r=4000,
col=JMAT$Cmat[1:(JMAT$Dcol[1]-1), 1:(JMAT$Dcol[2]-1)] , scale = FALSE,
      ltheta = 120, lphi=60, shade = 0.75, border = NA, expand=0.001, box = FALSE )



## End(Not run)

Description

Extract from ETOPO5 or ETOPO2 data a rectangular subset of the full data.

Usage

getETOPO(topo, glat = c(-90, 90), glon = c(0, 360))

Arguments

Details

ETOPO2 and ETOPO5 are stored in a strange way: the lons are okay the latitudes are upside down.

ETOPO5 or ETOPO2 can be downloaded from and installed using these links: http://leesj.sites.oasis.unc.edu/FETCH/GRAB/RPACKAGES/ETOPO2.RData and http://leesj.sites.oasis.unc.edu/FETCH/GRAB/RPACKAGES/ETOPO5.RData

Value

Returns a matrix with attributes in lat-lon that are correct for usage in image or other R imaging programs.

Author(s)

Jonathan M. Lees<[email protected]>

See Also

image

Examples

## Not run: 
library(geomapdata)
###  Download and install ETOPO Data
##    data(ETOPO5)
load(ETOPO5)
glat =c(45.4, 49)
glon = c(235, 243) 
b5 = getETOPO(ETOPO5, glat, glon)
image(x=attr(b5, 'lon'), y=attr(b5,'lat'), z=b5, col=terrain.colors(100) )
contour(   x=attr(b5, 'lon'), y=attr(b5,'lat'), z=b5, add=TRUE)

## End(Not run)

Description

Get Geomap from ascii files

Usage

getGEOmap(fn)

Arguments

Details

Files are stored as a pair: rootname.strks and rootname.pnts

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

plotGEOmapXY, boundGEOmap

Examples

## Not run: 
library(geomapdata)

 data(cosomap)
     data(faults)
     data(hiways)
     data(owens)

cosogeol = getGEOmap("/home/lees/XMdemo/GEOTHERM/cosogeol")

cosogeol = boundGEOmap(cosogeol)


 proj = cosomap$PROJ
    
plotGEOmapXY(cosomap, PROJ=proj,  add=FALSE, ann=FALSE, axes=FALSE)
  
 plotGEOmapXY(cosogeol, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  

  plotGEOmapXY(cosomap, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  
  plotGEOmapXY(faults, PROJ=proj,  add=TRUE, ann=FALSE, axes=FALSE)
  

## End(Not run)

Description

Get rectangular perimeter of region defined by set of Lat-Lon

Usage

getGEOperim(lon, lat, PROJ, N)

Arguments

Details

perimeter is used for antipolygon

Value

List:

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

Examples

###  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) )



perim= getGEOperim(PLOC$LON, PLOC$LAT, PROJ, 50)

Description

Get points along great circle between two locations

Usage

getgreatarc(lat1, lon1, lat2, lon2, num)

Arguments

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

getgreatarc, distaz

Examples

PARIS = c(48.8666666666667, 2.33333333333333)
RIODEJANEIRO =c( -22.9,  -43.2333333333333)

g = getgreatarc(PARIS[1],PARIS[2], RIODEJANEIRO[1], RIODEJANEIRO[2],
100)
library(geomapdata)
data(worldmap)

plotGEOmap(worldmap, add=FALSE, shiftlon=180)

lines(g$lon+180, g$lat)

Description

Estimate a size for plotting earthqukes recorded as a logarithmic scale

Usage

getmagsize(mag, minsize = 1, slope = 1, minmag = 0, maxmag = 8, style = 1)

Arguments

Details

The idea is to have a scale reflect the size of the earthquake. The default style (1) has a few parameters left over from old program geotouch.

Value

vector of sizes for plotting

Author(s)

Jonathan M. Lees<[email protected]>

Examples

mag = 0:9

x = runif(10, 1, 100)
y  = runif(10, 1, 100)

g = getmagsize(mag)

plot(c(0, 100), c(0, 100), asp=1, type='n')

points(x, y, pch=1, cex=g)

Description

Given a set of lat lon pairs, return a new set of tic marks

Usage

getnicetix(lats, lons)

Arguments

Value

Author(s)

Jonathan M. Lees<[email protected]>

See Also

niceLLtix

Examples

proj = setPROJ(7, LAT0 = 0 , LON0= -93)
rx = c(652713.4, 656017.4)
ry = c(1629271, 1631755)

   gloc = XY.GLOB(rx, ry, proj)
    
    G = getnicetix(gloc$lat, gloc$lon)
 
print(G)

Description

Get a spline curve along a set of points

Usage

getspline(x, y, kdiv)

Arguments

Value

LIST:

Author(s)

Jonathan M. Lees<[email protected]>

Examples

plot(c(-5,5), c(-5,5), asp=1, type='n' )
ff=list()
ff$x=c(-4.850,-4.700,-3.934,-2.528, 0.603, 2.647, 3.861, 2.626)
ff$y=c(-4.045,-2.087,-0.710, 0.172, 1.291, 2.087,-0.753,-4.131)

lines(ff, col='red')
G =getspline(ff$x, ff$y, kdiv=20)

lines(G, col='blue')

Description

Get a spline curve along a set of points

Usage

getsplineG(x, y, kdiv)

Arguments

Value

LIST:

Author(s)

Jonathan M. Lees<[email protected]>

Examples

plot(c(-5,5), c(-5,5), asp=1, type='n' )
ff=list()
ff$x=c(-4.850,-4.700,-3.934,-2.528, 0.603, 2.647, 3.861, 2.626)
ff$y=c(-4.045,-2.087,-0.710, 0.172, 1.291, 2.087,-0.753,-4.131)

lines(ff, col='red')
G =getsplineG(ff$x, ff$y, kdiv=20)

lines(G, col='blue')

Description

Example of how to use RPMG button functions. This example shows how to plot a DEM and interactively change the plot and find projected cross-sections through a surface.

Usage

GETXprofile(jx, jy, jz, LAB = "A", myloc = NULL, PLOT = FALSE, NEWDEV=TRUE,  asp=1)

Arguments

Details

The program uses a similar input format as image or contour, with structure from the locator() function of x and y coordinates that determine where the cross section is to be extracted.

Value

Returns a list of x,z values representing the projected values along the cross section.

Note

The program is an auxiliary program provided to illustrate the RPMG interactive R analysis.

Author(s)

Jonathan M. Lees<[email protected]>

See Also

locator, image

Examples

## Not run:  
#######  get data 
   data(volcano)
####  extract dimensions of image
   nx = dim(volcano)[1]
   ny = dim(volcano)[2]

###  establish units of image
   jx = 10*seq(from=0, to=nx-1)
   jy = 10*seq(from=0, to=ny-1)

####  set a letter for the cross section
   LAB = LETTERS[1]

###  coordinates of cross section on image
###  this is normally set by using the locator() function
   x1 = 76.47351
   y1 = 231.89055
   x2 = 739.99746
   y2 = 464.08185

## extract and plot cross section

 GETXprofile(jx, jy, volcano, myloc=list(x=c(x1, x2), y=c(y1, y2)), LAB=LAB, PLOT=TRUE)
            
## End(Not run)

Description

Convert from GLOBAL LAT-LON to X-Y

Usage

GLOB.XY(LAT, LON, PROJ.DATA)

Arguments

Details

Units should be given according to the projection. This is the inverse of XY.GLOB.

Value

Author(s)

Jonathan M. Lees<jonathan.lees.edu>

References

Snyder, John P., Map Projections- a working manual, USGS, Professional Paper, 1987.

See Also

XY.GLOB

Examples

proj = setPROJ(type = 2, LAT0 =23, LON0 = 35)

### get lat-lon
LL = XY.GLOB(200, 300, proj)


##  find x-y again, should be the same
XY = GLOB.XY(LL$lat, LL$lon, proj)
XY

Description

Plot globe with orthogonal

Usage

GLOBE.ORTH(lam0, phi1, R = 1, plotmap = TRUE, plotline=TRUE, add=FALSE,
 map = coastmap, mapcol = grey(0.2), linecol = grey(0.7), fill=FALSE)

Arguments

Details

Plots whole globe with grid.

Value

Graphical Side effects

Author(s)

Jonathan M. Lees<[email protected]>

References

Snyder, John P., Map Projections- a working manual, USGS, Professional Paper, 1987.

See Also

setPROJ, projtype, plotGEOmap

Examples

######  simple map of world viewed at 40 degrees latitude
R = 1
R.MAPK = 6378.2064

phi1=40

viewlam = seq(from=0, to=340, by=2)

data(coastmap)

K=1
GLOBE.ORTH(viewlam[K], phi1, R=1, plotmap=TRUE)
######


OLIM = c(20, 40, 10, 40)
TLIM = c(-20, -10, -30, -10)


 olat = runif(1, OLIM[1], OLIM[2])
         olon = runif(1, OLIM[3], OLIM[4] )

          tlat = runif(1,TLIM[1], TLIM[2] )
         tlon = runif(1, TLIM[3], TLIM[4])

GLOBE.ORTH(olon, olat, 1,plotmap=FALSE )

 XYorg = ortho.proj(olat, olon, olon, olat, 1)
 XYtarg = ortho.proj(tlat, tlon, olon, olat, 1)

points( XYorg , col='red')
points(XYtarg , col='blue')
 da = distaz(olat, olon, tlat, tlon)
 ed2 =  Ellipsoidal.Distance(olat, olon, tlat, tlon, a=R.MAPK*1000, b=R.MAPK*1000)

  A = along.great(olat*pi/180, olon*pi/180,
seq(from=0, to=da$del, by=2)*pi/180,  da$az*pi/180)

     lat=A$phi*180/pi
     lon = A$lam*180/pi

 XYalong = ortho.proj(lat, lon, olon, olat, 1)

lines(XYalong , col='purple')

M = merid(tlon, lat1=tlat, phi1=olat, lam0=olon, R=1, by=2)

lines(M$x, M$y, col='blue' )

M2 = merid(olon, lat1=olat,  phi1=olat, lam0=olon,R=1, by=2)

lines(M2$x, M2$y, col='red' )

 leg = c( paste("del=", round(da$del)), paste("DA=", round(da$az),
 round(da$baz) ),
 paste("ED=", round(ed2$az) ,  round(ed2$revaz) ))
 
legend("topleft", legend=leg)

Description

Plot global view of the earth

Usage

GlobeView(phicen, lamcen, worldmap, MAXR, SEL = 1,
circol = rgb(1, 0.8, 0.8), innercol = "white", linecol = rgb(0, 0, 0),
 mapcol = rgb(0, 0, 0), backcol = "white", add=FALSE, antip=TRUE)

Arguments