lemur_fin_analysis/fig_plot_raw_data.R at master · densitymodelling/lemur_fin_analysis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
# plot transects, observations, bathymetry

library(ggplot2)
library(mapdata)
library(cmocean)
library(mgcv)


load("RData/0_format_aux_data.RData")
load("RData/1_model_and_data.RData")

# setup depth data
depthdat <- read.csv("data/Depth_ocean_CCEgrid.csv")
x <- depthdat$lon180
y <- depthdat$lat
ind <- inSide(cce_poly, x, y)
depthdat <- depthdat[ind,]

# make some breaks
depth_breaks <- quantile(depthdat$Depth_ETOPO1, prob=seq(0,1, len=8),na.rm=TRUE)

# setup palette for plotting
pal <- cmocean('deep')(length(depth_breaks)+2)

# map
w <- map_data("worldHires", ylim = range(cce_poly$y), xlim = range(cce_poly$x))

# plot bathymetry with effort and observations on map
p_dat <- ggplot() +
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_contour_filled(data=depthdat, aes(x=lon180, y=lat, z=Depth_ETOPO1),
                      breaks=depth_breaks) +
  scale_fill_discrete(type=pal) +
  geom_point(aes(y=mlat, x=mlon), size=0.2, data=fin) +
  geom_point(aes(y=mlat, x=mlon), colour="#d95f0e", size=0.4,
             data=subset(fin, count>0)) +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.text = element_text(size=12),
        legend.title = element_text(size=12),
        axis.text = element_text(size=12)) +
  labs(x="", y="", fill="Depth\n(m)") +
  coord_map(ylim = c(30.2, 48.0), xlim = c(-131.0, -117.3))
print(p_dat)

ggsave(p_dat, file="figures/rawdat.pdf", width=7, height=9)


# make another plot for figure that has e.g. prediction grid
source("support_scripts/prepare_preds.R")
pp <- prepare_preds("data/CCE_0.1deg_2008-09-13.csv", cce_poly, pred_areas)

# just plot mixed layer depth
p_ild <- ggplot() +
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_tile(data=pp, aes(x = mlon, y = mlat, fill = ild))+
  scale_fill_gradientn(colours=pal) +
  theme_minimal() +
  labs(x="", y="", fill="MLD") +
  coord_fixed(1.3, ylim = range(cce_poly$y), xlim = range(cce_poly$x))
print(p_ild)

ggsave(p_ild, file="figures/raster_ex.pdf", width=5, height=7)