# Load World Ocean Atlas 2D fields for PO₄, Si(OH)₄, and NO₄
using WorldOceanAtlasTools
product_year = 2018
tracer = "p"
period = 0
resolution = "1°"
ds = WorldOceanAtlasTools.WOA_Dataset(product_year, tracer, period, resolution, verbose=false) ;
field = "an"
field3D, lat, lon, depth = WorldOceanAtlasTools.get_gridded_3D_field(ds, tracer, field) ;
map_2D = field3D[:,:,1] ;

# Use Cartopy
ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
clf()


# WOA colors RGB (sort-of interpolated from their figure for phosphate 2018)
using Colors, Plots
WOA_red    = RGB(0.957, 0.271, 0.188)
WOA_yellow = RGB(0.984 , 0.929, 0.298)
WOA_blue   = RGB(0.231, 0.631, 0.922)
# colrmap
C(g::ColorGradient, levels) = [(g[z].r, g[z].g, g[z].b) for z=range(0.0,1.0,length=length(levels))]
g = cgrad([WOA_blue, WOA_yellow, WOA_red])
levels = 0:0.2:2

# Define the projection
# Sydney: 33.8688° S, 151.2093° E
ax = subplot(projection=ccrs.NearsidePerspective(central_latitude=-33.8688, central_longitude=151.2093, satellite_height=100000000.0))

# Add black land
cfeature = pyimport("cartopy.feature")
ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#AAAAAA")      # gray land

# Make the data cyclic
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
map_cyc = hcat(map_2D, map_2D[:,1])
map_cyc = map(x -> ismissing(x) ? NaN : min(x, levels[end]), map_cyc)

# contour
p = PyPlot.contourf(lon_cyc, lat, map_cyc, levels=levels, transform=ccrs.PlateCarree(), zorder=-1, colors=C(g,levels))

#ax.gridlines(xlocs=collect(-180:3:180), ylocs=collect(-90:3:90), linewidth = 2, color="black", alpha = 1.0)

PyPlot.savefig("WOAT_logo.svg")