Caroline4
Caroline4

Reputation: 157

How to convert UTM coordinates to lat and long in R

I am trying to run a species distribution model and need to create background points to run my logistic regression model. I have just created 500 randomPoints but they are in UTM coordinates and I need lat and long. Is there a way to convert them to lat and long in R? If so, can you share the code with me? I am fairly new to R. Thanks!

Upvotes: 13

Views: 36748

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47091

If you need long/lat you should probably generate the random points using that coordinate reference system. But otherwise you can do something like the below.

I first generate an example set of coordinates (points)

 set.seed(1)
 x <- runif(10, -10000, 10000)   
 y <- runif(10, -10000, 10000)   
 points <- cbind(x, y)

Now, assuming you know the coordinate reference system (CRS) of your points, you can create a spatial object and assign that known CRS. In my example, the points are in the UTM, zone 10, datum=WGS84 projection.

library(terra)
v <- vect(points, crs="+proj=utm +zone=10 +datum=WGS84  +units=m")
v
# class       : SpatVector 
# geometry    : points 
# dimensions  : 10, 0  (geometries, attributes)
# extent      : -8764.275, 8893.505, -6468.865, 9838.122  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 

Now we can transform these to another CRS, for example to longitude/latitude

y <- project(v, "+proj=longlat +datum=WGS84")
y
# class       : SpatVector 
# geometry    : points 
# dimensions  : 10, 0  (geometries, attributes)
# extent      : -127.5673, -127.4091, -0.05834327, 0.08873723  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
 

And you can extract the coordinates like this

lonlat <- geom(y)[, c("x", "y")]
head(lonlat, 3)
#             x             y
#[1,] -127.5308 -0.0530354276
#[2,] -127.5117 -0.0583432750
#[3,] -127.4757  0.0337371933

You can of course also do the inverse

back <- project(y, "+proj=utm +zone=10 +datum=WGS84 +units=m")

The same thing can be done with the sf package, or with the old sp package. With sp, create a SpatialPoints object and use spTransform.

library(sp)
sputm <- SpatialPoints(points, proj4string=CRS("+proj=utm +zone=10 +datum=WGS84")) 
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
lnlt <- coordinates(spgeo)
 

I used UTM zone 10 in the example. But note that there are 60 UTM zones and you have to pick one. Each covers a strip of (360/60=) 6 degrees. You should not use UTM if your data spans a lot of longitude or crosses UTM zones. For a longitude between [-180, 180) You can compute the zone you need like this

utm_zone <- function(longitude) { trunc((180 + longitude) / 6) + 1 }

longs <- c(-122,-119, -118)
utm_zone(min(longs))
# [1] 10
utm_zone(max(longs))
# [1] 11

utm_zone(max(longs))

Or have a look at a map like this one

A point of confusion with UTM in the use of "false northings" for locations in the Southern hemisphere, to avoid negative coordinates. This is done by adding 10,000,000 to the y coordinates, as illustrated below using the +south element.

s <- vect(cbind(174, -44), crs="+proj=longlat +datum=WGS84")
geom(project(s, "+proj=utm +zone=59"))[, c("x", "y")]
#       x          y 
#740526.3 -4876249.1 

geom(project(s, "+proj=utm +south +zone=59"))[, c("x", "y")]
#       x         y 
#740526.3 5123750.9 

Also note that I use the "PROJ4" notation to define the CRS. That works fine if the datum used (a datum is a model of the shape of the earth's surface) is WGS84 or NAD83. If that is not the case you would need to use an "EPSG" code or a "WKT2" description of your CRS.

Upvotes: 31

Related Questions