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!
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.
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.
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.
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.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With