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

Question

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)

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?

Is there some 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. — dww, Feb 11 '17 at 16:04
@dww There is no reason, I just do not know how to do the above using `sp`. So an `sp` solution is also welcome, if understandable for `sp` newbees :) — Mark Heckmann, Feb 11 '17 at 16:31
Ok, I added an answer already showing how to do it for mapproject. I can add also another answer for how to do the whole thing in `spTransform` without mapproject - if you also require this — dww, Feb 11 '17 at 16:43

dww · Accepted Answer · 2017-02-12T20:58:24.443

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

you're a star. I figure that in the long run its worthwhile to get used to sp but the learning curve appeared too steep for just a few tasks... Thanks a lot. :) — Mark Heckmann, Feb 11 '17 at 17:18
I just realised an error that i've now corrected. I do a lot of climate science, where we ususally use longitudes from zero to 360. Of course, normally people use -180 to +180. So I've chnaged the rotation to use `my.points$x[my.points$x > 180] = my.points$x[my.points$x > 180] - 360` to work with standard longitudes — dww, Feb 12 '17 at 21:01

score 2 · Answer 2 · answered Feb 11 '17 at 14:13

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']))
}

Basically, you suggest a lookup table for the projected coordinates. While a simple idea, it is very approximate, leading to slightly imprecise back-conversions. Of course, one could have a finder grained sequence. Maybe this will work. Props for this simple hands-on idea :) — Mark Heckmann, Feb 11 '17 at 14:36

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

2 Answers2