[R-sig-Geo] LL to UTM

[Roger Bivand]

On Tue, 6 Mar 2007, Gabor Grothendieck wrote:

> Great!  Your LL2UTM function was exactly what I wanted.
> 
> My data sets are each about a city block in size and I mostly work
> in UTM on this project but I do get a few points in lat long. Its small
> enough that I have been using that web site I mentioned but it gets
> tedious after a while and this should be a big help.
> 
> I did notice that the R function below is consistent with the web site
> output only to 5 digits after the decimal place in x and 2 digits after
> the decimal place in y.  I don't know how to determine whether R or
> the web site is the more accurate but at any rate either seem likely
> to be well within the accuracy I need since, in truth, I don't need
> high precision.

The units are metres, so centimetre accuracy in y (measured from the 
Equator) and much better in x isn't bad! The web calculator may have a 
different number of digits of precision in its ellipsoid definitions.

Roger


> 
> Thanks.
> 
> 
> 
> On 3/6/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > On Tue, 6 Mar 2007, Gabor Grothendieck wrote:
> >
> > > Thanks. I looked at spTransform prior to posting but don't understand it.
> > > What I want is to create a function that returns UTM coordinates:
> > >
> > > LL2UTM <- function(lat, long, zone = 18) {
> > > ...
> > > }
> >
> > library(rgdal)
> > LL2UTM <- function(lat, long, zone=18) {
> >  project(cbind(long, lat), paste("+proj=utm +zone=", zone, sep=""))
> > }
> >
> > gets you there, for example for:
> >
> > lat <- seq(30,45,1)
> > long <- seq(-78,-72,1)
> > grd <- expand.grid(lat, long)
> > LL2UTM(grd[,1], grd[,2])
> >
> > The project() interface to PROJ4 does not support datum transfromation,
> > the spTransform() method does. The cost is having to convert the data to
> > a Spatial* object:
> >
> > SP_grd <- SpatialPoints(cbind(grd[,2], grd[,1]),
> >  proj4string=CRS("+proj=longlat +datum=NAD83"))
> > SP_grd1 <- spTransform(SP_grd, CRS("+proj=utm +zone=18 +datum=NAD83"))
> >
> > returning a SpatialPoints object - use the coordinates() method to
> > retrieve as a matrix.
> >
> > Assuming that your input geographical coordinates are in NAD84/WGS84, and
> > the output UTM coordinates are in the same datum, project() will be OK. If
> > you need datum transformation, which may lead to errors of hundreds of
> > metres if not used, the PROJ4 string needs more detail. If this is at the
> > block scale, datum transformation will make a difference if the input and
> > output specifications vary (say placing points on a UTM map not in
> > WGS84/NAD83). At the continental scale the differences are typically not
> > great, but then you wouldn't use UTM anyway.
> >
> > Hope this helps,
> >
> > Roger
> >
> > >
> > > If necessary, zone=18 can be hardcoded in the function and that
> > > arg removed.
> > >
> > > I assume that using spTransform its just a one line body. Can
> > > you give me the specific line that it should be?
> > >
> > > Thanks.
> > >
> > >
> > > On 3/6/07, Edzer J. Pebesma <e.pebesma at geo.uu.nl> wrote:
> > > > Gabor,
> > > >
> > > > package rgdal provides an interface to the PROJ.4 library for projection
> > > > of geographical data. Look for the function (or rather method)
> > > > spTransform. It takes any of the spatial classes provided by package sp,
> > > > and afaik any of the known projection systems, including UTM and
> > > > ellipsoids (LL reference "model").
> > > > --
> > > > Edzer
> > > >
> > > > Gabor Grothendieck wrote:
> > > > > I am currently using this web page to convert LL to UTM (and UTM to LL).
> > > > >
> > > > >    http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html
> > > > >
> > > > > All my points are in zone 18.  How would I accomplish the same
> > > > > thing using R?
> > > > >
> > > > > Thanks.
> > > > >
> > > > > _______________________________________________
> > > > > R-sig-Geo mailing list
> > > > > R-sig-Geo at stat.math.ethz.ch
> > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > >
> > > >
> > > >
> > >
> > > _______________________________________________
> > > 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
> >
> >
> 

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