Skip to content

Instantly share code, notes, and snippets.

@J-Cleeland

J-Cleeland/polar_projection_polygons.Rmd

Last active February 12, 2017 18:08
Show Gist options
  • Save J-Cleeland/51ee85cd830ed54c31470ed04bdbf22a to your computer and use it in GitHub Desktop.
Save J-Cleeland/51ee85cd830ed54c31470ed04bdbf22a to your computer and use it in GitHub Desktop.
Download ZIP
For plotting circumpolar polygons on raster maps.
Raw
polar_projection_polygons.Rmd
---
title: "Polar projection polygons for light-mantled albatross distributions"
author: "Jaimie Cleeland, Michael Sumner"
date: "5/13/2016"
output: html_document
---
Load libraries.
```{r setup, message=FALSE}
library(graticule)
library(png)
library(raster)
library(RColorBrewer)
library(rgeos)
library(rworldmap)
```
Load world map and re-project to polar centric view.
```{r, echo=TRUE, message=FALSE, warning=FALSE}
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))
```
Option to add graticule lines to map.
```{r}
xx <- c(0,90,180,270,360); yy <- c(-80,-60,-40,-20)
g3 <- graticule(xx,yy,proj=CRS(pprj))
g4 <- graticule(xx,-20,proj=CRS(pprj))
g3labs1 <- graticule_labels(lons=180,xline=180,yline=-20,proj=CRS(pprj))
g3labs2 <- graticule_labels(lons=0,xline=180,yline=-20,proj=CRS(pprj))
```
Function to add longitudinal labels
```{r}
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)
par(p)
}
```
Create a marker for Maquarie Island (the breeding colony).
```{r}
colony <- as.data.frame(t(c(158.9452,-54.4957)))
colnames(colony) <- c("lon","lat")
coordinates(colony) <- c("lon","lat")
projection(colony) <- "+proj=longlat +datum=WGS84"
macca <- spTransform(colony,CRS(pprj))
```
Set plotting window
```{r}
par(mfrow=c(1,2), mar=c(0,1,0,0), family="serif")
```
Plot polygon on map (representing the distribution of breeding Macquarie Island light-mantled albatross)
```{r}
plot(w,col="white",border=FALSE)
plot(g3,add=TRUE,lty=3)
plot(macca,col="black",border=FALSE,add=TRUE,pch=19,cex=0.3)
g <- graticule(c(135, 200), c(-70, -40), proj = "+proj=laea +lat_0=-90 +datum=WGS84")
pg <- spTransform(g, CRS(pprj))
plot(gPolygonize(pg), col = "#1B9E774D", add=T, border=FALSE)
plot(w,col="darkgrey",border=FALSE,add=TRUE)
plot(g4,add=TRUE)
pltg()
mtext('a) Breeding',side=3,cex=1,line=-1,col="grey40")
```
And if you need to know what a light-mantled albatross looks like:
```{r}
ima <- readPNG("/user/jaimiec/MI_Albatross/plots/Albatross_heads/Light-mantled.png")
par(new=TRUE, mar=c(0,0,0,0))
plot(1,type='n',xlim=0:1,ylim=0:1,main="",axes = FALSE,xlab="", ylab="")
# xleft,ybottom,xright,ytop
rasterImage(ima,0,0.8,0.3,1)
```
Finally plot nonbreeding distribution (circumpolar polygon)
```{r, include=TRUE, eval=TRUE}
plot(w,col="white",border=FALSE)
plot(g3,add=TRUE,lty=3)
plot(macca,col="black",border=FALSE,add=TRUE,pch=19,cex=0.3)
g <- graticule(c(0, 180), c(-60, -40), proj = "+proj=laea +lat_0=-90 +datum=WGS84")
pg <- spTransform(g, CRS(pprj))
plot(gPolygonize(pg), col = "#1B9E774D", add=T, border=FALSE)
g <- graticule(c(180, 360), c(-40,-60), proj = "+proj=laea +lat_0=-90 +datum=WGS84")
pg <- spTransform(g, CRS(pprj))
plot(gPolygonize(pg), col = "#1B9E774D", add=T, border=FALSE)
plot(w,col="darkgrey",border=FALSE,add=TRUE)
plot(g4,add=TRUE)
pltg()
mtext('b) Nonbreeding',side=3,cex=1,line=-1,col="grey40")
```
@J-Cleeland
Copy link
Author

lma_distributions

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment