[R-sig-Geo] plotting latlon axes of sp object of different projection

[Roger Bivand]

On Tue, 5 Jul 2011, Horacio Samaniego wrote:

> Thanks a lot Roger,
>
> This is a good recipe, I'ld never have figured by my own!
>
> To fine tune this, I need to place the labels on the other sides
> (north, east) of the graph. Modifying the pos column of the  grdat
> object (output of gridat()) does not do the trick. I imagine that any
> modification should happen in the coordinates of that object.
>
> So I wonder if that is achievable earlier in the process.
>
> Any helpfull thought on doing that? It would be much appreciate.

I've added an additional argument to gridat(), and committed to SVN on the 
R-forge rspatial project, module sp. Both the coordinates and the pos need 
to be changed, pos to 3 for 1 and 4 for 2. The coordinates are more 
difficult, because they must be offset by the lengths of the grid lines 
before projection. If you can install from source, please try out the 
R-forge code, otherwise wait for a new release, or try to work things out 
from for example:

sapply(slot(slot(grd, "lines")[[1]], "Lines"), LineLength)
sapply(slot(slot(grd, "lines")[[2]], "Lines"), LineLength)

and work out how to shift the coordinates, before spTransform() to 
projected coordinates.

Hope this helps,

Roger

PS. Edzer has added a function to the next rgdal release to help in 
overplotting long-lat grids.


>
> thanks!
>
> H
>
> On Tue, Jul 5, 2011 at 5:20 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> On Mon, 4 Jul 2011, Horacio Samaniego wrote:
>>
>>> I need to add a lat/lon axes to an equal area plot of a
>>> SpatialLinesDataFrame object currently in units of meters. Is that
>>> achievable using sp/maptools?
>>>
>>> The .prj file where the object comes from is:
>>> PROJCS["Lambert_Azimuthal_Equal_Area",GEOGCS["GCS_unnamed
>>>
>>> ellipse",DATUM["D_unknown",SPHEROID["Unknown",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
>>>
>>> I'ld hate to install arcinfo to generate 1 plot! (plus the plot is
>>> already nice as it is!!! just need the axes!)
>>
>> No, the axes are in the same coordinate system as the plot, and would be
>> curves here. You can overplot a long-lat grid on a projected plot:
>>
>> library(sp)
>> data(meuse)
>> coordinates(meuse) <- ~x+y
>> proj4string(meuse) <- CRS("+init=epsg:28992")
>> plot(meuse)
>> library(rgdal)
>> meuse_ll <- spTransform(meuse, CRS("+proj=longlat +datum=WGS84"))
>> grd <- gridlines(meuse_ll)
>> grd_x <- spTransform(grd, CRS("+init=epsg:28992"))
>> plot(grd_x, add=TRUE, lty=2)
>> grdat_ll <- gridat(meuse_ll)
>> grdat_x <- spTransform(grdat_ll, CRS("+init=epsg:28992"))
>> text(coordinates(grdat_x), labels=parse(text=grdat_x$labels),
>>  pos=grdat_x$pos, offset=grdat_x$offset)
>>
>> Hope this helps,
>>
>> Roger
>>
>>>
>>> thanks,
>>> I'll summarize!
>>>
>>> H
>>>
>>>
>>
>> --
>> Roger Bivand
>> Department of Economics, NHH Norwegian School of Economics,
>> 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
>>
>>
>
>
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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