Convert longitude latitude coordinates (WGS) to grid with equidistant axes (in given area)

Question

I have a lot of geocoordinates in tow data sets and want to run a nearest-neighbor-search. I came across the package 'RANN' and the function nn2(x,y) runs really fast.

Now there is the problem, that of course in the area of London a degree to the north is a quite a longer way then a degree to the west.

My idea now was to convert the location coordinates to some grid where one step in the direction of x is nearly the same as one step in the direction of y. The area is London (Center -0.1045,51.489). How can I perform this transformation?

library(RANN)

xyunf <- structure(c(-0.19117, -0.173862, -0.187623, -0.187623, -0.192366, 
-0.176224, 51.489096, 51.482442, 51.50226, 51.50226, 51.491632, 
51.495429), .Dim = c(6L, 2L), .Dimnames = list(c("1", "2", "3", 
"4", "6", "7"), c("Longitude", "Latitude")))
xyosm <- structure(c(-0.1966434, -0.1097162, -0.2023061, -0.198467, -0.4804301, 
-0.4286548, 51.6511198, 51.6134576, 51.6042042, 51.5186019, 51.3757395, 
51.3351355), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("lon", 
"lat")))

res <- nn2(data=xyunf, query=xyosm, k=1)
res$nn.dists
res$nn.idx

Answer 1

If you read the R Spatial Task view you can find out all about Spatial objects - these are points, grids, lines, or polygons that can have a coordinate reference system associated.

Once you've got those, you can use spTransform to convert between coordinate systems. So to convert a dataframe with lat/lon items from lat-long to Ordnance Survey Grid Coordinates:

coordinates(ptsLL) = ~Longitude+Latitude # turns a dataframe into a SpatialPointsDataFrame
proj4string(ptsLL) = CRS("+init=epsg:4326")  # tells it to be lat-long WGS84
ptsOS = spTransform(ptsLL, CRS("+init=epsg:27700")) # converts to UK Grid system
ptsOS = pts@coords

Now the thing about epsg:27700 is that it is a square grid in metres, so work with that. If you need to convert back to lat-long, spTransform again.

There are other square grid coordinate systems for other parts of the world, so don't use this for Australia!

Answer 2

You could get a pretty good estimate for north-south distances for temperate latitudes using just 111.1 km/degree. For longitude (east-west):

mtrs_degr_long <- function(long) { a=6378137.0 ; b=6356752.3142 
         e.sqr <- a^2/b^2 -1; long=long*2*pi/360
         pi*a*cos(long)/( 180*(1-e.sqr*sin(long)^2)^(1/2) ) }

{See the Wikipedia article.)

Now you can multiply your longitude values by mtrs_degr_lat(51.5)/111100 and then you will have adjusted the scale so relative distances are correct. If you convert both to meters you can get absolute distances. Given the relatively small range, the errors induced by the non-rectilinear coordinate system will be very small.