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

[Roger Bivand]

On Fri, 12 Aug 2005, Barry Rowlingson wrote:

> Dag Johan Steinskog wrote:
> 
> > 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.
> 
>   I think you'd have to loop over every cell in your matrix, get the 
> projected coordinates of the cell and intermediate points along the 
> edges, then draw a filled polygon in a colour relating to the value you 
> want to display.
> 
>   The code in map.grid might help, since it shows how to compute 
> projected coordinates, but the 'grid' is only a set of latitude and 
> longitude lines. You'd need to work out coordinates for each cell in 
> order to fill it in.
> 
>   I could probably bash out a function to do it in half a day, sadly I 
> dont have half a day free at the moment!
> 

I'd got as far as:

This is about warping a raster (note that 180 and -180 in lon overlap) to
the polar stereographic projection, or possibly faking a warp by
projecting the raster cells as polygons (this doesn't give the two curved
edges, but avoids resampling the data).

A further issue is that map() from the maps package imposes its own
display coordinate system.

And:

> map("world", proj="stereographic", xlim=c(-180,180), ylim=c(50,90), 
+ interior=FALSE, lwd=1)
> polygon(mapproject(c(0,0,2.5,2.5,0), c(87.5,90,90,87.5,87.5), 
+ proj="stereographic"), col="red", border=NA)
> polygon(mapproject(c(0,0,2.5,2.5,0), c(57.5,60,60,57.5,57.5), 
+ proj="stereographic"), col="blue", border=NA)

does work, but the polygons would need plotting first, so as not to 
overwrite the coastlines. Curiously, 

> mapproject(c(60,60,62.5,62.5,60), c(57.5,60,60,57.5,57.5),proj="stereographic")

gives the same coordinates as x=c(0,0,2.5,2.5,0), so maybe it isn't so 
easy.

as("GridTopology", "SpatialRings"), anyone? 

Dag, when do you need this (apart from yesterday)?

This from project() in spproj (not yet on CRAN):

> project(cbind(c(60,60,62.5,62.5,60), c(57.5,60,60,57.5,57.5)), 
"+proj=stere +lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84")
 +proj=stere +lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84 
+ellps=WGS84 +towgs84=0,0,0
        [,1]     [,2]
[1,] 3227409 -1863346
[2,] 2967384 -1713220
[3,] 3039289 -1582154
[4,] 3305616 -1720795
[5,] 3227409 -1863346
> project(cbind(c(-60,-60,-62.5,-62.5,-60), c(57.5,60,60,57.5,57.5)), 
"+proj=stere +lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84")
 +proj=stere +lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84 
+ellps=WGS84 +towgs84=0,0,0
         [,1]     [,2]
[1,] -3227409 -1863346
[2,] -2967384 -1713220
[3,] -3039289 -1582154
[4,] -3305616 -1720795
[5,] -3227409 -1863346
> project(cbind(c(0,0,2.5,2.5,0), c(57.5,60,60,57.5,57.5)), "+proj=stere 
+lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84")
 +proj=stere +lat_ts=-90 +lat_0=90 +lon_0=0 +k_0=1.0 +datum=WGS84 
+ellps=WGS84 +towgs84=0,0,0
         [,1]     [,2]
[1,]      0.0 -3726691
[2,]      0.0 -3426439
[3,] 149459.2 -3423178
[4,] 162556.0 -3723144
[5,]      0.0 -3726691

looks promising (project() is based on Barry's code).

Roger



> Barry
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no