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

[Dag Johan Steinskog]

Dear Roger and Barry (and rest of you)

I need this as soon as possible I think. I am working with climate 
analysis, and this need is quite urgent in general, not only for me. I 
will do some testing on the things you have written down, but I feel 
this is something that should be done.

All the best
Dag Johan

Roger Bivand skrev:

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