Skip to content

Instantly share code, notes, and snippets.

@SimonGoring
Last active September 12, 2015 04:35
Show Gist options
  • Save SimonGoring/e1698c35b43227a08ad0 to your computer and use it in GitHub Desktop.
Save SimonGoring/e1698c35b43227a08ad0 to your computer and use it in GitHub Desktop.
# 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])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment