[R-sig-Geo] Intersection of two coordinates (triangulation)

[Barry Rowlingson]

On Fri, Sep 11, 2009 at 7:43 PM, Breitbach, Nils <breitbach at uni-mainz.de> wrote:
> Hi Community,
>
> at the moment I am starting to analyse my telemetry data trying to be independent from commercial software. So far I worked with ESRI's ArcInfo to project data from one system into another. No I want to do this with R and so I am looking for a function that is able to do it.
>
> My dataset consists of coordinates in longitude and latitude using WGS1984 and decimal degrees. Now I want to project these coordinates a) into UTM zone 32 north coordinates and b) into Gauss Krueger zone 3 coordinates. Is it possible using R functions or do I need to use GRASS GIS for example? If it is possible, how? Since with the pevious software used I did this type of operations by mouse clicks and ticking options, I would probably need a somewhat detailed explanation. ;)

 How far have you got at the moment? What form is your data in?
Shapefiles? CSV files?

 The way to do coordinate transforms in R is with the spTransform
function from the rgdal package. Here I'll make some fake data and
show you how it's done:

 map=data.frame(x=runif(100,0,5),y=runif(100,20,22),s=1:100)
 plot(map$x,map$y)

this is just a plain R data frame. We need to make it spatial, so we
give it coordinates:

 coordinates(map)=~x+y
 plot(map)

 now we see a different view - the aspect ratio is now set to 1.
Currently the points don't have a spatial reference. We'll give them
one, by setting its 'proj4string':

  proj4string(map)=CRS("+init=epsg:4326")

 PROJ4 is a library for doing map transformations, and the sp package
uses it. The "epsg:4326" bit is the epsg code for WGS84 (see
www.epsg.org for more info on epsg codes).

 Now do:

 plot(map);axis(1);axis(2)

 to see the coordinates.

So how do we transform? spTransform to a new CRS!

Look up your EPSG code either from www.epsg.org or maybe here:

http://spatialreference.org/ref/epsg/32632/

 and transform....

 map2 = spTransform(map,CRS("+init=epsg:32632"))
 plot(map2);axis(1);axis(2)

 (Note I'm not sure if the numbers I chose are valid in that zone).

Gauss-Krueger zone 3 seems to be 2166:

http://spatialreference.org/ref/epsg/2166/

Hope this helps!

Barry