library(rgdal)
library(rworldxtra)
library(mapplots)
library(raster)
library(RColorBrewer)
#
# shapefiles etc
#
rfa <-
readOGR(file.path(
system.file("shapefiles", package = "IBTSindices"),
"Roundfish_shapefiles"
))
rfa_labelpos <- coordinates(rfa)
rfa_labelpos[9,1] <- rfa_labelpos[9,1]+0.4 #adjust position for area 3
rfa_labelpos[6,1] <- rfa_labelpos[6,1]+0.2 #adjust position for area 5
rfa_labelpos[6,2] <- rfa_labelpos[6,2]+0.3 #adjust position for area 5
rfa_labelpos[6,2] <- rfa_labelpos[6,2]+0.2 #adjust position for area 5
rfa_labelpos[5,2] <- rfa_labelpos[5,2]+0.2 #adjust position for area 6
rfa_labelpos[9,1] <- rfa_labelpos[9,1]+0.3 #adjust position for area 3
rfa_labelpos[2,1] <- rfa_labelpos[2,1] #adjust position for area 9
rfa_labelpos[2,2] <- rfa_labelpos[2,2]+0.1 #adjust position for area 9
ices_rectangles <- readOGR(file.path(system.file("shapefiles", package = "IBTSindices"), "ICES_StatRec_mapto_ICES_Areas"), "StatRec_map_Areas_Full_20170124")
ices_midpoints <- SpatialPoints(as.data.frame(coordinates(ices_rectangles)), proj4string=CRS(proj4string(ices_rectangles)))
ibts_rectangles <- ices_rectangles[!is.na(unlist(sp::over(ices_midpoints, rfa))),]
ns <- readRDS(file.path(system.file("bathymetry", package = "IBTSindices"), "wgs84_raster_north_sea_crop", "ns_halfres.rds"))*-1
data(countriesHigh)
map <- countriesHigh
#Needs fiddling when final output dimensions are decided
map_labels <- c("Denm.", "France", "Germany", "Netherl.", "Norway", "UK Scotl.", "Sweden", "UK Engl.")
map_labelpoints <- t(matrix(c(
10.9, 55.6,
2.5, 49.3,
10, 51,
6.9, 52.28,
9.4, 60.7,
-5, 57,
14, 59,
-1.9, 52.4), ncol=length(map_labels), nrow=2))
#
# colors
#
batpal <- brewer.pal(9, "Blues")
batpal <- c(batpal[1:(length(batpal)-3)], rep(batpal[7],4),rep(batpal[8],4), rep(batpal[9],4))
plot_ibts_map <- function(polygons = rfa, labelcol="AreaName", labelpos=rfa_labelpos, polygonscol="black", cex.polygons=0.8, rectangles=ibts_rectangles, rectanglecol="lightgray", landmap=map, landcol="lightgrey", landbordercol="darkgrey", landlabels=map_labels, land_labelpoints=map_labelpoints, land_labelcol="black", cex.landlabels=0.6, seacol="white", bat=ns, batcol=batpal, xlim=c(-8,18), ylim=c(49,62.5), xlab="Longitude", ylab="Latitude", cex.legend=1, main=""){
if (!is.null(bat)){
ext <- as(extent(xlim[1], xlim[2], ylim[1], ylim[2]), 'SpatialPolygons')
plot(crop(bat, ext), col=batcol, xlab=xlab, ylab=ylab, legend.args=list(text="depth (m)", cex=cex.legend), main=main)
}
else{
basemap(
xlim = xlim,
ylim = ylim,
bg = seacol,
xlab=xlab, ylab=ylab,
main=main
)
}
plot(landmap, col = landcol, add = T, border="darkgrey")
if (!is.null(polygons)) {
if (!is.null(rectangles)){
plot(rectangles, add=T, border=rectanglecol)
}
plot(polygons, add = T, border=polygonscol)
if (!is.null(labelcol)){
if (is.null(labelpos)){
labelpos <- coordinates(polygons)
}
text(labelpos, as.character(polygons@data[[labelcol]]), col=polygonscol, cex=cex.polygons)
}
if (!is.null(land_labelpoints)){
text(land_labelpoints, landlabels, cex=cex.landlabels, col=land_labelcol)
}
}
}
plot_on_pdf <- function(filename){
pdf(filename, width=3.35*2, height=3.35)
old.par <- par(no.readonly = T)
par(mfrow=c(1,2), mar=c(5.1, 4.1, 4.1, 4.1))
plot_ibts_map(bat = NULL, land_labelpoints = NULL)
par(mar=c(5.1, 1.1, 4.1, 7.1))
plot_ibts_map(rectangles = NULL, ylab="", labelcol = NULL)
par(old.par)
dev.off()
}
#plot_on_pdf("surveyarea.pdf")
# could not figure out how to make R set meta information correct.
# The right number of pixels are exported, but resolution is set to 72 dpi, so it displays too wide in viewers.
# may edit this in imaging program after export:
# E.g, with Gimp: Image > Scale Image
plot_on_tiff <- function(filename){
tiff(filename, width=17, height=7.1, res = 1200, units="cm")
old.par <- par(no.readonly = T)
par(mfrow=c(1,2), mar=c(5.1, 4.1, 1, 3.1))
plot_ibts_map(bat = NULL, land_labelpoints = NULL)
par(mar=c(5.1, 1.1, 1, 7.1))
plot_ibts_map(rectangles = NULL, ylab="", labelcol = NULL)
par(old.par)
dev.off()
}
plot_on_tiff("surveyarea_1200dpi_two_column.tiff")
oldpar <- par(no.readonly = T)
par(cex.axis=1.5, cex.lab=1.5, mgp=c(3, 2, 0))
plot_ibts_map(cex.legend=1.5, bat=NULL)
par(oldpar)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.