Skip to content

Commit

Permalink
updating scripts from 2 years ago lol
Browse files Browse the repository at this point in the history
  • Loading branch information
@anthony-meza
anthony-meza committed Feb 16, 2024
1 parent 15ac5b2 commit e44fc97
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 136 deletions.
Binary file modified .DS_Store
View file
Binary file not shown.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,5 @@
*.jl.cov
*.jl.mem
/docs/build/
/plots/
*~
68 changes: 34 additions & 34 deletions scripts/deltaC13_regen.jl
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
#=
% Example: Compute the spatial distribution of δ¹³C %
=#

using Revise
using TMI, PyPlot, PyCall, Statistics
using TMI, PyPlot, PyCall, Statistics, DrWatson

pygui(false)
mkpath(plotsdir())

TMIversion = "modern_90x45x33_G14_v2"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
Expand All @@ -23,60 +22,61 @@ d_NATL = getRegion(TMIversion,"d_NATL")
#plot the distribution at the surface
fig = figure()
ax = fig.add_subplot(projection = TMI.cartopy.crs.PlateCarree())
ax.contourf(γ.lon,γ.lat,δ¹³C[:, :,lvl]' .* d_NATL' ,
cax = ax.contourf(γ.lon,γ.lat,δ¹³C[:, :,lvl]' .* d_NATL' ,
transform=TMI.cartopy.crs.PlateCarree()
)
ax.coastlines() #show coastlines
gl = ax.gridlines(draw_labels=true, dms=true, x_inline=false, y_inline=false)
gl.top_labels = false
gl.right_labels = false
ax.set_title("δ¹³C d_NATL")
fig.colorbar(cax, orientation = "horizontal")
fig
fig.savefig(plotsdir("δ¹³C_d_NATL.png"))


x = δ¹³C[:, :,lvl]' .* d_NATL'
δ¹³C_NA = mean(filter(!isnan,x))
println("δ¹³C_NA: ", δ¹³C_NA)
d_TROPPAC = getRegion(TMIversion,"d_TROPPAC")

fig = figure()
ax = fig.add_subplot(projection = TMI.cartopy.crs.PlateCarree())
ax.contourf(γ.lon,γ.lat,δ¹³C[:, :,lvl]' .* d_SUBANTPAC' ,
ax.contourf(γ.lon,γ.lat,δ¹³C[:, :,lvl]' .* d_TROPPAC' ,
transform=TMI.cartopy.crs.PlateCarree()
)
ax.coastlines() #show coastlines
gl = ax.gridlines(draw_labels=true, dms=true, x_inline=false, y_inline=false)
gl.top_labels = false
gl.right_labels = false
ax.set_title("δ¹³C d_TROPPAC")
x = δ¹³C[:, :,lvl]' .* d_SUBANTPAC'
δ¹³C_P = mean(filter(!isnan,x))
fig.savefig(plotsdir("δ¹³C_d_NATL.png"))
fig.colorbar(cax, orientation = "horizontal")
fig
fig.savefig(plotsdir("δ¹³C_d_TROPPAC.png"))


δ¹³C_AS = δ¹³C_P - δ¹³C_NA
# δ¹³C_AS = δ¹³C_P - δ¹³C_NA

println("δ¹³C_AS ", δ¹³C_AS)
# println("δ¹³C_AS ", δ¹³C_AS)

x = PO₄ᴿ[:, :,lvl]' .* d_NATL'
PO₄ᴿ_NA = mean(filter(!isnan,x))
R = -3 / 2.5
# x = PO₄ᴿ[:, :,lvl]' .* d_NATL'
# PO₄ᴿ_NA = mean(filter(!isnan,x))
# R = -3 / 2.5

pred_δ¹³C = ones(length(axes(δ¹³C,3)))
true_δ¹³C = similar(pred_δ¹³C)
for deep_lvl in axes(δ¹³C,3)
x = PO₄ᴿ[:, :,deep_lvl]' .* d_SUBANTPAC'
PO₄ᴿ_DP = mean(filter(!isnan,x))
RΔPO₄ᴿ = R * (PO₄ᴿ_DP - PO₄ᴿ_NA)
pred_δ¹³C[deep_lvl] = δ¹³C_NA + RΔPO₄ᴿ + δ¹³C_AS
x = δ¹³C[:, :,deep_lvl]' .* d_SUBANTPAC'
true_δ¹³C[deep_lvl] = mean(filter(!isnan,x))
end
# pred_δ¹³C = ones(length(axes(δ¹³C,3)))
# true_δ¹³C = similar(pred_δ¹³C)
# for level in size(δ¹³C,3)
# x = PO₄ᴿ[:, :,level]' .* d_SUBANTPAC'
# PO₄ᴿ_DP = mean(filter(!isnan,x))
# RΔPO₄ᴿ = R * (PO₄ᴿ_DP - PO₄ᴿ_NA)
# pred_δ¹³C[level] = δ¹³C_NA + RΔPO₄ᴿ + δ¹³C_AS
# x = δ¹³C[:, :,level]' .* d_SUBANTPAC'
# true_δ¹³C[level] = mean(filter(!isnan,x))
# end

figs, axs = subplots(1,3, figsize = (10, 5))
axs[1].plot(pred_δ¹³C, -depth)
axs[1].set_title("pred_δ¹³C")
axs[2].plot(true_δ¹³C, -depth)
axs[2].set_title("true_δ¹³C")
axs[3].plot(pred_δ¹³C .- true_δ¹³C, -depth)
axs[3].set_title("pred_δ¹³C .- true_δ¹³C")
# figs, axs = subplots(1,3, figsize = (10, 5))
# axs[1].plot(pred_δ¹³C, -depth)
# axs[1].set_title("pred_δ¹³C")
# axs[2].plot(true_δ¹³C, -depth)
# axs[2].set_title("true_δ¹³C")
# axs[3].plot(pred_δ¹³C .- true_δ¹³C, -depth)
# axs[3].set_title("pred_δ¹³C .- true_δ¹³C")

println("Compared to: ", mean(filter(!isnan,x)))
# println("Compared to: ", mean(filter(!isnan,x)))
148 changes: 46 additions & 102 deletions scripts/deltaC13vsPO4.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,16 @@
#=
% Example: Compare average δ¹³C and PO₄ values for different ocean basins %
% Example: Compare average δ¹³C and PO₄ values for different ocean regions %
=#
# import Pkg
# Pkg.activate("TMI.jl")#do this if using script interactively

using Revise
using TMI, PyPlot, PyCall, Colors
using TMI, PyPlot, PyCall, Statistics, DrWatson, PyCall, Colors

# pygui(true) #needed for Atom, not sure what it will do in other places
@pyimport matplotlib.colors as mcolors
@pyimport matplotlib.colors as plt_colors
mkpath(plotsdir())

# A, Alu, γ, inputfile = config(url,inputdir)
# ΔPO₄ = readtracer(inputfile,"qpo4")
TMIversion = "modern_90x45x33_G14_v2"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)
list = ("NATL","NPAC","ARC",
Expand All @@ -20,124 +22,66 @@ list = ("NATL","NPAC","ARC",
# surface ocean
depth = getVar(TMIversion, "depth")

#initialize a grid containing the water mass distriubtion
#for each ocean in "list"

which_water_mass = identifywatermass(TMIversion,Alu,γ, list)

# plot a section at 330 east longitude (i.e., 30 west)
lon_section = 330; # only works if exact
# wm_section = section(which_water_mass,lon_section,γ)
#initialize a grid containing the water mass sources
sources = identifywatermass(TMIversion,Alu,γ, list)

lims = 1:1.0:length(list)

cm = get_cmap(:tab20)
colorrange = (0:(length(list)-1)) ./ (length(list))

# oceans_mask = similar(wm_section[:,33:-1:1])

XC = γ.lon' .* ones(length(γ.lat))
YC = ones(length(γ.lon))' .* γ.lat

fig, ax = subplots(1,1)
counter = []
for ff in eachindex(XC)
val = which_water_mass[:, :, 20]'[ff]
if isnan(val)
# pass, dont plot anything
else
ax.scatter(XC[ff],YC[ff], color=cm(colorrange[Int(val)]))
if Int(val) ∉ counter
ax.scatter(XC[ff],YC[ff], color=cm(colorrange[Int(val)]), label = list[Int(val)])
counter = push!(counter, Int(val))
end
end
#oceans_mask[:, :, 1]', lims, cmap=cmap_seismic
end
ax.legend()
level = 20;
sources_in_layer = sources[:, :, level]
source_values = sort(filter(!isnan, unique(sources_in_layer)))
nsource = length(source_values)

#iterate through maskes for Ind, NPac, SPac, Na and Sa
# possibley include equatorial waters
# O13_values .* mask
# PO_values .* mask
# scatter(O13_values[:], PO_values[:])
#making a color bar
cmap = cm(colorrange[Int.(source_values)]);
my_cmap = plt_colors.ListedColormap(cmap)
bounds = source_values # define the bins and normalize
norm = plt_colors.BoundaryNorm(bounds, my_cmap.N)

#plot the distribution of source waters within a layer
fig, ax = subplots(1,1)
good_pts = (!isnan).(sources_in_layer) # filter for nonnan values
[ax.scatter(XC'[good_pts][1], YC'[good_pts][1], c = val, cmap = my_cmap, label = list[val], norm=norm, zorder = 0) for val in source_values] #plot the labels
fig.legend(ncols = 3)
ax.scatter(XC'[good_pts][:],YC'[good_pts][:], c = sources_in_layer[good_pts][:], cmap = my_cmap, norm=norm) #plot all points
fig.savefig(plotsdir("RegionalMasks.png"))
fig

#Plot PO vs C13 for various regions, shade the dots by their source region
fig, ax = subplots(1,1)
basin_list = ["","SUBANTIND", "NPAC",
"SUBANTPAC","NATL","SUBANTATL"]
# "TROPATL","TROPPAC","TROPIND")

δ¹³C = getVar(TMIversion,"δ¹³C")
PO₄ = getVar(TMIversion,"PO₄")
fig, axs = plt.subplots(2, 3, figsize = (13, 8), sharex = true, sharey = true)
nb = length(basin_list)
[a.set_xlabel("PO₄") for a in axs]
[a.set_ylabel("δ¹³C") for a in axs]

fig, axs = plt.subplots(2, 3, figsize = (13, 8))
# lvl = 2 # right below surface
# lvl = 14
lvl = 19
for i in 2:length(basin_list)
for i in 2:nb
ax = axs[i]
basin_mask = getRegion(TMIversion,"d_" * basin_list[i])
og_pts = sum(filter(!isnan, basin_mask))
println("Original mask has ",og_pts , " points")
trop_basin_mask = getRegion(TMIversion,
"d_" * "TROP" * basin_list[i][end-2:end])

if (basin_list[i][end-2:end] == "IND")
basin_list[i] = "IND"
hemisp = (γ.lat .> -1000)
elseif basin_list[i][1] == 'N'
hemisp = (γ.lat .> 0)
basin_list[i] = "N" * basin_list[i][end-2:end]
elseif basin_list[i][1] == 'S'
hemisp = (γ.lat .< 0)
basin_list[i] = "S" * basin_list[i][end-2:end]
end

ax.set_title(basin_list[i])

ax.set_xlim(0, 3.5)
ax.set_xticks(collect(0:0.5:3.5))
ax.set_ylim(-1, 3)
ax.set_yticks(collect(-1:0.5:3))

ax.set_xlabel("PO₄")
ax.set_ylabel("δ¹³C")

equatorpoints = findall(trop_basin_mask[:, hemisp] .== 1.0)
hem_basin_mask = @view basin_mask[:, hemisp]
hem_basin_mask[equatorpoints] .= 1.0
println("Changed mask has ",sum(filter(!isnan, basin_mask)) - og_pts, " addtl points")

mask_water_mass = which_water_mass[:, :, lvl:end] .* basin_mask
good_pts = findall(x->(x!= NaN) && ( x >0.0),mask_water_mass)
masked_sources = sources[:, :, level] .* basin_mask
masked_sources[iszero.(masked_sources)] .= NaN

good_pts = (!isnan).(mask_water_mass)
println("good points are: ", size(good_pts), "(should be geq to og + addtl pts)")
mask_water_mass = mask_water_mass[good_pts]
mask_δ¹³C = (δ¹³C[:, :, lvl:end] .* basin_mask)[good_pts]
mask_PO₄ = (PO₄[:, :, lvl:end] .* basin_mask)[good_pts]
counter = []

for ff in eachindex(mask_water_mass)
val = mask_water_mass[ff]
ax.scatter(mask_PO₄[ff],mask_δ¹³C[ff], color=cm(colorrange[Int(val)]), s = 1)
if Int(val) ∉ counter
ax.scatter(mask_δ¹³C[ff], mask_δ¹³C[ff], color=cm(colorrange[Int(val)]),
label = list[Int(val)],s = 5)
counter = push!(counter, Int(val))
end

end
ax.legend()

# good_idxs = findall( NaN .!= basin_mask .== 1.0)
# findall(x->(x!= NaN)&&(x>1),x)
# X = ones(length(δ¹³C[good_idxs]), 2)
# X[:, 2] .= δ¹³C[good_idxs][:]
# y = PO₄[good_idxs]
# β̂ = (X'*X)\X'*y
# ŷ=X*β̂
# ax.plot(δ¹³C[good_idxs][:], ŷ, color = "black")
masked_sources = masked_sources[good_pts]
mask_δ¹³C = (δ¹³C[:, :, level] .* basin_mask)[good_pts]
mask_PO₄ = (PO₄[:, :, level] .* basin_mask)[good_pts]
ax.scatter(mask_PO₄[:],mask_δ¹³C[:], c = masked_sources[:], cmap = my_cmap, norm=norm)
end
depth_range = string(depth[lvl]) * "m" * " to " * string(depth[end])
fig.suptitle(depth_range)

depth_range = string(depth[level])
fig.suptitle("PO₄vsδ¹³C for various regions at " * depth_range * "m")
fig.tight_layout()
fig.savefig("PO₄vsδ¹³C_" * depth_range * ".png")
close("all")
fig

0 comments on commit e44fc97

Please sign in to comment.