[R-sig-Geo] Plotting an image on a stereographic grid orientation

[White.Denis at epamail.epa.gov]

Here is a script using two functions.  I hope this helps.  Note that
this version of the stereographic projection uses a spherical model of
the earth not an ellipsoidal, and is set up for the north polar aspect
of the projection, since that appears to be most appropriate for the
kind of data you have in your example.  (I also included a function for
the equatorial aspect of the projection in case you have use for that.
Both the north polar and the equatorial are special cases of the more
general oblique aspect which is what mapproject and most projection
functions would use.)

best wishes,
Denis White

-------------------------
stereo <- function (lon, lat) # north polar aspect
{
    R <- 6371 # ~ GRS80 Authalic sphere
    v <- tan (pi / 4 - lat * pi / 180 / 2)
    x <- 2 * R * v * sin (lon * pi / 180)
    y <- -2 * R * v * cos (lon * pi / 180)
    x[is.na (lat)] <- lon[is.na (lat)]
    xy <- cbind (x=x, y=y)
    xy
}

stereo.equator <- function (lon, lat) # equatorial aspect
{
    R <- 6371 # ~ GRS80 Authalic sphere
    coslat <- cos (lat * pi / 180)
    k <- 2 / (1 + coslat * cos (lon * pi / 180))
    x <- R * k * coslat * sin (lon * pi / 180)
    y <- R * k * sin (lat * pi / 180)
    x[is.na (lat)] <- lon[is.na (lat)]
    xy <- cbind (x=x, y=y)
    xy
}

graticule <- function (east=180, west=-180, north=90, south=-90,
    ginc=10, dinc=1)
{
    ew <- east - west
    ns <- north - south
    nr <- (ew / dinc + 1) * (ns / ginc + 1) +
        (ns / dinc + 1) * (ew / ginc + 1) +
        ew / ginc + ns / ginc + 2
    geo <- matrix (0, nrow=nr, ncol=2)
    i <- n <- 0

    for (lon in seq (west, east, by=ginc))
    {
        i <- i + 1
        n <- n + 1
        geo[i, 1] <- n
        geo[i, 2] <- NA
        for (lat in seq (south, north, by=dinc))
        {
            i <- i + 1
            geo[i, 1] <- lon
            geo[i, 2] <- lat
        }
    }
    for (lat in seq (south, north, by=ginc))
    {
        i <- i + 1
        n <- n + 1
        geo[i, 1] <- n
        geo[i, 2] <- NA
        for (lon in seq (west, east, by=dinc))
        {
            i <- i + 1
            geo[i, 1] <- lon
            geo[i, 2] <- lat
        }
    }
    geo
}

east <- 180
west <- -180
north <- 90
south <- 50
ginc <- 2.5
dinc <- 0.5

grat <- graticule (east, west, north, south, ginc, dinc)
grat <- stereo (grat[,1], grat[,2])

lon <- seq (west, east, ginc)
lat <- seq (south, north, ginc)
pts <- expand.grid (lon, lat)
pts <- stereo (pts[,1], pts[,2])
val <- rnorm (nrow (pts), sd=10)
val <- cut (val, breaks=7, include.lowest=TRUE, labels=FALSE)
kol <- rainbow (7)

plot.map (grat[,1], grat[,2])
lines (grat)
points (pts, col=kol[val], pch=20)
-------------------------

r-sig-geo-bounces at stat.math.ethz.ch wrote on 2005-08-12 00:45:14:

> Dear all
>
> First of all, I am sorry for all inconvinience I am causing you, but I

> have a problem on plotting some data on stereographic grids.
>
> Say I have data given to be (for simplicity):
> lon<-seq(-180,180,2.5)
> lat<-seq(50,90,2.5)
> data<-matrix(rnorm(145*17,sd=10),ncol=17,nrow=145)
>
> And this I want to plot on a map like this (maps and mapproj package
> needed):
> map("world",proj="stereographic",xlim=c(-180,180),ylim=c(50,90),
> interior=FALSE,lwd=1)
>
map.grid(lim=c(-180,180,50,90),nx=5,ny=3,pretty=TRUE,col=1,label=FALSE)
>
> I would like to have a colour presentation of the data. Before I have
> used the image command to plot a datamatrix when there is no specific
> projection, and I would like to use the image command on this problem
too.
>
> I hope someone have time to help me.
>
> Best Wishes
> Dag Johan Steinskog
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo