## Monday, May 30, 2011

### map.xyz(): interpolation of XYZ data and projection onto a map

I am still struggling to get a grasp of R's mapping capabilities. Part of my frustration lies in the fact that I often work on areas near the poles, which complicates interpolation across the 180 degree line. For smaller areas, interpolation can be done using the interp() function in the package akima. I have taken the results from interp and projected the image onto a map. You will need the akima, maps, and mapproj packages and the functions new.lon.lat(), earth.dist(), and pos2coord().

As an example I have mapped the distance from Mecca, Saudi Arabia:

The function...

### Array position to matrix coordinates conversion

```#A function that is sometimes useful in determining the
#coordinate(i.e. row and column number) of a matrix position
#(and vice-versa).
#Either a vector of positions ("pos")
#OR a 2 column matrix of matrix coordinates, ("coord", i.e. cbind(row,col)),
#AND the matrix dimentions must be supplied (dim.mat, i.e. c(nrow,ncol)).
pos2coord<-function(pos=NULL, coord=NULL, dim.mat=NULL){
if(is.null(pos) & is.null(coord) | is.null(dim.mat)){
stop("must supply either 'pos' or 'coord', and 'dim.mat'")
}
if(is.null(pos) & !is.null(coord) & !is.null(dim.mat)){
pos <- ((coord[,2]-1)*dim.mat)+coord[,1]
return(pos)
}
if(!is.null(pos) & is.null(coord) & !is.null(dim.mat)){
coord <- matrix(NA, nrow=length(pos), ncol=2)
coord[,1] <- ((pos-1) %% dim.mat) +1
coord[,2] <- ((pos-1) %/% dim.mat) +1
return(coord)
}
}```
Created by Pretty R at inside-R.org

## Sunday, May 29, 2011

### R functions for Earth geographic coordinate calculations

Here are some functions that I regularly use for geographic data (e.g. binning, filtering, calculation of new positions etc.).

```#distance in kilometers between two long/lat positions (from "fossil" package)
earth.dist <- function (long1, lat1, long2, lat2)
{
dlon <- b2 - a2
dlat <- b1 - a1
a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
c <- 2 * atan2(sqrt(a), sqrt(1 - a))
R <- 6378.145
d <- R * c
return(d)
}```
Created by Pretty R at inside-R.org

```#degree bearing between two long/lat positions (from "fossil" package)
earth.bear <- function (long1, lat1, long2, lat2)
{
dlon <- b2 - a2
bear <- atan2(sin(dlon) * cos(b1), cos(a1) * sin(b1) - sin(a1) *
cos(b1) * cos(dlon))
deg <- (bear%%(2 * pi)) * (180/pi)
return(deg)
}```
Created by Pretty R at inside-R.org

```new.lon.lat <-
function (lon, lat, bearing, distance)
{
d <- distance/6378.145
nlat <- asin(sin(a1) * cos(d) + cos(a1) * sin(d) * cos(tc))
dlon <- atan2(sin(tc) * sin(d) * cos(a1), cos(d) - sin(a1) *
sin(nlat))
nlon <- ((a2 + dlon + pi)%%(2 * pi)) - pi
return(npts)
}```
Created by Pretty R at inside-R.org

```#tells which lon lat positions are within the defined limits to the west, east, north, and south
lon.lat.filter <-
function (lon_vector, lat_vector, west, east, north, south)
{
if(west>east) {
lon_vector_new=replace(lon_vector, which(lon_vector<0), lon_vector[which(lon_vector<0)]+360)
east_new=east+360
} else {
lon_vector_new=lon_vector
east_new=east
}
hits=which(lon_vector_new < east_new & lon_vector_new > west & lat_vector < north & lat_vector > south)
return(hits)
} ```
Created by Pretty R at inside-R.org