MadisonDatumChanges.R

#  This executes the map figure shown for GEOG378 Section 2.1 using R.
#  The figure demonstrates the impact that changing datum may have on lat/long coordinates.

#  Load the package 'rgdal'.  If this package is missing on your system you can load it
#  using the call 'install.packages('rgdal')'

library(rgdal)

# This is a built-in GDAL command, implemented in R using the rgdal package loaded above.
datum_list  <- projInfo(type='datum')

#  These are the coordinates 
madison_coords <- data.frame(long = 89.38269, lat = 43.05476)

#  To define a SpatialPoints object in R we must first indicate the coordinates and the
#  projection.

coordinates(madison_coords) <- ~ long + lat
proj4string(madison_coords) <- CRS('+proj=longlat +datum=WGS84')

#  This is the list into which we stick the new coordinates once reprojected.
output_coords <- data.frame(datum = rep(NA, 10),
                            long  = rep(NA, 10),
                            lat   = rep(NA, 10))

#  Here is our for loop:
for(i in 1:10){
  # We take the new datum and use the R function spTransform to convert the coordinates to the new datum.
  # the CRS - coordinate reference system - is made up of the projection and the datum, pasted together,
  # this is then transformed, and then we take the coordinates and put them into columns 2 and three of the 'i'th row
  # of the data.frame output_coords.
  output_coords[i,1]   <- datum_list[i,1]
  output_coords[i,2:3] <-coordinates(spTransform(madison_coords, CRSobj = CRS(paste0('+proj=longlat +datum=',datum_list[i,1]))))
}

#  I need to convert the output_coordinates from a "W" to a 180^o^ system.
output_coords[,2] <- -output_coords[,2]

#  This is sneaky, I like this leaflet package, but it's a bit unweildy to use.You can install it as any other package:
#  install.packages('leaflet')

library(leaflet)
leaflet(output_coords) %>% addTiles() %>% addTiles() %>% addCircles(radius = 1, color = 'red', opacity = 1) %>% addMarkers(popup=output_coords[,1])