CODEWITHSUNDEEP
rgeomap-projections
Mark Heckmann
Mark Heckmann

Reputation: 11431

reverse map-projection: how to get lat/lon coordinates from projected coordinates

I have a set of lat/lon coords which I can project using, for example, Mollweide projection.

library(mapproj)
set.seed(0)
n <- 100
s <- data.frame(lon = rnorm(n, 0, 60),
                lat = rnorm(n, 0, 40))
p <- mapproject(s$lon, s$lat, proj="mollweide", par=NULL, 
                orientation=c(90,200,0))                 

# plot projected coors
plot(p$x, p$y, type="n", asp=1/1, bty="n")
map.grid(c(-180, 180, -90, 90), nx=20, ny=20, 
         font=1, col=grey(.6), labels=F)
points(p$x, p$y, pch="x", cex = .8)

# a point to reverse project
points(1,0, pch=16, col="red", cex=2)

enter image description here

Now, I have a scenario where I need to do some calculations on the projected coordinates and reverse project the results back to lat/lon coords. For example, how can I reverse project the red point [1,0]?

Any ideas how that can be done?

Upvotes: 3

Views: 2046

Answers (2)

dww
dww

Reputation: 31452

I don't know if there's a reason you need to use mapproject to project in the fist place. If you can use spTransform instead, then this becomes easier because you can also use spTransform to reverse the same process.

Assuming that you do need to use mapproject, we can still use spTransform to convert points from your projected coordinate system into lat-long coordinates, but a little more fiddling is required to deal with the non-standard format of the mapproject projections I.e. the points are normalised to lie between -1 to 1 latitude and -2 to 2 longitude. In more standard map projections the lat/long are expressed in distances (usually meters).

So, first we can use spTransform to find out the conversion factors we need to convert the normalised mapproject coordinates into actual distances:

library(rgdal)
my.points = data.frame(x=c(0,180),y=c(90,0))
my.points = SpatialPoints(my.points, CRS('+proj=longlat'))
my.points = spTransform(my.points, CRS('+proj=moll'))
# SpatialPoints:
#             x       y
# [1,]        0 9020048
# [2,] 18040096       0
# Coordinate Reference System (CRS) arguments: +proj=moll +ellps=WGS84

Now we can use these references to convert from normalised mapproject coordinates into distances in meters:

my.points = data.frame(x=p$x * 18040096/2 , y=p$y * 9020048)
my.points = SpatialPoints(my.points, CRS('+proj=moll'))

And reproject these into lat/long geographic coordinates:

my.points = as.data.frame(spTransform(my.points, CRS('+proj=longlat')))

Finally we need to rotate these points by longitude, to undo the rotation that was performed in mapproject.

my.points$x = my.points$x + 200
my.points$x[my.points$x > 180] = my.points$x[my.points$x > 180] - 360

And lets check that it worked:

head(my.points)
#           x          y
# 1  75.77725  31.274368
# 2 -19.57400 -31.071065
# 3  79.78795 -24.639597
# 4  76.34576   1.863212
# 5  24.87848 -45.215432
# 6 -92.39700  23.068752

head(s)
#         lon        lat
# 1  75.77726  31.274367
# 2 -19.57400 -31.071065
# 3  79.78796 -24.639596
# 4  76.34576   1.863212
# 5  24.87849 -45.215431
# 6 -92.39700  23.068751

Upvotes: 3

Taeke
Taeke

Reputation: 321

If nothing out of the box is available, you could write a function of your own, along the lines of:

ref_project <- function(x, y) {
    long <- tibble(
        long = seq(-180, 180, 1),
        x = mapproject(long, rep(0, length(long)), projection = 'mollweide', orientation = c(90, 200, 0))$x
    )
    lat <- tibble(
        lat = seq(-90, 90, 1),
        x = mapproject(rep(0, length(lat)), lat, projection = 'mollweide', orientation = c(90, 200, 0))$y
    )

    return(c(long[which(abs(long$x - x) == min(abs(long$x - x))), 'long'],
             lat[which(abs(lat$x - y) == min(abs(lat$x - y))), 'lat']))
}

Upvotes: 2

Related Questions