Create a map detailing study site and oceanographic features

[nmbool]

Hi Mike,

I'm posting this one for someone else (they're having issues with GitHub). How can we create a really nice map to show a study site, including oceanographic features etc?

thanks

[mdsumner]

These days I would start with leaflet. https://rstudio.github.io/leaflet/

First though: "what oceanographic features"? Unfortunately there's no obvious collection of the main things, but we can look at oce and marmap and there are some others on CRAN.

Note we have orsifronts in a package, and leaflet usings "magittr/dplyr" style piping.

library(orsifronts)  ## install. packages("orsifronts")
library(leaflet)        ## install.packages("leaflet")

## orsifronts is a SpatialPolygonsDataFrame (could be read via shapefile, or contoured from other data)
data(orsifronts)
fronts_as_points <- as(orsifronts, "SpatialPointsDataFrame")
fronts_as_points <- fronts_as_points[sample(nrow(fronts_as_points), 20), ]

leaflet() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% 
  addPolylines(data = rmapshaper::ms_simplify(orsifronts), 
               color = ~colorFactor(viridis::viridis(nrow(orsifronts)), front)(front)) %>% 
  addMarkers(data = fronts_as_points, label = ~front,  labelOptions = labelOptions(noHide = TRUE))

(I had to simplify the fronts layer to get it to work, and I don't know why. )

We can talk about how to plot this more traditionally, but I think the key is getting familiarity with plotting map stuff generally. In R we need sp and raster, everything else comes from those kinds of objects.

[J-Cleeland]

---

title: "Making a general Southern Ocean map" author: "Jaimie Cleeland, jaimie.cleeland@utas.edu.au" date: "3/4/2017" output: html_document

Load libraries.

rm(list=ls())
library(graticule)
library(raadtools)
library(raster)
library(RColorBrewer)
library(rworldmap)
library(sp)

Load world map, crop it and re-project to polar centric view.

data(countriesLow)
raster.map <- raster(xmn=-180, xmx=180, ymn=-90, ymx=-20)
mp <- crop(countriesLow, extent(raster.map))
pprj <- "+proj=laea +lat_0=-60 +lon_0=180 +datum=WGS84 +ellps=WGS84 +no_defs +towgs84=0, 0, 0"
w <- spTransform(mp, CRS(pprj))

Grab topo for bathymetric contours, remove contours on land (>0m) and re-project to polar centric view.

topo1 <- readtopo("etopo2", xylim=extent(raster.map))
topo1[topo1 > 0 ] <- 0
cl1 <- rasterToContour(aggregate(topo1, fact=16, fun=mean))
pcl1 <- spTransform(cl1, CRS(pprj))

Load fronts and re-project to polar centric view.

front <- frontsmap(map="orsi")
front <- spTransform(front, CRS(pprj))

Creat graticule lines to add to map.

xx <- c(0, 90, 180, 270, 360); yy <- c(-80, -60, -40, -20)
g3 <- graticule(xx, yy, proj=projection(pprj))
g4 <- graticule(xx, -20, proj=projection(pprj))
g3labs1 <- graticule_labels(lons=180, xline=180, yline=-20, proj=projection(pprj))
g3labs2 <- graticule_labels(lons=0, xline=180, yline=-20, proj=projection(pprj))
g4labs <- graticule_labels(lats=c(-40, -60, -20), yline = -20, xline=0, proj=projection(pprj))

Function to add latitude and longitudinal labels.

pltg <- function() {
  p <- par(xpd=NA)
  text(coordinates(g3labs1[g3labs1$islon, ]), lab=parse(text=g3labs1$lab[g3labs1$islon]), pos=3, cex=0.8)
  text(coordinates(g3labs2[g3labs2$islon, ]), lab=parse(text=g3labs2$lab[g3labs2$islon]), pos=1, cex=0.8)
  text(coordinates(g4labs[!g4labs$islon, ]), lab=parse(text=g4labs$lab[!g4labs$islon]), pos=3, cex=0.8)
 par(p)
}

Create a dataframe with all the relevant geographic locations you would like to include on your map, including a logical (Y/N) as to whether we want a marker added or not and a left, centre or right adjustment (adj).

colony <- data.frame(lon=c(158.945, 160.431, 160.5, -120, -25.673, 81.826, -71.383, -155.847, -90, 175, 169.16, -38.03, 37.866, 9.85, -62.03, 69.16, 51.21, 73.50, -59.52), lat=c(-54.495, -40.858, -50.6, -31, -31, -33, -59.914, -40, -40, -49, -52.51, -54.00, -46.88, -25.21, -48.00, -49.25, -46.322, -53.08, -51.79), name=c("Macquarie Is.", "Tasman\nSea", "Macquarie\nRidge", "Pacific\nOcean", "Atlantic\nOcean", "Indian\nOcean", "Drake Passage", "South-west\n Pacific Basin", "South-east\n Pacific Basin", "Campbell\nPlateau", "Campbell Is.", "Bird Is.\n(South Georgia)", "Marion and\nPrince Edward Is.", "Benguela\nCurrent", "Patagonian\nShelf", "Kerguelen Is.", "Crozet Is.", "Heard and\nMacDonald Is.", "Falkland Is."), marker=c("y", "n", "n", "n", "n", "n", "n", "n", "n", "n", "y", "y", "y", "n", "n", "y", "y", "y", "y"), adj = c(0, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 1))

Now here is the ugly bit. To offset the labels from the markers, we adjust the lat/lons of the subsetted labels. Then re-project them to have a polar-centric projection.

lab_pos <- colony[grep("Is.", colony$name), ]
lab_pos$lon[lab_pos$name %in% c("Macquarie Is.", "Campbell Is.")] <- lab_pos$lon[lab_pos$name %in% c("Macquarie Is.", "Campbell Is.")] + 2
lab_pos$lat[lab_pos$name %in% c("Macquarie Is.", "Campbell Is.")] <- lab_pos$lat[lab_pos$name %in% c("Macquarie Is.", "Campbell Is.")] - 0.5
lab_pos$lon[lab_pos$name %in% c("Marion and\nPrince Edward Is.")] <- lab_pos$lon[lab_pos$name %in% c("Marion and\nPrince Edward Is.")]  - 1.0
lab_pos$lat[lab_pos$name %in% c("Marion and\nPrince Edward Is.")] <- lab_pos$lat[lab_pos$name %in% c("Marion and\nPrince Edward Is.")]  - 0.5
lab_pos$lon[lab_pos$name %in% c("Crozet Is.")] <- lab_pos$lon[lab_pos$name %in% c("Crozet Is.")]  - 1.0
lab_pos$lat[lab_pos$name %in% c("Crozet Is.")] <- lab_pos$lat[lab_pos$name %in% c("Crozet Is.")]  - 1.0
lab_pos$lon[lab_pos$name %in% c("Kerguelen Is.")] <- lab_pos$lon[lab_pos$name %in% c("Kerguelen Is.")]  - 2.5
lab_pos$lat[lab_pos$name %in% c("Kerguelen Is.")] <- lab_pos$lat[lab_pos$name %in% c("Kerguelen Is.")]  - 0.5
lab_pos$lat[lab_pos$name %in% c("Heard and\nMacDonald Is.")] <- lab_pos$lat[lab_pos$name %in% c("Heard and\nMacDonald Is.")]  - 1.0
lab_pos$lon[lab_pos$name %in% c("Bird Is.\n(South Georgia)")] <- lab_pos$lon[lab_pos$name %in% c("Bird Is.\n(South Georgia)")]  - 1.0
lab_pos$lat[lab_pos$name %in% c("Bird Is.\n(South Georgia)")] <- lab_pos$lat[lab_pos$name %in% c("Bird Is.\n(South Georgia)")]  + 1.0
lab_pos$lon[lab_pos$name %in% c("Falkland Is.")] <- lab_pos$lon[lab_pos$name %in% c("Falkland Is.")]  + 1.0
lab_pos$lat[lab_pos$name %in% c("Falkland Is.")] <- lab_pos$lat[lab_pos$name %in% c("Falkland Is.")]  - 0.5

Make labels and markers SpatialPoints. Then re-project them to have a polar-centric projection.

coordinates(colony) <- c("lon", "lat")
projection(colony) <- "+proj=longlat +datum=WGS84"
geog <- spTransform(colony, CRS(pprj))
coordinates(lab_pos) <- c("lon", "lat")
projection(lab_pos) <- "+proj=longlat +datum=WGS84"
lab_pos <- spTransform(lab_pos, CRS(pprj))

We're going to save this plot out as a high res tiff file.

tiff(file=paste0("~/MI_Albatross/plots/SuppFigure2_geographic_locations.tiff"), width=7.5, height=7.5, units="in", res=300)

Ok now for the plot

#Set the plotting parameters
par(family="serif", bty="n", mar=c(1, 1, 1, 0), font=2)
#plot blank map to appropriately set bounds
plot(w, col="white", border=FALSE)
#add contours
plot(pcl1, add=TRUE, col=grey(0.7, alpha=0.3))
#Add graticule lines
plot(g3, add=TRUE, lty=3)
#Choose some colours for the Southern Ocean fronts
fcol <- brewer.pal(n=9, name="Greys")[4:8]
#Plot fronts
plot(front[1, 1], add=TRUE, lty=1, lwd=1, col=fcol[1]) #STF
plot(front[2, 1], add=TRUE, lty=1, lwd=1, col=fcol[2]) #SAF
plot(front[3, 1], add=TRUE, lty=1, lwd=1, col=fcol[3]) #PF
plot(front[4, 1], add=TRUE, lty=1, lwd=1, col=fcol[4]) #Southern Antarctic circumpolar current front
plot(front[5, 1], add=TRUE, lty=1, lwd=1, col=fcol[5]) #Southern boundary of ACC
#Choose some colours for the labels we will add to the plot
col <- brewer.pal(3, "Dark2")
#Add world map
plot(w, col="darkgrey", border=FALSE, add=TRUE)
#Add markers
plot(geog[geog$marker=="y", ], col=col[3], border=FALSE, add=TRUE, pch=19, cex=0.5)
#Add labels for all those we want left justified
text(lab_pos[lab_pos$adj==0, ], labels=lab_pos$name[lab_pos$adj==0], cex= 0.75, adj=0, col=col[3])
#Add labels for all those we want right justified
text(lab_pos[lab_pos$adj==1, ], labels=lab_pos$name[lab_pos$adj==1], cex= 0.75, adj=1, col=col[3])
#Add a label that we want in a different colour and right justified
text(geog[geog$name %in% c("Macquarie\nRidge"), ], labels=geog$name[geog$name %in% c("Macquarie\nRidge")], cex= 0.75, adj=1, col=col[2])
#Add in the rest of the feature we want in black text
text(geog[geog$name %in% c("Tasman\nSea", "Drake Passage", "South-west\n Pacific Basin", "South-east\n Pacific Basin", "Campbell\nPlateau", "Benguela\nCurrent", "Patagonian\nShelf"), ], labels=geog$name[geog$name %in% c("Tasman\nSea", "Drake Passage", "South-west\n Pacific Basin", "South-east\n Pacific Basin", "Campbell\nPlateau", "Benguela\nCurrent", "Patagonian\nShelf")], cex= 0.75, adj=0.5, col="black")
#Add the ocean names in green
text(geog[geog$name %in% c("Pacific\nOcean", "Atlantic\nOcean", "Indian\nOcean"), ], labels=geog$name[geog$name %in% c("Pacific\nOcean", "Atlantic\nOcean", "Indian\nOcean")], cex= 0.75, adj=0.5, col=col[1])
#Add in a line to signify Macquarie Ridge
arrows(x0=-1341759, x1=-1016381, y0=500292.8, y1=1408583, length=0, angle=0, code=1, col=col[2], lty=1, lwd=0.8)
#Add a graticule border
plot(g4, add=TRUE)
#Add labels
pltg()

Lets make a legend to finish it off

#Set up par(new=T) to plot over the top of your map
par(new=T, mar=c(0, 0, 0, 0), font=1)
#Add empty plot
plot(1, type="n", axes=F, xlab="", ylab="")
#Add headed
mtext('Fronts', side=1, cex=1, line=-5.3, col="grey40", at = 0.64)
#Add legend, referencing the coordinates of the empty plot
legend(x = 0.57, y = 0.67, c("Subtropical Front", "Subantarctic Front", "Polar Front", "Southern ACC Front", "Southern Boundary of ACC"), lty=c(1, 1, 1, 1, 1), col=fcol[1:5], bty="n", cex=0.7, inset=0.5)

Close plotting window to save.

dev.off()

[J-Cleeland]

[J-Cleeland]

I've had a crack at a southern ocean map that may be helpful. But I think the leaflet option is going to be a better way to go.

[mdsumner]

This is wonderful

[J-Cleeland]

Saturday night R fever! Hopefully, a few in the lab can use it...and improve it. Time for a beer.

[mdsumner]

I made some small changes and added it as .Rmd file to this repo, you can see the changes using the history feature:

https://github.com/DataScienceHobart/2017-02-03-graphics-imas/commit/6d97298405cde09f3e5f8968b77efd89adbc41a4

They are

• use orsifronts package rather than frontsmap (I'll remove this from raadtools eventually)
• plot the fronts in vectorized form, easiest is to add the colours as data to the object
• I started to convert the label-mods to "tribble" form, but it's not that important so I abandoned it

We should capture this as a tool to trot out more easily, and I'm happy to do that. I think this is a good example of the kinds of shared-needs we can work on together.

We might be able to get the labels out of antanym, and @raymondben will have better ideas about that.

Nice one!

[mdsumner]

Notes for myself, use this for rewriting the lab-mods: https://gist.github.com/hrbrmstr/8c7862092681d06e6535ab38d03074bf (it's not relevant to the map, but it's a tool to help rewrite the code in a more readable form). tribble ftw

[J-Cleeland]

Hey @mdsumner, Mind my noob status here. Do I have write access to this code? Still learning GitHub. I'd like to add a polygon to it. This is a common thing we do in the lab and could be included in this example.

[J-Cleeland]

@mdsumner Couldn't change the projection as per recommendations without it falling over...

probably use lat_0=-90 for more general maps, it's a bit off-centre

pprj <- "+proj=laea +lat_0=-60 +lon_0=180 +datum=WGS84 +ellps=WGS84 +no_defs +towgs84=0,0,0"

For others- Up to date map is here: https://github.com/DataScienceHobart/2017-02-03-graphics-imas/blob/master/making-general-socean-map.rmd

Can we directly raise issues on that code?

[J-Cleeland]

Hey @mdsumner, just wondering if you could offer up your fix to the cropping problem when world map is reprojected into polar stereographic form? Here: https://github.com/DataScienceHobart/2017-02-03-graphics-imas/blob/master/making-general-socean-map.rmd

[J-Cleeland]

Case in point: annoying crop on central Aus

[mdsumner]

Try this, crop after reproject

https://gist.github.com/mdsumner/3da004f9df5ad30811641bbcbf3fb9e0

The crop should insert a segmentized arc (from the circular buffer) to carry the parallel at the northern edge.