Making nice pdf maps is giving me trouble. I want to export as a pdf
using. pdf("mariana_chlor_log2.pdf", width=4, height = 3,
layout(matrix(c(1,2,3,4), nrow=1, ncol=2), widths=c(3,1), heights=4) ).
I end up getting an empty 2 panel layout. Any ideas what I can do to
correct this?
pdf("mariana_chlor_log2.pdf", width=4, height = 3,
layout(matrix(c(1,2,3,4), nrow=1, ncol=2), widths=c(3,1), heights=4) )
#layout(matrix(c(1,2,3,4), nrow=1, ncol=2), widths=c(3,1),
heights=4)
layout.show(2)
par(mar=c(3,3,3,1))
land <- maps::map("world2",
fill = TRUE, col='grey',xlim = xlim, ylim = ylim, plot=FALSE,
xaxs='i',yaxs='i',asp=1, cex.axis = 3, proj4string =
CRS("+proj=longlat +datum=WGS84"))
ids <- sapply(strsplit(land$names, ":"), function(x)
x[1])
bPols <- map2SpatialPolygons(land, IDs = ids,
proj4string = CRS("+proj=longlat +datum=WGS84"))
plot(bPols, col = "grey", axes = FALSE, xlim = xlim, ylim =
ylim, xaxs='i',yaxs='i',asp=1, cex.axis = 3, main="Chlor
a")
x = seq(140, 150, 1)
axis(1, x, sapply(paste(x,
"E", sep = ""), function(x) x[1]))
y = seq(12, 22, 2)
axis(2, las = 1, y, sapply(paste(y,
"N", sep = ""), function(x) x[1]))
box()
temp.range=list()
for (j in 1:10) {
print(j)
I=which(mariana[,4]==j) #select the feature #
poly=mariana[I,1:2]
#write.csv(mariana,'mariana.csv')
#We are only interested in column 2 and 3 of pnmn.csv,
which correspond to the longitudes and latitudes of the boundary,
respectively.
#poly=read.csv('mariana.csv',header=TRUE)[,2:3]
poly=data.frame(poly)
names(poly) <-
c("lon","lat")
#We need to convert the
longitudes to 0-360? to make it easier to work across the
dateline.
I=which(poly$lon<0)
if (length(I)>0)
poly$lon[I]=poly$lon[I]+360
ERDDAP_Node="https://oceanwatch.pifsc.noaa.gov/erddap/"
dataInfo <- rerddap::info('esa-cci-chla-2009-2018-clim-v4-2',
url=ERDDAP_Node)
#Let's look at the variables in our
dataset:
dataInfo$variable$variable_name
#[1]
"analysed_sst" "sea_ice_fraction"
#We are
interested in sea surface temperature data, so we need to select the
fist variable:
#analysed_sst
parameter=dataInfo$variable$variable_name[1]
xcoord
<- poly$lon
ycoord <- poly$lat
chl
<- rxtractogon (dataInfo, parameter=parameter, xcoord=xcoord,
ycoord=ycoord, tcoord=tcoord)
chl$time ##only two time steps
here! why not all 12
#We need to transform our
extracted chl data for that time step into a matrix for the image
function to work properly. This chl dataset is also in Kelvin degrees,
so let-s convert it to Celsius degree:
test=as.matrix(chl$chlor_a[,,i])
#want to get annual mean of
chl for 2019
#mean_chl=as.matrix(apply(chl$chlor_a,c(1,2),mean,na.rm=TRUE))
#chl.yr=apply(chl[,,1:12],c(1,2),mean,na.rm=TRUE)
test=as.matrix(log(chl$chlor_a[,,i]))
#Let's
make pretty axes:
#We'll add the color plot on top of the
map:
par(new=TRUE)
image(chl$longitude,chl$latitude,test,col=c,breaks=breaks,xlab='',ylab='',
axes=FALSE, xlim = xlim, ylim = ylim,xaxs='i',yaxs='i',asp=1,
las=1)
#plot the boundary lines
lines(poly$lon,poly$lat,lwd=1)
temprange<-c(min=min(chl$chlor_a, na.rm=TRUE),max=max(chl$chlor_a,
na.rm=TRUE),mean=mean(chl$chlor_a, na.rm=TRUE), L=j) #tryed
total=sum(chl$chlor_a, na.rm=TRUE) here did not work
temp.range[[j]]<-temprange
}
#Plot color scale using the scale.R file:
source('scale.R')
#par(mar=c(3,1,3,3))
image.scale(test,
col=c, breaks=breaks, horiz=FALSE,
yaxt="n",xlab='',ylab='',main='log chlor a')
axis(4,las=1)
box()
dev.off()
