[R-sig-Geo] LL to UTM

Gabor Grothendieck ggrothendieck at gmail.com
Tue Mar 6 15:37:04 CET 2007


Thanks for the clarification.

On 3/6/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
>




More information about the R-sig-Geo mailing list