Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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!

like image 445
Caroline4 Avatar asked May 03 '15 19:05

Caroline4


People also ask

How do you convert XY coordinates to latitude and longitude?

Calculate latitude and longitude using the formula: latitude = asin (z/R) and longitude = atan2 (y,x). In this formula, we have the values of x, y, z and R from step 2. Asin is arc sin, which is a mathematical function, and atan2 is a variation of the arc tangent function. The symbol * stands for multiplication.

Is UTM the same as lat long?

They are just two different ways of positioning a point. Many experienced users prefer UTM over latitude/longitude when using 7.5' topographic quadrangle maps. Ocean-going sailors and other marine users almost always use latitude/longitude because navigation charts are optimized for this method.

How do I express my UTM coordinates?

For locations north of the equator, give the zone number, the letter “N”, and the distance from the central meridian in meters. For example, Port-au-Prince, Haiti, which is located at lat/long 18.5425*, - 72.3386*, has UTM coordinates of 18N 780950E 2052283N.

Why is UTM better than lat long?

UTM Provides a constant distance relationship anywhere on the map. In angular coordinate systems like latitude and longitude, the distance covered by a degree of longitude differs as you move towards the poles and only equals the distance covered by a degree of latitude at the equator.


1 Answers

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)

 library(terra)
 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(rgdal)
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.

like image 100
Robert Hijmans Avatar answered Oct 07 '22 05:10

Robert Hijmans