Reputation: 157
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
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