auto-detect coordinate reference system, based on coordinates in GPX file?

0
votes

I'm building a Shiny app that color-codes GPX track logs based on the local slope at each point.

It's based extensively on https://rpubs.com/chrisbrunsdon/hiking

To calculate the "run" part of slope = rise/rune, I'm converting from latitude/longitude data into a X,Y grid (in meters) with sf::st_transform. One of the arguments for that function is crs, or "coordinate reference system".

Up to now, I've been testing with GPX files I gathered near my home in southeastern Pennsylvania, so I've been using EPSG:2272 as my CRS.

To make this useful to anyone with logs form anywhere in the world, I'd like to auto-detect the most appropriate CRS based on the centroid of the points in the plotted track. Is there some canned function for doing that?

2
rgisgeospatial

2 Answers

1
votes

You can use the UTM projection.
Basically, retrieve the appropriate zone number and letter for the centroid, convert the track to that zone and perform the calculations.

//pseudocode
utm.from_latlon(51.2, 7.5)
//EASTING, NORTHING, ZONE NUMBER, ZONE LETTER
395201.3103811303, 5673135.241182375, 32, "U"

The zone letters are actually latitude bands and won't be needed if you're working with EPSG codes.
To "manually" calculate them use:

zone_num <- floor((longitude + 180) / 6) + 1  #each zone is 6 degrees wide
hemisphere <- if (latitude >= 0) "northern" else "southern"
epsg <- 32600 + zone_num
if (hemisphere == "southern") {
    epsg <- epsg + 100
}
0
votes

I found how-to-get-appropriate-crs-for-a-position-specified-in-lat-lon-coordinates (for the related EPSG crs, not UTM zone per se), and re-wrote it in R, assuming the user wants to analyze the first track segment in file gpx.file. 

rg.result <- readGPX(gpx.file)

outer.track.list <- rg.result$tracks
inner.track.list <- outer.track.list[[1]]
track.frame <- inner.track.list[[1]]

tf.avg.lat <- mean(track.frame$lat)
tf.avg.lon <- mean(track.frame$lon)

EPSG <- 32700-round((45+tf.avg.lat)/90,0)*100+round((183+tf.avg.lon)/6,0)

# make a spatial frame, based on GPX's use of WGS84/EPSG 4326
coords <-
  st_as_sf(track.frame,
           coords = c("lon", "lat"),
           crs = 4326)

# project that according to the EPSG crs determined above
st_transformed <-
  st_transform(coords$geometry, crs = EPSG)

# `xy` will be a matrix of positions on a grid, in meters
xy <- st_coordinates(st_transformed)